Hello R users,
I'm analyzing an experiment in a balanced incomplet block design (BIB). The
effect of blocks are assumed to be random, so I'm using nlme::lme for this.
I'm analysing another more complex experiments and I notice some diferences
from doBy::popMeans() compared multcomp::glht() and contrast::contrast().
In my example, glht() and contrast() were equal I suspect popMeans() are
doing something diferent. This a proprosital difference? Does popMeans
account the variance of random effects or something similar? The code below
is reproducible, for easy evaluation the principal results are appended.
Thank you.
# reading BIB data
da <- read.table("http://www.leg.ufpr.br/~walmes/data/pimentel_bib3.txt";,
+ header=TRUE, sep="\t",
colClasses=c("factor","factor",NA))
str(da)
## 'data.frame': 30 obs. of 3 variables:
## $ bloco : Factor w/ 10 levels "1","10","2","3",..: 1 1 1 3 3 3 4 4 4
5 ...
## $ variedade: Factor w/ 5 levels "1","2","3","4",..: 1 2 3 1 2 4 1 2 5 1
...
## $ y : int 35 28 27 30 20 22 28 16 18 36 ...
# fixed effect model
m0 <- lm(y~bloco+variedade, da)
summary(m0)
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 35.7556 1.0004 35.742 < 2e-16 ***
## bloco10 -0.3778 1.2729 -0.297 0.770447
## ...
## variedade2 -9.2000 0.9262 -9.933 3.01e-08 ***
## variedade3 -8.0667 0.9262 -8.710 1.81e-07 ***
## ...
require(contrast)
contr <- contrast(m0, type="average", list(bloco=levels(da$bloco),
variedade="1"))
contr
## Contrast S.E. Lower Upper t df Pr(>|t|)
## 1 32.76667 0.6438886 31.40168 34.13165 50.89 16 0
require(doBy)
pop <- popMeans(m0, effect="variedade")
pop
## beta0 Estimate Std.Error t.value DF Pr(>|t|) Lower Upper
variedade
## 1 0 32.76667 0.6438886 50.88872 16 0 31.40168
34.13165 1
summary(glht(m0, linfct=matrix(contr$X, nrow=1)))
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## 1 == 0 32.7667 0.6439 50.89 <2e-16 ***
# ok, results are the same!
# random effects model
require(nlme)
mm0 <- lme(y~variedade, random=~1|bloco, da)
summary(mm0)
## Fixed effects: y ~ variedade
## Value Std.Error DF t-value p-value
## (Intercept) 32.67807 1.5887794 16 20.568037 0
## variedade2 -9.08144 0.9236918 16 -9.831679 0
contr <- contrast(mm0, list(variedade="1"))
contr
## Contrast S.E. Lower Upper t df Pr(>|t|)
## 32.67807 1.588779 29.56412 35.79202 20.57 23 0
pop <- popMeans(mm0, effect="variedade")
pop
## beta0 Estimate Std.Error t.value DF Pr(>|t|) Lower Upper
variedade
## 1 0 30.50000 1.918043 15.90162 25 0 26.54972
34.45028 1
summary(glht(mm0, linfct=matrix(contr$X, nrow=1)))
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## 1 == 0 32.678 1.589 20.57 <2e-16 ***
# sdt error calculated by hand
sqrt(contr$X%*%vcov(mm0)%*%t(contr$X))
## 1
## 1 1.588779
==========================================================================
Walmes Marques Zeviani
LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W)
Departamento de Estatística - Universidade Federal do Paraná
fone: (+55) 41 3361 3573
VoIP: (3361 3600) 1053 1173
e-mail: wal...@ufpr.br
skype: walmeszeviani
twitter: @walmeszeviani
homepage: http://www.leg.ufpr.br/~walmes
linux user number: 531218
==========================================================================
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
