[R-sig-ME] contrasts in lmer models

[Frei Esther]

Hello Fernando

Thank you for your answer! I am not sure whether I understand correctly what you suggest. I have read the book of Faraway and have already calculated LRTs in order to test the significance of the fixed and random factors. Following that, now, I would like to test contrasts between the individual levels of the fixed factors (e.g. my factor Alt has three levels 600m, 1200m, and 1800m) and also because there are significant interactions in my model, I would like to estimate and test  various contrasts between different combinations of the factor levels. 

Best regards,
Esther


-----Original Message-----
From: fernando barbero [mailto:fbarbero at bariloche.inta.gov.ar] 
Sent: Donnerstag, 7. Juni 2012 15:04
To: Frei Esther; mailman, r-sig-mixed-models
Subject: RE: [R-sig-ME] contrasts in lmer models

Why don t you test those contrasts using the LRT aproach? A full model against a reduced one with all the effects except the one you are trying to test. The book "Extending the linear model with R" (by Julian Faraway) surely will help you Best regards Fernando

-----Mensaje original-----
De: r-sig-mixed-models-bounces at r-project.org
[mailto:r-sig-mixed-models-bounces at r-project.org] En nombre de Frei Esther Enviado el: miércoles, 06 de junio de 2012 03:40 p.m.
Para: mailman, r-sig-mixed-models
Asunto: [R-sig-ME] contrasts in lmer models

Dear R experts,

I am fitting a linear mixed model with lmer. Some of the a-priori-contrasts I want to test involve interaction terms. I ran into serious trouble when I tried to test these contrasts:

1. For testing interactions in the fixed factors, I tried using
fit.contrast() from the gmodels-package. In principle this works, but the results I get are very suspicious. They vary quite a lot if I execute fit.contrasts several times. I have increased the number of samples from the default  1000 up to 200,000. However, the results don't seem to converge. It seems also strange that sometimes p-values are exactly 0. Does it mean that the contrasts are highly significant (i.e. smaller than the number of decimals displayed in the result)? Or is there something wrong in the way I use the function? Maybe I simply have to run more samples with my complex model?

> require("gmodels")
> ffin <- lmer(y ~ CoV + Alt + Orig + Soil + (1 | Bed) + (1 | Site) + (1 
> |
Pair),  data = TETri,
+             na.action = na.exclude, REML=TRUE)
> cmat <- rbind("First Contrast (600m vs.1200m)"= c(1,0,-1),"Second 
> Contrast
(1200m vs.1800m)"=
 +                         c(1,-1,0))

> fit.contrast(ffin,"Alt",cmat,conf.int=0.95,sim.mer=TRUE,n.sim=5000)
 
Estimate Std. Error   p value  Lower.CI    Upper.CI
AltFirst Contrast (600m vs.1200m)   -4.879821   0.979278 0.0000000 -6.317718
-3.86426640
AltSecond Contrast (1200m vs.1800m) -2.867170   1.796602 0.3333333 -4.593912
0.03296012


2. Some of the contrasts that I want to test involve interaction terms. As I understand these cannot directly be approached with fix.contrasts(). Instead I tried to recode the contrast terms into a new combined variable using factorA.factorB = paste(factorA,factorB,sep="."). This only works when I also refit the model with only the new combined factor as fixed effect but results in suspicious p-values of '0'. As I understand, this is not the same model anymore and so it is not very useful.

> TETri$O.A.Soil  <- as.factor(paste(TETri$Orig, TETri$Alt, TETri$Soil, 
> sep=
levels(TETri$O.A.Soil) 
> cmat_1 <- rbind("HA Orig12 (O12.T12.SOL vs. O18.T18.GRA)"=
c(0,1,0,0,0,0,0,0,-1,0,0,0),"LF Orig12 + (O12.T12.SOL vs. O18.T12.SOL)"= c(0,1,0,0,0,0,0,-1,0,0,0,0))
> f1c <- lmer(y ~ CoV + O.A.Soil + (1 | Bed) + (1 | Site) + (1 | Pair),
data = TETri,
+            na.action = na.exclude, REML=TRUE)
> fit.contrast(f1c,"O.A.Soil",cmat_1,conf.int=0.95,sim.mer=TRUE,n.sim=10
> 00)
 
Estimate Std. Error p value   Lower.CI   Upper.CI
O.A.SoilHA Orig12 (O12.T12.SOL vs. O18.T18.GRA) -3.958218   2.110968       0
-6.6922765 -0.7359227
O.A.SoilLF Orig12 (O12.T12.SOL vs. O18.T12.SOL)  2.714403   1.387434       0
0.5868685  5.3193961


3. Using the function estimable() from gmodels, I have trouble specifying the contrast matrix correctly. The help entry does not really give a good example how to indicate the factor to which I want to refer to. Could anyone give me a hint (example) how to define the contrasts correctly, please?
> f0 <- lmer(y ~ CoV + Alt * Orig * Soil + (1 | Bed) + (1 | Site) + (1 |
Pair),  data = TETri,
+             na.action = na.exclude, REML=TRUE)
> cmat <- rbind('Alt600m vs.Alt1200m'= c(1,0,-1),'Alt1200m vs.Alt1800m'=
c(1,-1,0))
> estimable(f0,cmat)
Error in FUN(newX[, i], ...) : 'param' has no names and does not match number of coefficients of model. Unable to construct coefficient vector


4. My last question is a general question not specifically related to R: If one of the fixed terms in my mixed model is not significant overall, is it possible that one of the contrasts between levels of that term is still significant?

Thanks a lot for your help!
Esther

_______________________________________________
R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
-----
Se certificó que el correo no contiene virus.
Comprobada por AVG - www.avg.es
Versión: 2012.0.2177 / Base de datos de virus: 2433/5053 - Fecha de la
versión: 06/06/2012