Solved – Controlling nuisance covariate in lmer

ancovamixed modelr

I want to control for a nuisance covariate in a linear model. Since the covariate interacts significantly with one of the fixed factors, the homogeneity of regression slopes assumption is violated for an ANCOVA approach. My understanding is that this can be handled in a multi-level, or mixed effects, model.

My data:

> str(seeds)
'data.frame':   186 obs. of  4 variables:
 $ fixed.A : Factor w/ 31 levels "A1","A2","A3",..: 7 7 7 7 7 7 10 10 10 10 ...
     $ fixed.B : Factor w/ 2 levels "Y","N": 2 2 2 1 1 1 2 2 2 1 ...
 $ cov     : num  10.3 10.5 11 12.8 12.9 ...
     $ response: num  10.8 11 11.1 14.7 15.3 ...

The covariate cov has a significant interaction with fixed.A. To fit a random intercept and slope for cov conditioned on fixed.A, it seems an approach using lmer (in lme4) might be:

> lmer(response ~ fixed.A*fixed.B + (1 + cov | fixed.A), data = seeds)

I'm aware that this is not the usual approach in lmer since grouping factors tend to be random and therefore don't also appear as fixed factors in the model formula. However, since I'm interested in the interaction between my two fixed factors, I don't see how else to proceed. Any help is greatly appreciated.

Best Answer

Your model,

response ~ fixed.A*fixed.B + (1 + cov | fixed.A)

translates to

$$ \begin{split} y_{ijk} & \sim \textrm{N}(\mu_{ijk},\sigma^2_r) \\ \mu_{ijk} & = \beta_{ij} + b_{i1} + b_{i2} x_{ijk} \\ (b_{i1}, b_{i2}) & \sim \textrm{MVN}(\boldsymbol 0,\Sigma(\boldsymbol \theta)) \end{split} $$

The problem is that some of the $\beta_{ij}$ parameters and the $b_{i1}$ parameters are jointly unidentifiable.

If you instead use

response~fixed.A*fixed.B + cov + (0 + cov|fixed.A)

you eliminate the random intercept and get

$$ \begin{split} y_{ijk} & \sim \textrm{N}(\mu_{ijk},\sigma^2_r) \\ \mu_{ijk} & = \beta_{ij} + \beta_x x_{ijk} + b_{i2} x_{ijk} \\ b_{i2} & \sim \textrm{N}(0,\theta^2) \end{split} $$

(For a sensible model you need to include an overall fixed effect of $x$ (cov), because the random effects are zero-centered.)

Related Solutions

Mixed Model – Getting Fixed-Effect Only Predictions on New Data in R

Maybe this is cheating, but

fixedformula <- as.formula(lme4.0:::nobars(formula(m))[-2])
model.matrix(fixedformula,newdata=cake2)

note:

I am using lme4.0 here, which is the r-forge version of "old" (CRAN) lme4: you can substitute lme4 for lme4.0 in the code above
the new (r-forge/development) version of lme4 has a predict method: in that case
```
predict(m,re.form=NA,newdata=cake2)
```

works fine (re.form=NA sets all random effects to zero, equivalent to level=0 in the old predict.lme)

Solved – Mixed effects model hypothesis testing

As I understand it, this issue is not at all specific to mixed models, but is a variant of "the difference between significant and non-significant is not significant" (see also other blog posts by Gelman); see also Pockock et al.2002. Consider the following statements:

"within C1, we can reject the null hypothesis of no $A \times B$ interaction; the effect of A on responses differs significantly across levels of B (and vice versa)"
"within C2, we cannot reject the null hypothesis; the effect of A on responses does not differ significantly ..."
"we cannot reject the hypothesis that the $A \times B$ interaction is the same in both levels of C"

All three are true. As @Wolfgang points out below, the test of the three-way interaction is the formal test that there is a difference between the strength of the interaction across different levels.

The corresponding picture:

enter image description here

Best Answer

Related Solutions

Mixed Model – Getting Fixed-Effect Only Predictions on New Data in R

Solved – Mixed effects model hypothesis testing

Related Question