[R-sig-ME] [R-sig-eco] question on using ghlt (multcomp) for mixed models with lme (nlme)
Kay Cecil Cichini
Kay.Cichini at uibk.ac.at
Tue Jun 14 12:20:30 CEST 2011
Zitat von Hem Nalini Morzaria Luna <hmorzarialuna at gmail.com>:
> Hi:
>
> I am comparing biomass response in a an experiment with three
factors
> (N-two levels, W-three levels, and S-two levels) using lme. Unit is
> just specified as a block effect. (Data example below)
>
> This is my model
>
>> w.lme <- lme(fixed=log(Ry.above) ~ W * S * N, data= Mono, random=~1
| Unit/W/N)
>
> I'd like to do multiple comparisons.
> In this case, interaction terms of W:Species is slightly significant
> (p 0.07) so if I want to only test for W within each Species.
>
> I tried the following code to get a contrast matrix for those
factors
>
>> library(contrast)
>> contrast(w.lme, list(W=c('High','Medium','Low'), Species = 'Tc'),
list(W=c('High','Medium','Low'), Group = 'Sv'), type='average')
>
> But I get
>
>> Error in gendata.default(fit = list(modelStruct = list(reStruct =
list( :
> not enough factors
>
a blind guess: maybe W, S, N are vectors and not factors, in that
case W<-as.factor(W); etc. may solve your problem..
> I was also trying to understand what the correct way to specify
> multiple comparions for the whole model is:
>
maybe you can use this as a starter:
http://thebiobucket.blogspot.com/2011/06/glmm-with-custom-multiple-comparisons.html
though, i never checked if the method used in this comparisons is
suitable for
(generalised) mixed models.
i'd be very glad about any useful comments on this topic!
> One factor is straightforward
>> comp.Species <- glht(w.lme, linfct=mcp(Species="Tukey"))
>
> To specify all comparisons I know mcp will not take interactions but
> the linfct argument can use a matrix of comparisons but I am unclear
> what is the correct approach to set up that matrix.
>
> I tried:
>
> 1) matrix from fixed effects
>
>> K <- diag(length(fixef(w.lme)))
>> glht(w.lme, linfct=K)
>
> 2) matrix using contrasts
>
>> matrix.con = contrasts(Mono$S:Mono$N:Mono$W)
> #matrix dimension is [11,12]
> #so I added one more column, results in this matrix being equal to
K(above)
>> missing.column = matrix(0,12,1)
>> missing.column[1,1] = 1
> con.matrix = cbind(missing.column,matrix.con)
>> glht(w.lmer, linfct = con.matrix)
>
> 1 & 2 give the same result from glht. Is it correct to assume the
rows
> in K are the same as in the contrast matrix?
>
> 3) Model matrix using model.matrix
>> lme.matrix <- model.matrix(w.lme, data=Mono)
> glht(w.lmer, linfct = lme.matrix)
>
> In 3) lme.matrix dimension is [36,12] and the results of glht repeat
> three times each. I am having a hard time understanding why this is
> and what the rows and columns in the matrix represent. Colnames are
an
> incomplete set of contrasts.
>
> What is the correct approach to specify all comparisons?
>
>
> Any help is appreciated. Thank you!
>
> Hem Nalini Morzaria Luna
>
> Data (Mono)
> N W Unit Species Ry.above
> 1 High High 1 Sv 0.6156346
> 2 High High 2 Sv 0.6235052
> 3 High High 3 Sv 0.6985453
> 4 Low High 4 Sv 1.7452381
> 5 Low High 5 Sv 0.8424228
> 6 Low High 6 Sv 0.6643357
> 7 High Medium 7 Sv 0.8251745
> 8 High Medium 8 Sv 0.7424867
> 9 High Medium 9 Sv 0.9362663
> 10 Low Medium 10 Sv 1.0038536
> 11 Low Medium 11 Sv 0.7767575
> 12 Low Medium 12 Sv 1.5007492
> 13 High Low 13 Sv 0.9536121
> 14 High Low 14 Sv 0.9799515
> 15 High Low 15 Sv 0.5596921
> 16 Low Low 16 Sv 0.8706301
> 17 Low Low 17 Sv 1.0032415
> 18 Low Low 18 Sv 1.3620130
> 19 High High 1 Tc 0.6514929
> 20 High High 2 Tc 0.4309266
> 21 High High 3 Tc 0.2800695
> 22 Low High 4 Tc 0.7883481
> 23 Low High 5 Tc 0.4541833
> 24 Low High 6 Tc 0.6007958
> 25 High Medium 7 Tc 0.6987778
> 26 High Medium 8 Tc 0.2240063
> 27 High Medium 9 Tc 0.1477910
> 28 Low Medium 10 Tc 0.2017751
> 29 Low Medium 11 Tc 0.3621421
> 30 Low Medium 12 Tc 0.2590503
> 31 High Low 13 Tc 0.4265925
> 32 High Low 14 Tc 0.5525204
> 33 High Low 15 Tc 0.4196785
> 34 Low Low 16 Tc 0.4507284
> 35 Low Low 17 Tc 1.1027065
> 36 Low Low 18 Tc 0.3676259
>
> _______________________________________________
> R-sig-ecology mailing list
> R-sig-ecology at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
>
>
More information about the R-sig-mixed-models
mailing list