[R] multcomp::glht() doesn't work for an incomplete factorial using aov()?

[Joshua Wiley]

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/