I've fit a linear mixed effects model to some data in R with afex::mixed. I'm interested in the fixed effect and have no expectations for the random effects structure; I just want the most parsimonious model specification.
My data consists of a numeric outcome y, a two-level categorical predictor a (up vs down), and participant IDs p. I want to know the fixed effect if a on y, while accounting for random effects by p.
> str(DF)
'data.frame': 1521 obs. of 3 variables:
$ p: Factor w/ 100 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 2 2 ...
$ a: Factor w/ 2 levels "Down","Up": 2 2 2 2 2 2 2 2 2 2 ...
$ y: num 12 0 13 0 0 10 5 0 0 5 ...
First I fit the maximally complex random effects structure:
> m1 <- mixed(y ~ a + (a|p), # maximal random effects structure
+ data = DF,
+ expand_re = TRUE,
+ method = "S",
+ return = "merMod")
Contrasts set to contr.sum for the following variables: a, p
boundary (singular) fit: see ?isSingular
> summary(m1)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: y ~ a + (1 + re1.a1 | p)
Data: data
REML criterion at convergence: 11625.7
Scaled residuals:
Min 1Q Median 3Q Max
-4.1982 -0.4479 -0.1815 0.2700 4.9898
Random effects:
Groups Name Variance Std.Dev. Corr
p (Intercept) 11.359 3.370
re1.a1 8.653 2.942 1.00
Residual 112.938 10.627
Number of obs: 1521, groups: p, 100
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 6.8167 0.4414 106.1053 15.44 <2e-16 ***
a1 4.6929 0.4093 112.5658 11.46 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr)
a1 0.607
convergence code: 0
boundary (singular) fit: see ?isSingular
It's singular. So I reduce the random effects structure by suppressing the correlation terms:
> m2 <- mixed(y ~ a + (a||p), # reduced complexity spec
+ data = DF,
+ expand_re = TRUE,
+ method = "S",
+ return = "merMod")
Contrasts set to contr.sum for the following variables: a, p
> summary(m2)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: y ~ a + (1 + re1.a1 || p)
Data: data
REML criterion at convergence: 11693.8
Scaled residuals:
Min 1Q Median 3Q Max
-4.0236 -0.5564 -0.1455 0.2483 5.0291
Random effects:
Groups Name Variance Std.Dev.
p (Intercept) 7.726 2.78
p.1 re1.a1 4.841 2.20
Residual 118.777 10.90
Number of obs: 1521, groups: p, 100
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 6.8375 0.4039 100.1387 16.93 <2e-16 ***
a1 4.6081 0.3660 99.1904 12.59 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr)
a1 0.028
No more singular fit warning. Looks good. But I'd like treatment contrasts to more easily interpret the fixed effect:
> m2 <- mixed(y ~ a + (a||p), # reduced complexity spec
+ data = DF,
+ check_contrasts = FALSE, # use treatment contrasts
+ expand_re = TRUE,
+ method = "S",
+ return = "merMod")
boundary (singular) fit: see ?isSingular
> summary(m2)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: y ~ a + (1 + re1.aUp || p)
Data: data
REML criterion at convergence: 11705.4
Scaled residuals:
Min 1Q Median 3Q Max
-3.8373 -0.6266 -0.1392 0.2322 5.0521
Random effects:
Groups Name Variance Std.Dev.
p (Intercept) 6.561 2.561
p.1 re1.aUp 0.000 0.000
Residual 124.179 11.144
Number of obs: 1521, groups: p, 100
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 11.3858 0.4955 242.6631 22.98 <2e-16 ***
aUp -9.1344 0.5842 1514.8102 -15.63 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr)
aUp -0.619
convergence code: 0
boundary (singular) fit: see ?isSingular
Now the singular fit is back... What's happening here? Why does changing the contrasts determine whether the model is singular or not?
