R help - Nov 2012 - Diference in results from doBy::popMeans, multcomp::glht and contrast::contrast for a lme model

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: walmes@ufpr.br
skype: walmeszeviani
twitter: @walmeszeviani
homepage: http://www.leg.ufpr.br/~walmes
linux user number: 531218
=========================================================================
	[[alternative HTML version deleted]]