In this note I show how to suppress estimation of covariance random effects when estimating random effects for a categorical predictor.
The following data come from Ch 5 of Linear Mixed Models by West et al.
An experiment that examined nucleotide activation in 5 different rats. Each rat was subjected to two different treatments (basal vs carbachol) and then nucleotide activation was measured in 3 different brain regions. Therefore each rat has 6 different measures.
library(lme4)library(afex)library(nlme)ratbrain <-read.csv("ratbrain.csv")head(ratbrain, n =8)
Model activate as a function of region and treatment (and their interaction) with random effects on the intercept and region. Notice the singular fit message, This is due to the ill-conditioned random-effects covariance matrix, which suggests the model may be over-parameterized. Specifically we may not need to model covariance between the levels of the region random effects. The summary shows them all as -1 or 1.
m1 <-lmer(activate ~ region * treatment + (region | animal), data = ratbrain)
boundary (singular) fit: see help('isSingular')
summary(m1, corr = F)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: activate ~ region * treatment + (region | animal)
Data: ratbrain
REML criterion at convergence: 275.2
Scaled residuals:
Min 1Q Median 3Q Max
-1.38732 -0.58863 0.08712 0.41762 1.56355
Random effects:
Groups Name Variance Std.Dev. Corr
animal (Intercept) 5287.56 72.716
regionLS 32.38 5.690 -1.00
regionVDB 12.26 3.502 -1.00 1.00
Residual 2443.21 49.429
Number of obs: 30, groups: animal, 5
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 428.506 39.321 5.694 10.898 5.00e-05 ***
regionLS -190.762 31.365 19.630 -6.082 6.56e-06 ***
regionVDB -216.212 31.301 19.855 -6.908 1.08e-06 ***
treatmentcarbachol 98.204 31.262 20.000 3.141 0.00514 **
regionLS:treatmentcarbachol 99.322 44.210 20.000 2.247 0.03612 *
regionVDB:treatmentcarbachol 261.822 44.210 20.000 5.922 8.60e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
How can we fit a model without the covariance between random effects? Unfortunately, we can’t do that in {lme4}. However we can do so with the {afex} or {nlme} packages.
The {afex} package provides the lmer_alt() function that allows the use of two bars between region and animal (region || animal) to suppress estimation of random effect covariances (i.e., they are held constant at 0). However we still get warning about singularities.
m2 <-lmer_alt(activate ~ region * treatment + (region || animal), data = ratbrain)
boundary (singular) fit: see help('isSingular')
summary(m2, corr = F)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: activate ~ region * treatment + (1 + re1.regionLS + re1.regionVDB ||
animal)
Data: data
REML criterion at convergence: 275.3
Scaled residuals:
Min 1Q Median 3Q Max
-1.4643 -0.5762 0.0802 0.4072 1.5509
Random effects:
Groups Name Variance Std.Dev.
animal (Intercept) 4850 69.64
animal.1 re1.regionLS 0 0.00
animal.2 re1.regionVDB 0 0.00
Residual 2450 49.50
Number of obs: 30, groups: animal, 5
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 428.506 38.210 7.484 11.214 6.04e-06 ***
regionLS -190.762 31.307 20.000 -6.093 5.91e-06 ***
regionVDB -216.212 31.307 20.000 -6.906 1.04e-06 ***
treatmentcarbachol 98.204 31.307 20.000 3.137 0.00519 **
regionLS:treatmentcarbachol 99.322 44.275 20.000 2.243 0.03636 *
regionVDB:treatmentcarbachol 261.822 44.275 20.000 5.914 8.76e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
In the {nlme} package we can use the pdDiag() to specify we want to fit a diagonal random effects matrix with covariances fixed at 0. Notice this requires a named list object, where the name is the grouping factor. This converges with no warning messages. This may be due to lme() using a different optimizer. It uses “BFGS” while lmer() uses “nloptwrap”.
library(nlme)m3 <-lme(activate ~ region * treatment, random =list(animal =pdDiag(~ region)), data = ratbrain)summary(m3)
Linear mixed-effects model fit by REML
Data: ratbrain
AIC BIC logLik
295.2822 307.0627 -137.6411
Random effects:
Formula: ~region | animal
Structure: Diagonal
(Intercept) regionLS regionVDB Residual
StdDev: 69.64057 0.002818092 0.003204284 49.50045
Fixed effects: activate ~ region * treatment
Value Std.Error DF t-value p-value
(Intercept) 428.506 38.21022 20 11.214434 0.0000
regionLS -190.762 31.30684 20 -6.093302 0.0000
regionVDB -216.212 31.30684 20 -6.906223 0.0000
treatmentcarbachol 98.204 31.30684 20 3.136823 0.0052
regionLS:treatmentcarbachol 99.322 44.27455 20 2.243320 0.0364
regionVDB:treatmentcarbachol 261.822 44.27455 20 5.913600 0.0000
Correlation:
(Intr) regnLS rgnVDB trtmnt rgnLS:
regionLS -0.410
regionVDB -0.410 0.500
treatmentcarbachol -0.410 0.500 0.500
regionLS:treatmentcarbachol 0.290 -0.707 -0.354 -0.707
regionVDB:treatmentcarbachol 0.290 -0.354 -0.707 -0.707 0.500
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-1.46425392 -0.57621824 0.08020305 0.40718429 1.55091561
Number of Observations: 30
Number of Groups: 5
References
West, B., Welch, K., & Galecki, A. (2015) Linear Mixed Models. Chapman and Hall/CRC.
