R help - Aug 2011 - multcomp::glht() doesn't work for an incomplete factorial using aov()?

Hi R users,

I sent a message yesterday about NA in model estimates (
http://r.789695.n4.nabble.com/How-set-lm-to-don-t-return-NA-in-summary-td3722587.html).
If I use aov() instead of lm() I get no NA in model estimates and I use
gmodels::estimable() without problems. Ok!
Now I'm performing a lot of contrasts and I need correcting for
multiplicity. So, I can use multcomp::glht() for this. However, glht()
return an error message that is not compatible with my expectations. Someone
know I or has a suggestion for? Below some reproducible code.

# toy data
adi <- expand.grid(cult=gl(1,3,la=LETTERS[1]), fert=101)
fat <- expand.grid(cult=gl(2,3,la=LETTERS[2:3]), fert=seq(50,150,50))
da <- rbind(adi, fat)
da$y <- da$fert+rnorm(nrow(da),0,10)

# plot
require(lattice)
xyplot(y~fert|cult, da)

# adjust complete factorial
m0 <- aov(y~cult*fert, subset(da, cult!="A"))
summary(m0)
coef(m0)

# adjust incomplete factorial
m1 <- aov(y~cult*fert, da)
summary(m1)
coef(m1)

require(multcomp)
glht(m0, linfct=matrix(c(1,1,10,0), nrow=1))   # work
glht(m1, linfct=matrix(c(1,1,0,10,0), nrow=1)) # don't work

Thank you.
Walmes.

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

Hi Walmes,

The problem goes back to the singularity caused by the incomplete
design.  The summary() and coef() methods for aov class objects hide
it rather than printing NA as with lm(), but look at just:

m1

it says at the end "1 out of 6 effects not estimable".  Now when you
use glht() you pass a 1 x 5 matrix, but there are really *6* effects.
glht drops inestimable parameters, and it knows the 6th should be
dropped, so it is trying to subset your linear hypothesis matrix to
only those effects that are estimable.  The code basically does:
> c(rep(TRUE, 5), FALSE) # TRUE being estimable effects
[1]  TRUE  TRUE  TRUE  TRUE  TRUE FALSE
> matrix(c(1,1,0,10,0), nrow = 1) # your hypothesis matrix
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    1    0   10    0
> matrix(c(1,1,0,10,0), nrow = 1)[, c(rep(TRUE, 5), FALSE)] # how glht tries
to subset it
Error: (subscript) logical subscript too long # the error you get

The solution is relatively straight forward, just include a value in
your H matrix for the effect that cannot be estimated.  For example:

glht(m1, linfct=matrix(c(1,1,0,10,0,0), nrow=1))

Cheers,

Josh

On Sat, Aug 6, 2011 at 7:48 AM, Walmes Zeviani <walmeszeviani at
gmail.com> wrote:
> Hi R users,
>
> I sent a message yesterday about NA in model estimates (
>
http://r.789695.n4.nabble.com/How-set-lm-to-don-t-return-NA-in-summary-td3722587.html).
> If I use aov() instead of lm() I get no NA in model estimates and I use
> gmodels::estimable() without problems. Ok!
> Now I'm performing a lot of contrasts and I need correcting for
> multiplicity. So, I can use multcomp::glht() for this. However, glht()
> return an error message that is not compatible with my expectations.
Someone
> know I or has a suggestion for? Below some reproducible code.
>
> # toy data
> adi <- expand.grid(cult=gl(1,3,la=LETTERS[1]), fert=101)
> fat <- expand.grid(cult=gl(2,3,la=LETTERS[2:3]), fert=seq(50,150,50))
> da <- rbind(adi, fat)
> da$y <- da$fert+rnorm(nrow(da),0,10)
>
> # plot
> require(lattice)
> xyplot(y~fert|cult, da)
>
> # adjust complete factorial
> m0 <- aov(y~cult*fert, subset(da, cult!="A"))
> summary(m0)
> coef(m0)
>
> # adjust incomplete factorial
> m1 <- aov(y~cult*fert, da)
> summary(m1)
> coef(m1)
>
> require(multcomp)
> glht(m0, linfct=matrix(c(1,1,10,0), nrow=1)) ? # work
> glht(m1, linfct=matrix(c(1,1,0,10,0), nrow=1)) # don't work
>
> Thank you.
> Walmes.
>
>
=========================================================================>
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 at ufpr.br
> twitter: @walmeszeviani
> homepage: http://www.leg.ufpr.br/~walmes
> linux user number: 531218
>
=========================================================================>
> ? ? ? ?[[alternative HTML version deleted]]
> ______________________________________________
> R-help at 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.
-- 
Joshua Wiley
Ph.D. Student, Health Psychology
Programmer Analyst II, ATS Statistical Consulting Group
University of California, Los Angeles
https://joshuawiley.com/