[R-sig-ME] contrasts in lmer models

[Frei Esther]

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=1000)
                                                                                            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