[R] Interaction term in lmer

[Douglas Bates]

On 6/1/07, emine özgür Bayman <eozgur10 at hotmail.com> wrote:
> Dear R users,
> I'm pretty new on using lmer package. My response is binary and I have fixed
> treatment effect (2 treatments) and random center effect (7 centers). I want
> to test the effect of treatment by fitting 2 models:
>
> Model 1: center effect (random) only
> Model 2: trt (fixed) + center (random) + trt*center interaction.
>
> Then, I want to compare these 2 models with Likelihood Ratio Test. Here are
> my lmer codes that I don't feel comfortable about their correctness.
>
> model1 <- try(lmer(cbind( yvect, nvect-yvect) ~ 1 + (1 | center),
>   family = binomial, niter = 25, method = "Laplace", control = list(usePQL =
> FALSE) ))
>
> model2 <- try(lmer(cbind( yvect, nvect-yvect) ~ trt*center  + ( 1 | center)
> ,
>   family = binomial, niter = 25, method = "Laplace", control = list(usePQL =
> FALSE) ))

As you have seen, that model includes center as both a fixed effect
and as a grouping factor for a random effect.  You don't want to do
that.

I think the model you want is either

cbind(yvect, nvect-yvect) ~ trt + (trt|center)

which allows for correlation between the random effect for the
intercept and the random effect for treatment within center, or

cbind(yvect, nvect-yvect) ~ trt + (1|center/trt)

which is equivalent to

cbind(yvect, nvect-yvect) ~ trt + (1|center) + (1|center:trt)

This model has a random effect for the center and a random effect for
each center:trt combination but the variance of the latter random
effects are assumed to be constant.  Thus there are a total of three
variance components (random effects for center, for center:trt
combinations and for the per-observation noise term) in this model.
The previous model will have 1 + (ntrt * (ntrt + 1))/2 variances and
covariances to be estimated (where ntrt is the number of levels of the
treatment factor).

I hope this helps.

> (I have attached outputs below)
> What I don't understand is; I thought in model2 I have defined center effect
> as a random effect. Then how come I got center effects and trt*center
> interactions under the fixed effects list on the output? Probably I didn't
> understand how to set up these models in lmer. Could anyone help me about
> this?
>
> Thanks a lot for your help...
>
> Emine
>
> model1 <- try(lmer(cbind( yvect, nvect-yvect) ~ 1 + (1 | center),
>   family = binomial, niter = 25, method = "Laplace", control = list(usePQL =
> FALSE) ))
>
> >summary(model1)
> Generalized linear mixed model fit using Laplace
> Formula: cbind(yvect, nvect - yvect) ~ 1 + (1 | center)
> Family: binomial(logit link)
>      AIC      BIC    logLik deviance
> 236.817 238.0951 -116.4085  232.817
> Random effects:
> Groups Name        Variance Std.Dev.
> center (Intercept) 0.088127 0.29686
> number of obs: 14, groups: center, 7
>
> Estimated scale (compare to 1)  0.2672612
>
> Fixed effects:
>             Estimate Std. Error z value Pr(>|z|)
> (Intercept) -0.32084    0.14709 -2.1812  0.02917 *
> ---
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>
>
> ##################
> model2 <- try(lmer(cbind( yvect, nvect-yvect) ~ trt*center  + ( 1 | center)
> ,
>   family = binomial, niter = 25, method = "Laplace", control = list(usePQL =
> FALSE) ))
>
> >summary(model2)
> Generalized linear mixed model fit using Laplace
> Formula: cbind(yvect, nvect - yvect) ~ trt * center + (1 | center)
> Family: binomial(logit link)
> AIC      BIC        logLik     deviance
>   30 39.58586 -1.547024e-07 3.094048e-07
> Random effects:
> Groups Name        Variance Std.Dev.
> center (Intercept) 5e-10    2.2361e-05
> number of obs: 14, groups: center, 7
>
> Estimated scale (compare to 1)  0.2672612
>
> Fixed effects:
>               Estimate Std. Error  z value Pr(>|z|)
> (Intercept)  -1.060869   0.065372 -16.2282  < 2e-16 ***
> trt2          1.118029   0.086842  12.8743  < 2e-16 ***
> center2      -0.325428   0.504256  -0.6454  0.51869
> center3      -0.325440   0.504258  -0.6454  0.51868
> center4       0.655407   0.413449   1.5852  0.11292
> center5      -0.325433   0.504256  -0.6454  0.51869
> center6      -0.325421   0.504255  -0.6454  0.51870
> center7       0.655422   0.413448   1.5853  0.11291
> trt2:center2  0.673737   0.651313   1.0344  0.30094
> trt2:center3 -0.137183   0.651314  -0.2106  0.83318
> trt2:center4 -0.307083   0.583845  -0.5260  0.59891
> trt2:center5 -0.137203   0.651314  -0.2107  0.83316
> trt2:center6  1.654555   0.712419   2.3224  0.02021 *
> trt2:center7  0.673705   0.651311   1.0344  0.30096
> ---
>
> _________________________________________________________________
>
> Office Live http://clk.atdmt.com/MRT/go/aub0540003042mrt/direct/01/
>
> ______________________________________________
> R-help at stat.math.ethz.ch 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.
>