Factor Contrasts and Factor Levels in Linear Mixed Models
Reinhold Kliegl
2018-05-27 13:10:46
- 1 Overview / Summary
- 2 Setup
- 3 Means and correlations of Machines data
- 4 Maximal LMMs
- 5 Zero-correlation parameter LMMs
- 6 Interaction between fixed and random factor (intLMMs)
- 7 Minimal LMMs varying only in intercept (minLMMs)
- 8 Is intLMM nested under zcpLMM or not?
- 9 Parsimonious LMMs: Treatment contrasts
- 10 Parsimonious LMMs: Sum contrasts
- 11 Bates (2008) on interaction LMM (intLMM)
- 12 Appendix
1 Overview / Summary
Linear mixed models (LMMs) vary in the complexity of their random-effect structure. Model comparisons with likelihood ratio tests (LRTs) assume that the LMMs to be compared are in a hierarchical order, that is simpler LMMs are nested under more complex ones. This script compares LMMs varying in complexity with respect to a basic feature of the R model formula, that is whether 1 + factor or 0 + factor is used in the fixed-effect (FE) or random-effect (RE) part of the model.
The script is also a follow-up to an RPub on reduction of LMM complexity with double-bar syntax. The update was motivated by an exchange in a thread at R-Sig-ME. The thread is about the question whether an LMM with an interaction between a fixed factor with three or more levels and a random factor (intLMM) is nested under the corresponding zero-correlation parameter LMM (zcpLMM). It turns out intLMM and zcpLMM are not nested, counter to what was claimed in the RPub.
1.1 Maximum LMMs (maxLMMs)
In the R formula syntax the default intercept 1 (usually omitted in the code) can be replaced with 0 in the fixed-effect and in the random-factor terms of the model. By default 1 causes the estimation of coefficients according to the contrast specification for the factor in the fixed-effect part and the estimation of corresponding correlation parameters in the random-effect terms. When 0 is used, the fixed-effect coefficients are estimates of the means of the factor levels and the correlation parameters estimate the correlations between the levels (i.e., the repeated measures), rather than the correlations relating to contrast-based effects (or means). There are 2 x 2 specifications for the same LMM:
- maxLMM_FE1_RE1:
y ~ 1 + factor + (1 + factor | subj) - maxLMM_FE0_RE1:
y ~ 0 + factor + (1 + factor | subj)
- maxLMM_FE1_RE0:
y ~ 1 + factor + (0 + factor | subj) - maxLMM_FE0_RE0:
y ~ 0 + factor + (0 + factor | subj)
These models work as specified and all of them generate correlation parameters as part of the variance components. The four versions yield identical goodness of fit statistics; they represent reparameterizations of a single LMM.
Alternatively, factor can be replaced with the corresponding columns of the model matrix, that is with numeric covariates. For example, for a factor with three levels A, B, and C with c1 and c2 representing the two contrasts for the three levels, the four completely equivalent model versions are:
- maxLMM_FE1_RE1:
y ~ 1 + c1 + c2 + (1 + c1 + c2 | subj) - maxLMM_FE0_RE1:
y ~ 0 + A + B + C + (1 + c1 + c2 | subj)
- maxLMM_FE1_RE0:
y ~ 1 + c1 + c2 + (0 + A + B + C | subj) - maxLMM_FE0_RE0:
y ~ 0 + A + B + C + (0 + A + B + C | subj)
From the second specification it is transparent, that fixed effects in version 1 and 3 estimate the usual intercept and effects associated with the two contrasts whereas in version 2 and 4 they estimate the means for the three factor levels (no intercept). Similarly, correlation parameters in version 1 and 2 are based on intercept and contrasts whereas in version 3 and version 4 they are based on the measures associated with the levels of the factor. A specification of 1 + A + B + C leads to a redundant model parameter which may cause rank deficiency or convergence problems. 0 + c1 + c2 leads to underspecification of the model.
1.2 Zero-correlation parameters LMMs (zcpLMMs)
Aside from an arguably higher transparency of model parameters, the specification based on numeric covariates is actually required if one wants to use the double-bar syntax to force correlation parameters to zero. Specifically, replacing the | with || does not generate the expected results with factor in the random effect terms (RE). (Alternative: afex::lmer_alt() expands factors in RE to numeric covariates.) A zcpLMM must be based on numeric covariates if one wants to carry out LRT-based comparisons with the corresponding maxLMM or the minimal LMM varying only in intercept (minLMMs).
- zcpLMM_FE1_RE1:
y ~ 1 + c1 + c2 + (1 + c1 + c2 || subj) - zcpLMM_FE0_RE1:
y ~ 0 + A + B + C + (1 + c1 + c2 || subj)
- zcpLMM_FE1_RE0:
y ~ 1 + c1 + c2 + (0 + A + B + C || subj) - zcpLMM_FE0_RE0:
y ~ 0 + A + B + C + (0 + A + B + C || subj)
According to the analyses reported below, LMMs with the same RE structure and different FE structure are reparameterizations, but LMMs with RE1 (5, 6) yield better goodness of fit than LMMs with RE0 (7, 8). Note that all four LMMs use the same number of model parameters (degrees of freedom).
1.3 LMMs with interaction of fixed and random factor (intLMMs)
An LMM which allows for an interaction between the fixed and the random factor is a special case of a zcpLMM. Whether or not the intLMM can be directly compared with the corresponding zcpLMM will be examined in section 8.
- intLMM_FE1:
y ~ 1 + factor + (1 | subj) + (1 | factor:subj) - intLMM_FE0:
y ~ 0 + factor + (1 | subj) + (1 | factor:subj)
The models represent a reparameterization with identical goodness of fit.
1.4 Minimal LMMs varying only in intercept (minLMMs)
For sake of completenss, we also define the two minimal LMMs varying only in intercepts. They also represent a reparameterization of the same LMM.
- minLMM_FE1:
y ~ 1 + factor + (1 | subj) - minLMM_FE0:
y ~ 0 + factor + (1 | subj)
1.5 Nested sequences for model comparison
The analyses reported in the following sections suggest that there are three unique nested sequences for model comparison.
maxLMM_FE1_RE1 -> zcpLMM_FE1_RE1 -> minLMM_FE1_RE1
maxLMM_FE1_RE1 -> intLMM_FE1 -> minLMM_FE1_RE1
maxLMM_FE1_RE0 -> zcpLMM_FE1_RE0 The models not mentioned here are reparamterizations of one of the models listed in the three nested sequences. For example, LMMs with FE0 are always mathematically equivalent to LMMs with FE1 (i.e., the LRT statistics will be identical). The selection of a specific LMM version for a sequence was chosen to allow for comparison of parameter estimates across models within a sequence, but it is important to be aware of possible mathematical equivalences with other versions.
2 Setup
knitr::opts_chunk$set(echo = TRUE)
suppressMessages(library("lme4"))
library("RePsychLing")
suppressMessages(library("tidyverse"))3 Means and correlations of Machines data
The script uses the Machines data which are part of the MEMSS package [efficiency scores of 6 workers (random factor) tested 3 times on each of 3 machines (fixed factor) ]. For later reference, I compute means and correlations between the 3 machine scores by averaging the 3 measures available for the 6 workers (fixed-effect correlations).
In LMMs, means for and/or differences between the 3 machines are estimated with fixed effects. Of course, exactly which means or which differences are estimated depends on the contrast specification for the factor and on the forumula specification. In this script, I use the treatment contrast (i.e., the first machine is used as a reference level for comparison with the second and third machine).
In balanced designs as the current one, the fixed-effect correlations are also listed as part of the summary of lmer()-based maxLMMs. Due to shrinkage LMM correlation parameters are usually somewhat larger than the fixed-effect correlations.
data("Machines", package = "MEMSS")
Machines$Measure <- factor(rep(1:3))
d <- as.tibble(Machines)
contrasts(d$Machine) <- contr.treatment(3)
means <- d %>%
group_by(Machine, Worker) %>% summarize(Score = mean(score))
means## # A tibble: 18 x 3
## # Groups: Machine [?]
## Machine Worker Score
## <fct> <fct> <dbl>
## 1 A 1 52.6
## 2 A 2 52.6
## 3 A 3 59.5
## 4 A 4 51.2
## 5 A 5 51.4
## 6 A 6 46.8
## 7 B 1 62.9
## 8 B 2 59.6
## 9 B 3 68.0
## 10 B 4 62.7
## 11 B 5 65.1
## 12 B 6 43.6
## 13 C 1 67.2
## 14 C 2 61.8
## 15 C 3 70.8
## 16 C 4 64.8
## 17 C 5 71.7
## 18 C 6 61.3# Means of levels and contrasts
l_c <- means %>%
spread(Machine, Score) %>%
mutate(dBA = B-A, dCA = C-A, GM = (A+B+C)/3) %>%
summarize(A = mean(A), B = mean(B), C=mean(C), dBA = mean(dBA), dCA = mean(dCA), GM=mean(GM))
round(l_c, 2)## # A tibble: 1 x 6
## A B C dBA dCA GM
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 52.4 60.3 66.3 7.97 13.9 59.6# Interaction
ggplot(means, aes(x=Worker, y=Score, group=Machine, color=Machine)) +
geom_point() + geom_line() + theme_bw() + ggtitle("Figure 1: Interaction of worker by machine")
# fixed-effext correlations
fe_corr <- means %>%
spread(Machine, Score) %>%
mutate(dBA = B-A, dCA = C-A, GM = (A+B+C)/3) %>%
select(A, B, C, dBA, dCA, GM) %>%
cor()
round(fe_corr, 2)## A B C dBA dCA GM
## A 1.00 0.79 0.61 0.46 -0.37 0.87
## B 0.79 1.00 0.76 0.91 0.03 0.97
## C 0.61 0.76 1.00 0.69 0.51 0.86
## dBA 0.46 0.91 0.69 1.00 0.30 0.81
## dCA -0.37 0.03 0.51 0.30 1.00 0.06
## GM 0.87 0.97 0.86 0.81 0.06 1.004 Maximal LMMs
4.1 maxLMM_FE1_RE1
By default, the LMM takes the contrast assigned with the factor to build the model matrix. For the maxLMM one can use the factor or the numeric covariates (extracted from the default design matrix of the model).
print(summary(m1a <- lmer(score ~ 1 + Machine + (1 + Machine | Worker), d, REML=FALSE)))## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: score ~ 1 + Machine + (1 + Machine | Worker)
## Data: d
##
## AIC BIC logLik deviance df.resid
## 236.4 256.3 -108.2 216.4 44
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.40773 -0.51890 0.03229 0.45599 2.54088
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Worker (Intercept) 13.8157 3.7170
## Machine2 28.6862 5.3559 0.49
## Machine3 11.2431 3.3531 -0.36 0.30
## Residual 0.9246 0.9616
## Number of obs: 54, groups: Worker, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 52.356 1.534 34.124
## Machine2 7.967 2.210 3.605
## Machine3 13.917 1.406 9.899
##
## Correlation of Fixed Effects:
## (Intr) Machn2
## Machine2 0.463
## Machine3 -0.374 0.301# or equivalently with numeric covariates
mm1 <- model.matrix(~ 1 + Machine, d)
dBA <- mm1[, 2]
dCA <- mm1[, 3]
print(summary(m1b <- lmer(score ~ 1 + dBA + dCA + (1 + dBA + dCA | Worker), d, REML=FALSE)))## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: score ~ 1 + dBA + dCA + (1 + dBA + dCA | Worker)
## Data: d
##
## AIC BIC logLik deviance df.resid
## 236.4 256.3 -108.2 216.4 44
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.40773 -0.51890 0.03229 0.45599 2.54088
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Worker (Intercept) 13.8157 3.7170
## dBA 28.6862 5.3559 0.49
## dCA 11.2431 3.3531 -0.36 0.30
## Residual 0.9246 0.9616
## Number of obs: 54, groups: Worker, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 52.356 1.534 34.124
## dBA 7.967 2.210 3.605
## dCA 13.917 1.406 9.899
##
## Correlation of Fixed Effects:
## (Intr) dBA
## dBA 0.463
## dCA -0.374 0.301anova(m1a, m1b)## Data: d
## Models:
## m1a: score ~ 1 + Machine + (1 + Machine | Worker)
## m1b: score ~ 1 + dBA + dCA + (1 + dBA + dCA | Worker)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m1a 10 236.42 256.31 -108.21 216.42
## m1b 10 236.42 256.31 -108.21 216.42 0 0 1# Pick the LMM based on numeric covariates
maxLMM_FE1_RE1 <- m1 <- m1bThus, the intercept estimates the mean for machine A; the two fixed effects estimate the differenes of machines B and C to A (dBA and dCA). The variance components and correlation parameters correspond to this parameterization. (Of course, with a different contrast, fixed effects and correlation parameters will change accordingly.)
They are only similar not identical to the within-subject correlations due to shrinkage. For balanced designs like the current one, the within-subject correlations are actually reported in the fixed effect correlations of the lmer summary (see corresponding correlations in the last chunk.)
4.2 maxLMM_FE0_RE1
We can reparameterize the fixed-effects part.
print(summary(m2a <- lmer(score ~ 0 + Machine + (1 + Machine | Worker), d, REML=FALSE)))## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: score ~ 0 + Machine + (1 + Machine | Worker)
## Data: d
##
## AIC BIC logLik deviance df.resid
## 236.4 256.3 -108.2 216.4 44
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.40773 -0.51890 0.03229 0.45599 2.54088
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Worker (Intercept) 13.8157 3.7170
## Machine2 28.6862 5.3559 0.49
## Machine3 11.2431 3.3531 -0.36 0.30
## Residual 0.9246 0.9616
## Number of obs: 54, groups: Worker, 6
##
## Fixed effects:
## Estimate Std. Error t value
## MachineA 52.356 1.534 34.12
## MachineB 60.322 3.221 18.73
## MachineC 66.272 1.649 40.19
##
## Correlation of Fixed Effects:
## MachnA MachnB
## MachineB 0.794
## MachineC 0.612 0.763# or equivalently with numeric covariates
mm0 <- model.matrix(~ 0 + Machine, d)
A <- mm0[, 1]
B <- mm0[, 2]
C <- mm0[, 3]
print(summary(m2b <- lmer(score ~ 0 + A + B + C + (1 + dBA + dCA | Worker), d, REML=FALSE)))## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: score ~ 0 + A + B + C + (1 + dBA + dCA | Worker)
## Data: d
##
## AIC BIC logLik deviance df.resid
## 236.4 256.3 -108.2 216.4 44
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.40773 -0.51890 0.03229 0.45599 2.54088
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Worker (Intercept) 13.8157 3.7170
## dBA 28.6862 5.3559 0.49
## dCA 11.2431 3.3531 -0.36 0.30
## Residual 0.9246 0.9616
## Number of obs: 54, groups: Worker, 6
##
## Fixed effects:
## Estimate Std. Error t value
## A 52.356 1.534 34.12
## B 60.322 3.221 18.73
## C 66.272 1.649 40.19
##
## Correlation of Fixed Effects:
## A B
## B 0.794
## C 0.612 0.763anova(m2a, m2b)## Data: d
## Models:
## m2a: score ~ 0 + Machine + (1 + Machine | Worker)
## m2b: score ~ 0 + A + B + C + (1 + dBA + dCA | Worker)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2a 10 236.42 256.31 -108.21 216.42
## m2b 10 236.42 256.31 -108.21 216.42 0 0 1# Pick the LMM based on numeric covariates
maxLMM_FE0_RE1 <- m2 <- m2b4.3 maxLMM_FE1_RE0
Correlation parameters are easier to interpret if they estimate the correlations between the levels, because this is very close to a traditional correlation matrix. Of course, this is of no use if one is interested in the correlation of contrast-based experimental effects.
print(summary(m3a <- lmer(score ~ 1 + Machine + (0 + Machine | Worker), d, REML=FALSE)))## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: score ~ 1 + Machine + (0 + Machine | Worker)
## Data: d
##
## AIC BIC logLik deviance df.resid
## 236.4 256.3 -108.2 216.4 44
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.40773 -0.51890 0.03229 0.45599 2.54088
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Worker MachineA 13.8157 3.7170
## MachineB 61.9450 7.8705 0.80
## MachineC 16.0049 4.0006 0.62 0.77
## Residual 0.9246 0.9616
## Number of obs: 54, groups: Worker, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 52.356 1.534 34.124
## Machine2 7.967 2.210 3.605
## Machine3 13.917 1.406 9.899
##
## Correlation of Fixed Effects:
## (Intr) Machn2
## Machine2 0.463
## Machine3 -0.374 0.301# or equivalently with numeric covariates
print(summary(m3b <- lmer(score ~ 1 + dBA + dCA + (0 + A + B + C | Worker), d, REML=FALSE)))## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: score ~ 1 + dBA + dCA + (0 + A + B + C | Worker)
## Data: d
##
## AIC BIC logLik deviance df.resid
## 236.4 256.3 -108.2 216.4 44
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.40773 -0.51890 0.03229 0.45599 2.54088
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Worker A 13.8157 3.7170
## B 61.9450 7.8705 0.80
## C 16.0049 4.0006 0.62 0.77
## Residual 0.9246 0.9616
## Number of obs: 54, groups: Worker, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 52.356 1.534 34.124
## dBA 7.967 2.210 3.605
## dCA 13.917 1.406 9.899
##
## Correlation of Fixed Effects:
## (Intr) dBA
## dBA 0.463
## dCA -0.374 0.301anova(m3a, m3b)## Data: d
## Models:
## m3a: score ~ 1 + Machine + (0 + Machine | Worker)
## m3b: score ~ 1 + dBA + dCA + (0 + A + B + C | Worker)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m3a 10 236.42 256.31 -108.21 216.42
## m3b 10 236.42 256.31 -108.21 216.42 0 0 1# Pick the LMM based on numeric covariates
maxLMM_FE1_RE0 <- m3 <- m3b4.4 maxLMM_FE0_RE0
With this specification, the fixed effects estimate the means of the factor levels and the correlation parameters estimate the corresponding correlations between the scores of the three machines (not correlations based the contrast-based effects). Correlation parameters and fixed-effect correlations are not identical due to shrinkage in the former.
print(summary(m4a <- lmer(score ~ 0 + Machine + (0 + Machine | Worker), d, REML=FALSE)))## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: score ~ 0 + Machine + (0 + Machine | Worker)
## Data: d
##
## AIC BIC logLik deviance df.resid
## 236.4 256.3 -108.2 216.4 44
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.40773 -0.51890 0.03229 0.45599 2.54088
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Worker MachineA 13.8157 3.7170
## MachineB 61.9450 7.8705 0.80
## MachineC 16.0049 4.0006 0.62 0.77
## Residual 0.9246 0.9616
## Number of obs: 54, groups: Worker, 6
##
## Fixed effects:
## Estimate Std. Error t value
## MachineA 52.356 1.534 34.12
## MachineB 60.322 3.221 18.73
## MachineC 66.272 1.649 40.19
##
## Correlation of Fixed Effects:
## MachnA MachnB
## MachineB 0.794
## MachineC 0.612 0.763# or equivalently with numeric covariates
print(summary(m4b <- lmer(score ~ 0 + A + B + C + (0 + A + B + C | Worker), d, REML=FALSE)))## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: score ~ 0 + A + B + C + (0 + A + B + C | Worker)
## Data: d
##
## AIC BIC logLik deviance df.resid
## 236.4 256.3 -108.2 216.4 44
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.40773 -0.51890 0.03229 0.45599 2.54088
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Worker A 13.8157 3.7170
## B 61.9450 7.8705 0.80
## C 16.0049 4.0006 0.62 0.77
## Residual 0.9246 0.9616
## Number of obs: 54, groups: Worker, 6
##
## Fixed effects:
## Estimate Std. Error t value
## A 52.356 1.534 34.12
## B 60.322 3.221 18.73
## C 66.272 1.649 40.19
##
## Correlation of Fixed Effects:
## A B
## B 0.794
## C 0.612 0.763anova(m4a, m4b)## Data: d
## Models:
## m4a: score ~ 0 + Machine + (0 + Machine | Worker)
## m4b: score ~ 0 + A + B + C + (0 + A + B + C | Worker)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m4a 10 236.42 256.31 -108.21 216.42
## m4b 10 236.42 256.31 -108.21 216.42 0 0 1# Pick the LMM based on numeric covariates
maxLMM_FE0_RE0 <- m4 <- m4b4.5 Equivalence of maxLMMs
All maxLMMs yield identical goodness of fit; they are reparameterizations of one LMM.
anova(m1a, m2a, m3a, m4a)## Data: d
## Models:
## m1a: score ~ 1 + Machine + (1 + Machine | Worker)
## m2a: score ~ 0 + Machine + (1 + Machine | Worker)
## m3a: score ~ 1 + Machine + (0 + Machine | Worker)
## m4a: score ~ 0 + Machine + (0 + Machine | Worker)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m1a 10 236.42 256.31 -108.21 216.42
## m2a 10 236.42 256.31 -108.21 216.42 0 0 <2e-16 ***
## m3a 10 236.42 256.31 -108.21 216.42 0 0 <2e-16 ***
## m4a 10 236.42 256.31 -108.21 216.42 0 0 1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(m1b, m2b, m3b, m4b)## Data: d
## Models:
## m1b: score ~ 1 + dBA + dCA + (1 + dBA + dCA | Worker)
## m2b: score ~ 0 + A + B + C + (1 + dBA + dCA | Worker)
## m3b: score ~ 1 + dBA + dCA + (0 + A + B + C | Worker)
## m4b: score ~ 0 + A + B + C + (0 + A + B + C | Worker)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m1b 10 236.42 256.31 -108.21 216.42
## m2b 10 236.42 256.31 -108.21 216.42 0 0 <2e-16 ***
## m3b 10 236.42 256.31 -108.21 216.42 0 0 <2e-16 ***
## m4b 10 236.42 256.31 -108.21 216.42 0 0 1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 15 Zero-correlation parameter LMMs
5.1 zcpLMM_FE1_RE1
One can force correlation parameters to zero with the double-bar syntax. Again, I illustrate the two versions. The double-bar syntax of LMM m5 is a short cut for the explicit specification of LMM m5b. As shown in m5c, the double-bar syntax does not work for the factor-based specification (see RPub for an elaborate illustration).
One recent option is to use afex::lmer_alt() (see m5d) Here the factors are expanded to numeric covariates and the result is identical to m5. (It also reports p-values for fixed effects KR or S correction.)
print(summary(m5 <- lmer(score ~ 1 + dBA + dCA + (1 + dBA + dCA || Worker), d, REML=FALSE))) ## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: score ~ 1 + dBA + dCA + ((1 | Worker) + (0 + dBA | Worker) +
## (0 + dCA | Worker))
## Data: d
##
## AIC BIC logLik deviance df.resid
## 235.1 249.0 -110.6 221.1 47
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.31194 -0.50699 0.00438 0.45299 2.47124
##
## Random effects:
## Groups Name Variance Std.Dev.
## Worker (Intercept) 13.7771 3.7118
## Worker.1 dBA 28.8275 5.3691
## Worker.2 dCA 10.9321 3.3064
## Residual 0.9265 0.9626
## Number of obs: 54, groups: Worker, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 52.356 1.532 34.170
## dBA 7.967 2.215 3.596
## dCA 13.917 1.387 10.031
##
## Correlation of Fixed Effects:
## (Intr) dBA
## dBA -0.015
## dCA -0.024 0.017print(summary(m5b <- lmer(score ~ 1 + dBA + dCA + (1 | Worker) + (0 + dBA | Worker) + (0 + dCA | Worker),
data=d, REML=FALSE)))## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: score ~ 1 + dBA + dCA + (1 | Worker) + (0 + dBA | Worker) + (0 +
## dCA | Worker)
## Data: d
##
## AIC BIC logLik deviance df.resid
## 235.1 249.0 -110.6 221.1 47
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.31194 -0.50699 0.00438 0.45299 2.47124
##
## Random effects:
## Groups Name Variance Std.Dev.
## Worker (Intercept) 13.7771 3.7118
## Worker.1 dBA 28.8275 5.3691
## Worker.2 dCA 10.9321 3.3064
## Residual 0.9265 0.9626
## Number of obs: 54, groups: Worker, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 52.356 1.532 34.170
## dBA 7.967 2.215 3.596
## dCA 13.917 1.387 10.031
##
## Correlation of Fixed Effects:
## (Intr) dBA
## dBA -0.015
## dCA -0.024 0.017print(summary(m5c <- lmer(score ~ 1 + Machine + (1 + Machine || Worker), d, REML=FALSE))) ## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: score ~ 1 + Machine + ((1 | Worker) + (0 + Machine | Worker))
## Data: d
##
## AIC BIC logLik deviance df.resid
## 238.4 260.3 -108.2 216.4 43
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.40773 -0.51890 0.03229 0.45599 2.54088
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Worker (Intercept) 0.7011 0.8373
## Worker.1 MachineA 13.1147 3.6214
## MachineB 61.2439 7.8259 0.81
## MachineC 15.3038 3.9120 0.61 0.77
## Residual 0.9246 0.9616
## Number of obs: 54, groups: Worker, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 52.356 1.534 34.124
## Machine2 7.967 2.210 3.605
## Machine3 13.917 1.406 9.899
##
## Correlation of Fixed Effects:
## (Intr) Machn2
## Machine2 0.463
## Machine3 -0.374 0.301print(summary(m5d <- afex::lmer_alt(score ~ 1 + Machine + (1 + Machine || Worker), d, REML=FALSE))) ## Linear mixed model fit by maximum likelihood . t-tests use
## Satterthwaite's method [lmerModLmerTest]
## Formula: score ~ 1 + Machine + (1 + re1.Machine2 + re1.Machine3 || Worker)
## Data: data
##
## AIC BIC logLik deviance df.resid
## 235.1 249.0 -110.6 221.1 47
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.31194 -0.50699 0.00438 0.45299 2.47124
##
## Random effects:
## Groups Name Variance Std.Dev.
## Worker (Intercept) 13.7771 3.7118
## Worker.1 re1.Machine2 28.8275 5.3691
## Worker.2 re1.Machine3 10.9321 3.3064
## Residual 0.9265 0.9626
## Number of obs: 54, groups: Worker, 6
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 52.356 1.532 6.007 34.170 4.12e-08 ***
## Machine2 7.967 2.215 6.004 3.596 0.0114 *
## Machine3 13.917 1.387 6.007 10.031 5.65e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Machn2
## Machine2 -0.015
## Machine3 -0.024 0.017anova(m5, m5b, m5c)## Data: d
## Models:
## m5: score ~ 1 + dBA + dCA + ((1 | Worker) + (0 + dBA | Worker) +
## m5: (0 + dCA | Worker))
## m5b: score ~ 1 + dBA + dCA + (1 | Worker) + (0 + dBA | Worker) + (0 +
## m5b: dCA | Worker)
## m5c: score ~ 1 + Machine + ((1 | Worker) + (0 + Machine | Worker))
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m5 7 235.11 249.04 -110.56 221.11
## m5b 7 235.11 249.04 -110.56 221.11 0.0000 0 1.00
## m5c 11 238.42 260.30 -108.21 216.42 4.6954 4 0.32# Pick the LMM based on numeric covariates
zcpLMM_FE1_RE1 <- m55.2 zcpLMM_FE0_RE1
print(summary(m6 <- lmer(score ~ 0 + A + B + C + (1 + dBA + dCA || Worker), d, REML=FALSE)))## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: score ~ 0 + A + B + C + ((1 | Worker) + (0 + dBA | Worker) +
## (0 + dCA | Worker))
## Data: d
##
## AIC BIC logLik deviance df.resid
## 235.1 249.0 -110.6 221.1 47
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.31194 -0.50699 0.00438 0.45299 2.47124
##
## Random effects:
## Groups Name Variance Std.Dev.
## Worker (Intercept) 13.7771 3.7118
## Worker.1 dBA 28.8275 5.3691
## Worker.2 dCA 10.9321 3.3064
## Residual 0.9265 0.9626
## Number of obs: 54, groups: Worker, 6
##
## Fixed effects:
## Estimate Std. Error t value
## A 52.356 1.532 34.17
## B 60.322 2.674 22.56
## C 66.272 2.042 32.45
##
## Correlation of Fixed Effects:
## A B
## B 0.560
## C 0.734 0.420zcpLMM_FE0_RE1 <- m65.3 zcpLMM_FE1_RE0
print(summary(m7 <- lmer(score ~ 1 + dBA + dCA + (0 + A + B + C || Worker), d, REML=FALSE)))## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: score ~ 1 + dBA + dCA + ((0 + A | Worker) + (0 + B | Worker) +
## (0 + C | Worker))
## Data: d
##
## AIC BIC logLik deviance df.resid
## 241.6 255.5 -113.8 227.6 47
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.28287 -0.54219 0.01523 0.47815 2.36982
##
## Random effects:
## Groups Name Variance Std.Dev.
## Worker A 13.8157 3.7170
## Worker.1 B 61.9450 7.8705
## Worker.2 C 16.0049 4.0006
## Residual 0.9246 0.9616
## Number of obs: 54, groups: Worker, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 52.356 1.534 34.124
## dBA 7.967 3.568 2.233
## dCA 13.917 2.252 6.179
##
## Correlation of Fixed Effects:
## (Intr) dBA
## dBA -0.430
## dCA -0.681 0.293zcpLMM_FE1_RE0 <- m75.4 zcpLMM_FE0_RE0
print(summary(m8 <- lmer(score ~ 0 + A + B + C + (0 + A + B + C || Worker), d, REML=FALSE)))## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: score ~ 0 + A + B + C + ((0 + A | Worker) + (0 + B | Worker) +
## (0 + C | Worker))
## Data: d
##
## AIC BIC logLik deviance df.resid
## 241.6 255.5 -113.8 227.6 47
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.28287 -0.54219 0.01523 0.47815 2.36982
##
## Random effects:
## Groups Name Variance Std.Dev.
## Worker A 13.8157 3.7170
## Worker.1 B 61.9450 7.8705
## Worker.2 C 16.0049 4.0006
## Residual 0.9246 0.9616
## Number of obs: 54, groups: Worker, 6
##
## Fixed effects:
## Estimate Std. Error t value
## A 52.356 1.534 34.12
## B 60.322 3.221 18.73
## C 66.272 1.649 40.19
##
## Correlation of Fixed Effects:
## A B
## B 0.000
## C 0.000 0.000zcpLMM_FE0_RE0 <- m85.5 (Non-)Equivalence of zcpLMMs
anova(m5, m6, m7, m8)## Data: d
## Models:
## m5: score ~ 1 + dBA + dCA + ((1 | Worker) + (0 + dBA | Worker) +
## m5: (0 + dCA | Worker))
## m6: score ~ 0 + A + B + C + ((1 | Worker) + (0 + dBA | Worker) +
## m6: (0 + dCA | Worker))
## m7: score ~ 1 + dBA + dCA + ((0 + A | Worker) + (0 + B | Worker) +
## m7: (0 + C | Worker))
## m8: score ~ 0 + A + B + C + ((0 + A | Worker) + (0 + B | Worker) +
## m8: (0 + C | Worker))
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m5 7 235.11 249.04 -110.56 221.11
## m6 7 235.11 249.04 -110.56 221.11 0 0 <2e-16 ***
## m7 7 241.63 255.55 -113.81 227.63 0 0 1
## m8 7 241.63 255.55 -113.81 227.63 0 0 1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The different random-effect terms yield different deviances with a better fit for random-effect terms based on contrasts (RE1) than on levels (RE0). Given the equivalence of all maxLMMs, this translates into larger differences in goodness of fit between zcpLMM and maxLMM for random-effect terms based on levels.
anova(m5, m1)## Data: d
## Models:
## m5: score ~ 1 + dBA + dCA + ((1 | Worker) + (0 + dBA | Worker) +
## m5: (0 + dCA | Worker))
## m1: score ~ 1 + dBA + dCA + (1 + dBA + dCA | Worker)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m5 7 235.11 249.04 -110.56 221.11
## m1 10 236.42 256.31 -108.21 216.42 4.6954 3 0.1955anova(m6, m2)## Data: d
## Models:
## m6: score ~ 0 + A + B + C + ((1 | Worker) + (0 + dBA | Worker) +
## m6: (0 + dCA | Worker))
## m2: score ~ 0 + A + B + C + (1 + dBA + dCA | Worker)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m6 7 235.11 249.04 -110.56 221.11
## m2 10 236.42 256.31 -108.21 216.42 4.6954 3 0.1955anova(m7, m3)## Data: d
## Models:
## m7: score ~ 1 + dBA + dCA + ((0 + A | Worker) + (0 + B | Worker) +
## m7: (0 + C | Worker))
## m3: score ~ 1 + dBA + dCA + (0 + A + B + C | Worker)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m7 7 241.63 255.55 -113.81 227.63
## m3 10 236.42 256.31 -108.21 216.42 11.208 3 0.01065 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(m8, m4)## Data: d
## Models:
## m8: score ~ 0 + A + B + C + ((0 + A | Worker) + (0 + B | Worker) +
## m8: (0 + C | Worker))
## m4: score ~ 0 + A + B + C + (0 + A + B + C | Worker)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m8 7 241.63 255.55 -113.81 227.63
## m4 10 236.42 256.31 -108.21 216.42 11.208 3 0.01065 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 16 Interaction between fixed and random factor (intLMMs)
The intLMM is of interest in the current context, because it estimates two VPs, one for the intercept of the random factor and one for the interaction of the fixed and random factor. There is no correlation parameter. Therefore, it looks like a special case of a zcpLMM, but Section 8 shows that intLMM is not nested under zcpLMM.
The interaction between a fixed and a random factor is estimated with the following LMMs for FE1 and FE0; there are no RE0 versions for intLMM.
# intLMM_FE1
print(summary(m9 <- lmer(score ~ 1 + Machine + (1 | Worker) + (1 | Machine:Worker), d, REML=FALSE)))## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: score ~ 1 + Machine + (1 | Worker) + (1 | Machine:Worker)
## Data: d
##
## AIC BIC logLik deviance df.resid
## 237.3 249.2 -112.6 225.3 48
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.25899 -0.55135 -0.00649 0.44670 2.54889
##
## Random effects:
## Groups Name Variance Std.Dev.
## Machine:Worker (Intercept) 11.5398 3.3970
## Worker (Intercept) 19.0487 4.3645
## Residual 0.9246 0.9616
## Number of obs: 54, groups: Machine:Worker, 18; Worker, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 52.356 2.269 23.072
## Machine2 7.967 1.987 4.009
## Machine3 13.917 1.987 7.003
##
## Correlation of Fixed Effects:
## (Intr) Machn2
## Machine2 -0.438
## Machine3 -0.438 0.500# ... or equivalently
print(summary(m9b <- lmer(score ~ 1 + Machine + (1 | Worker/Machine), d, REML=FALSE)))## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: score ~ 1 + Machine + (1 | Worker/Machine)
## Data: d
##
## AIC BIC logLik deviance df.resid
## 237.3 249.2 -112.6 225.3 48
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.25899 -0.55135 -0.00649 0.44670 2.54889
##
## Random effects:
## Groups Name Variance Std.Dev.
## Machine:Worker (Intercept) 11.5398 3.3970
## Worker (Intercept) 19.0487 4.3645
## Residual 0.9246 0.9616
## Number of obs: 54, groups: Machine:Worker, 18; Worker, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 52.356 2.269 23.072
## Machine2 7.967 1.987 4.009
## Machine3 13.917 1.987 7.003
##
## Correlation of Fixed Effects:
## (Intr) Machn2
## Machine2 -0.438
## Machine3 -0.438 0.500anova(m9, m9b)## Data: d
## Models:
## m9: score ~ 1 + Machine + (1 | Worker) + (1 | Machine:Worker)
## m9b: score ~ 1 + Machine + (1 | Worker/Machine)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m9 6 237.27 249.2 -112.64 225.27
## m9b 6 237.27 249.2 -112.64 225.27 0 0 1intLMM_FE1 <- m9
# intLMM_FE0
print(summary(m10 <- lmer(score ~ 0 + Machine + (1 | Worker) + (1 | Machine:Worker), d, REML=FALSE)))## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: score ~ 0 + Machine + (1 | Worker) + (1 | Machine:Worker)
## Data: d
##
## AIC BIC logLik deviance df.resid
## 237.3 249.2 -112.6 225.3 48
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.25899 -0.55135 -0.00649 0.44670 2.54889
##
## Random effects:
## Groups Name Variance Std.Dev.
## Machine:Worker (Intercept) 11.5398 3.3970
## Worker (Intercept) 19.0487 4.3645
## Residual 0.9246 0.9616
## Number of obs: 54, groups: Machine:Worker, 18; Worker, 6
##
## Fixed effects:
## Estimate Std. Error t value
## MachineA 52.356 2.269 23.07
## MachineB 60.322 2.269 26.58
## MachineC 66.272 2.269 29.20
##
## Correlation of Fixed Effects:
## MachnA MachnB
## MachineB 0.617
## MachineC 0.617 0.617intLMM_FE0 <- m10
anova(m9, m10)## Data: d
## Models:
## m9: score ~ 1 + Machine + (1 | Worker) + (1 | Machine:Worker)
## m10: score ~ 0 + Machine + (1 | Worker) + (1 | Machine:Worker)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m9 6 237.27 249.2 -112.64 225.27
## m10 6 237.27 249.2 -112.64 225.27 0 0 17 Minimal LMMs varying only in intercept (minLMMs)
Finally, for sake of completeness, here are the two minimal LMMs varying only in intercept.
# minLMM_FE1
print(summary(m11 <- lmer(score ~ 1 + Machine + (1 | Worker), d, REML=FALSE)))## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: score ~ 1 + Machine + (1 | Worker)
## Data: d
##
## AIC BIC logLik deviance df.resid
## 303.7 313.6 -146.9 293.7 49
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.8012 -0.5245 0.1352 0.6740 1.7760
##
## Random effects:
## Groups Name Variance Std.Dev.
## Worker (Intercept) 21.93 4.683
## Residual 9.58 3.095
## Number of obs: 54, groups: Worker, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 52.356 2.046 25.584
## Machine2 7.967 1.032 7.722
## Machine3 13.917 1.032 13.489
##
## Correlation of Fixed Effects:
## (Intr) Machn2
## Machine2 -0.252
## Machine3 -0.252 0.500minLMM_FE1 <- m11
# minLMM_FE0
print(summary(m12 <- lmer(score ~ 0 + Machine + (1 | Worker), d, REML=FALSE)))## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: score ~ 0 + Machine + (1 | Worker)
## Data: d
##
## AIC BIC logLik deviance df.resid
## 303.7 313.6 -146.9 293.7 49
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.8012 -0.5245 0.1352 0.6740 1.7760
##
## Random effects:
## Groups Name Variance Std.Dev.
## Worker (Intercept) 21.93 4.683
## Residual 9.58 3.095
## Number of obs: 54, groups: Worker, 6
##
## Fixed effects:
## Estimate Std. Error t value
## MachineA 52.356 2.046 25.58
## MachineB 60.322 2.046 29.48
## MachineC 66.272 2.046 32.38
##
## Correlation of Fixed Effects:
## MachnA MachnB
## MachineB 0.873
## MachineC 0.873 0.873minLMM_FE0 <- m12
anova(m11, m12)## Data: d
## Models:
## m11: score ~ 1 + Machine + (1 | Worker)
## m12: score ~ 0 + Machine + (1 | Worker)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m11 5 303.7 313.65 -146.85 293.7
## m12 5 303.7 313.65 -146.85 293.7 0 0 18 Is intLMM nested under zcpLMM or not?
Here I will recapitulate my intial analysis the RPub which I interpreted as support for a nested structure. Given the counterarguments by Rune Haubo in this R-Sig-ME thread, I look for counterexamples with the present data.
8.1 Initial analysis with zcpLMM_FE1_RE1
In the RPub (see LMM m1d vs. m1e), I suggested that intLMM is nested under zcpLMMs for a factors with three or more levels, bascially assuming that the VP for the interaction is divided into separate VPs (one for each contrast) in the zcpLMM. The LRT for the Machines data shows that including the variance component for the interaction significantly improves the fit. Figure 1 shows that the profiles (across workers) for the three machines are indeed far from parallel.
Given that there are two contrasts for the Machine factor, I thought that one can think of the Machine:Worker term as a combination/aggregation of two VPs. I also thought that zcpLMM m5 (which uses contrasts for fixed effects and random-effect terms) represents an LMM in which the two components are estimated separately. From this perspective, LMM m9 appears to be nested under zcpLMM m5.
- maxLMM_FE1_RE1 -> zcpLMM_FE1_RE1 -> intLMM_FE1 -> minLMM_FE1,
where “->” means reduction of LMM complexity.
VarCorr(intLMM_FE1) # m9## Groups Name Std.Dev.
## Machine:Worker (Intercept) 3.39703
## Worker (Intercept) 4.36448
## Residual 0.96158VarCorr(zcpLMM_FE1_RE1) # m5## Groups Name Std.Dev.
## Worker (Intercept) 3.71176
## Worker.1 dBA 5.36913
## Worker.2 dCA 3.30638
## Residual 0.96257anova(m11, m9, m5, m1)## Data: d
## Models:
## m11: score ~ 1 + Machine + (1 | Worker)
## m9: score ~ 1 + Machine + (1 | Worker) + (1 | Machine:Worker)
## m5: score ~ 1 + dBA + dCA + ((1 | Worker) + (0 + dBA | Worker) +
## m5: (0 + dCA | Worker))
## m1: score ~ 1 + dBA + dCA + (1 + dBA + dCA | Worker)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m11 5 303.70 313.65 -146.85 293.70
## m9 6 237.27 249.20 -112.64 225.27 68.4338 1 < 2e-16 ***
## m5 7 235.11 249.04 -110.56 221.11 4.1562 1 0.04148 *
## m1 10 236.42 256.31 -108.21 216.42 4.6954 3 0.19551
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1My interpretation of this result was that the interaction between worker and the first contrast of Machine (i.e., dBA) is significantly larger than the interaction between worker and the second contrast of Machine (i.e., dCA), that is variance component dBA = 5.37 > dCA = 3.31 in zcpLMM m5. Rune Haubo’s post at R-Sig-ME suggested that this may not be true. Therefore, I carried out a few additional checks.
8.2 Test with zcpLMM_FE1_RE0
One option is to test the sequence for the alternative zcpLMM_FE1_RE0 specification. The minLMM_FE1 is not included, because there is no intercept in this
VarCorr(zcpLMM_FE1_RE0) # m7## Groups Name Std.Dev.
## Worker A 3.71695
## Worker.1 B 7.87051
## Worker.2 C 4.00061
## Residual 0.96158anova(m9, m7, m3)## Data: d
## Models:
## m9: score ~ 1 + Machine + (1 | Worker) + (1 | Machine:Worker)
## m7: score ~ 1 + dBA + dCA + ((0 + A | Worker) + (0 + B | Worker) +
## m7: (0 + C | Worker))
## m3: score ~ 1 + dBA + dCA + (0 + A + B + C | Worker)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m9 6 237.27 249.20 -112.64 225.27
## m7 7 241.63 255.55 -113.81 227.63 0.000 1 1.00000
## m3 10 236.42 256.31 -108.21 216.42 11.208 3 0.01065 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1In this sequence, intLMM m9 fits better than zcpLMM m7, despite the extra degree of freedom for the latter. This is a problem, because in this case intLMM m9 cannot be nested under zcpLMM m5. This is strong evidence against a nested structure.
8.3 Test with two-level factor
For a two-level factor, intLMM and zcpLMM estimate the same number of model parameters. Therefore, the two models should represent a reparameterization, that is they should be mathematically equivalent. As a quick check, I used the data from only two machines.
d2 <- droplevels(subset(d, Machine != "C"))
contrasts(d2$Machine) <- contr.treatment(2)
mm_treatment <- model.matrix(~ 1 + Machine, d2)
dBA_d2 <- mm_treatment[,2]
print(summary(int1_d2 <- lmer(score ~ 1 + dBA_d2 + (1 | Worker) + (1 | Machine:Worker), d2, REML=FALSE)))## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: score ~ 1 + dBA_d2 + (1 | Worker) + (1 | Machine:Worker)
## Data: d2
##
## AIC BIC logLik deviance df.resid
## 169.8 177.7 -79.9 159.8 31
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.02351 -0.50430 0.06411 0.42752 2.29980
##
## Random effects:
## Groups Name Variance Std.Dev.
## Machine:Worker (Intercept) 14.26 3.777
## Worker (Intercept) 23.54 4.852
## Residual 1.16 1.077
## Number of obs: 36, groups: Machine:Worker, 12; Worker, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 52.356 2.523 20.753
## dBA_d2 7.967 2.210 3.605
##
## Correlation of Fixed Effects:
## (Intr)
## dBA_d2 -0.438print(summary(zcp1_d2 <- lmer(score ~ 1 + dBA_d2 + (1 + dBA_d2 || Worker), d2, REML=FALSE)))## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: score ~ 1 + dBA_d2 + ((1 | Worker) + (0 + dBA_d2 | Worker))
## Data: d2
##
## AIC BIC logLik deviance df.resid
## 165.2 173.1 -77.6 155.2 31
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.06919 -0.52887 0.05533 0.42967 2.22646
##
## Random effects:
## Groups Name Variance Std.Dev.
## Worker (Intercept) 13.992 3.741
## Worker.1 dBA_d2 29.057 5.390
## Residual 1.155 1.075
## Number of obs: 36, groups: Worker, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 52.356 1.548 33.823
## dBA_d2 7.967 2.230 3.573
##
## Correlation of Fixed Effects:
## (Intr)
## dBA_d2 -0.019anova(int1_d2, zcp1_d2)## Data: d2
## Models:
## int1_d2: score ~ 1 + dBA_d2 + (1 | Worker) + (1 | Machine:Worker)
## zcp1_d2: score ~ 1 + dBA_d2 + ((1 | Worker) + (0 + dBA_d2 | Worker))
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## int1_d2 5 169.76 177.68 -79.879 159.76
## zcp1_d2 5 165.17 173.09 -77.587 155.17 4.5845 0 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1This was not expected: intLMM and zcpLMM are not equivalent! Does the contrast specification matter? I replaced the treatment contrast with a sum contrast.
contrasts(d2$Machine) <- contr.sum(2)
mm_sum <- model.matrix(~ 1 + Machine, d2)
dA_GM <- mm_sum[,2]
print(summary(int_d2 <- lmer(score ~ 1 + dA_GM + (1 | Worker) + (1 | Machine:Worker), d2, REML=FALSE)))## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: score ~ 1 + dA_GM + (1 | Worker) + (1 | Machine:Worker)
## Data: d2
##
## AIC BIC logLik deviance df.resid
## 169.8 177.7 -79.9 159.8 31
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.02351 -0.50430 0.06411 0.42752 2.29980
##
## Random effects:
## Groups Name Variance Std.Dev.
## Machine:Worker (Intercept) 14.26 3.777
## Worker (Intercept) 23.54 4.852
## Residual 1.16 1.077
## Number of obs: 36, groups: Machine:Worker, 12; Worker, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 56.339 2.268 24.841
## dA_GM -3.983 1.105 -3.605
##
## Correlation of Fixed Effects:
## (Intr)
## dA_GM 0.000print(summary(zcp_d2 <- lmer(score ~ 1 + dA_GM + (1 + dA_GM || Worker), d2, REML=FALSE)))## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: score ~ 1 + dA_GM + ((1 | Worker) + (0 + dA_GM | Worker))
## Data: d2
##
## AIC BIC logLik deviance df.resid
## 169.8 177.7 -79.9 159.8 31
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.02351 -0.50430 0.06411 0.42752 2.29980
##
## Random effects:
## Groups Name Variance Std.Dev.
## Worker (Intercept) 30.670 5.538
## Worker.1 dA_GM 7.132 2.671
## Residual 1.160 1.077
## Number of obs: 36, groups: Worker, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 56.339 2.268 24.841
## dA_GM -3.983 1.105 -3.605
##
## Correlation of Fixed Effects:
## (Intr)
## dA_GM 0.000anova(int_d2, zcp_d2)## Data: d2
## Models:
## int_d2: score ~ 1 + dA_GM + (1 | Worker) + (1 | Machine:Worker)
## zcp_d2: score ~ 1 + dA_GM + ((1 | Worker) + (0 + dA_GM | Worker))
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## int_d2 5 169.76 177.68 -79.879 159.76
## zcp_d2 5 169.76 177.68 -79.879 159.76 0 0 1Indeed, contrast matters! With contr.sum(2), intLMM and zcpLMM are clearly identical. Thus, it appears that for two-level factors mathematical equivalence is possible for sum contrasts, but not for treatment contrasts. My hunch is that the intercept in the treatment contrast estimates one of the factor levels.
8.4 Test with modified Machines data
Figure 1 suggests that Workers 5 and 6 are mainly responsible for this difference. In response to Rune Haubo’s post and as a check, I changed their Machine B data to align parallel to Machine C.
For a similar exercise see also “Interactions of grouping factors and other covariates” in Bates (2009, slides 84 to 95). A synopsis of this section is given in Section 11 below.
d$score_2 <- ifelse(d$Worker == "6" & d$Machine == "B", d$score + 15, d$score)
d$score_2 <- ifelse(d$Worker == "5" & d$Machine == "B", d$score + 4, d$score_2)
means_2 <- d %>%
group_by(Machine, Worker) %>% summarize(Score = mean(score_2))
means_2## # A tibble: 18 x 3
## # Groups: Machine [?]
## Machine Worker Score
## <fct> <fct> <dbl>
## 1 A 1 52.6
## 2 A 2 52.6
## 3 A 3 59.5
## 4 A 4 51.2
## 5 A 5 51.4
## 6 A 6 46.8
## 7 B 1 62.9
## 8 B 2 59.6
## 9 B 3 68.0
## 10 B 4 62.7
## 11 B 5 69.1
## 12 B 6 58.6
## 13 C 1 67.2
## 14 C 2 61.8
## 15 C 3 70.8
## 16 C 4 64.8
## 17 C 5 71.7
## 18 C 6 61.3ggplot(means_2, aes(x=Worker, y=Score, group=Machine, color=Machine)) +
geom_point() + geom_line() + theme_bw() + ggtitle("Figure 2: Modified means of worker x machine")
contrasts(d$Machine) <- contr.treatment(3)
mm1 <- model.matrix(~ 1 + Machine, d)
dBA <- mm1[, 2]
dCA <- mm1[, 3]
print(summary(m9_2 <- lmer(score_2 ~ Machine + (1 | Worker) + (1 | Machine:Worker), d, REML=FALSE)))## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: score_2 ~ Machine + (1 | Worker) + (1 | Machine:Worker)
## Data: d
##
## AIC BIC logLik deviance df.resid
## 220.5 232.5 -104.3 208.5 48
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.34708 -0.52782 0.05395 0.50259 2.41783
##
## Random effects:
## Groups Name Variance Std.Dev.
## Machine:Worker (Intercept) 3.6593 1.9129
## Worker (Intercept) 11.2694 3.3570
## Residual 0.9246 0.9616
## Number of obs: 54, groups: Machine:Worker, 18; Worker, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 52.356 1.594 32.854
## Machine2 11.133 1.150 9.681
## Machine3 13.917 1.150 12.101
##
## Correlation of Fixed Effects:
## (Intr) Machn2
## Machine2 -0.361
## Machine3 -0.361 0.500print(summary(m5_2 <- lmer(score_2 ~ 1 + dBA + dCA + (1 + dBA + dCA || Worker), d, REML=FALSE))) ## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: score_2 ~ 1 + dBA + dCA + ((1 | Worker) + (0 + dBA | Worker) +
## (0 + dCA | Worker))
## Data: d
##
## AIC BIC logLik deviance df.resid
## 229.0 242.9 -107.5 215.0 47
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.37755 -0.53888 0.03744 0.49481 2.31515
##
## Random effects:
## Groups Name Variance Std.Dev.
## Worker (Intercept) 13.3031 3.6473
## Worker.1 dBA 9.9696 3.1575
## Worker.2 dCA 10.3992 3.2248
## Residual 0.9394 0.9692
## Number of obs: 54, groups: Worker, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 52.356 1.506 34.754
## dBA 11.133 1.329 8.378
## dCA 13.917 1.356 10.266
##
## Correlation of Fixed Effects:
## (Intr) dBA
## dBA -0.026
## dCA -0.026 0.029print(summary(m10_2 <- lmer(score_2 ~ Machine + (1 | Worker), d, REML=FALSE)))## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: score_2 ~ Machine + (1 | Worker)
## Data: d
##
## AIC BIC logLik deviance df.resid
## 254 264 -122 244 49
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.4522 -0.5530 0.0685 0.5359 1.7610
##
## Random effects:
## Groups Name Variance Std.Dev.
## Worker (Intercept) 12.184 3.491
## Residual 3.669 1.915
## Number of obs: 54, groups: Worker, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 52.3556 1.4948 35.02
## Machine2 11.1333 0.6385 17.44
## Machine3 13.9167 0.6385 21.80
##
## Correlation of Fixed Effects:
## (Intr) Machn2
## Machine2 -0.214
## Machine3 -0.214 0.500print(summary(m1_2 <- lmer(score_2 ~ 1 + Machine + (1 + Machine | Worker), d, REML=FALSE)))## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: score_2 ~ 1 + Machine + (1 + Machine | Worker)
## Data: d
##
## AIC BIC logLik deviance df.resid
## 215.8 235.7 -97.9 195.8 44
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.2654 -0.6364 0.1528 0.5789 2.4308
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Worker (Intercept) 13.8234 3.7180
## Machine2 10.8421 3.2927 -0.39
## Machine3 11.2807 3.3587 -0.36 1.00
## Residual 0.9016 0.9495
## Number of obs: 54, groups: Worker, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 52.356 1.534 34.124
## Machine2 11.133 1.381 8.062
## Machine3 13.917 1.407 9.889
##
## Correlation of Fixed Effects:
## (Intr) Machn2
## Machine2 -0.404
## Machine3 -0.373 0.974anova(m10_2, m9_2, m5_2, m1_2)## Data: d
## Models:
## m10_2: score_2 ~ Machine + (1 | Worker)
## m9_2: score_2 ~ Machine + (1 | Worker) + (1 | Machine:Worker)
## m5_2: score_2 ~ 1 + dBA + dCA + ((1 | Worker) + (0 + dBA | Worker) +
## m5_2: (0 + dCA | Worker))
## m1_2: score_2 ~ 1 + Machine + (1 + Machine | Worker)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m10_2 5 254.02 263.97 -122.012 244.02
## m9_2 6 220.53 232.46 -104.263 208.53 35.497 1 2.554e-09 ***
## m5_2 7 228.99 242.91 -107.494 214.99 0.000 1 1.0000000
## m1_2 10 215.79 235.68 -97.894 195.79 19.200 3 0.0002485 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(m10_2, m9_2, m1_2)## Data: d
## Models:
## m10_2: score_2 ~ Machine + (1 | Worker)
## m9_2: score_2 ~ Machine + (1 | Worker) + (1 | Machine:Worker)
## m1_2: score_2 ~ 1 + Machine + (1 + Machine | Worker)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m10_2 5 254.02 263.97 -122.012 244.02
## m9_2 6 220.53 232.46 -104.263 208.53 35.497 1 2.554e-09 ***
## m1_2 10 215.79 235.68 -97.894 195.79 12.740 4 0.01262 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(m10_2, m5_2, m1_2)## Data: d
## Models:
## m10_2: score_2 ~ Machine + (1 | Worker)
## m5_2: score_2 ~ 1 + dBA + dCA + ((1 | Worker) + (0 + dBA | Worker) +
## m5_2: (0 + dCA | Worker))
## m1_2: score_2 ~ 1 + Machine + (1 + Machine | Worker)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m10_2 5 254.02 263.97 -122.012 244.02
## m5_2 7 228.99 242.91 -107.494 214.99 29.037 2 4.951e-07 ***
## m1_2 10 215.79 235.68 -97.894 195.79 19.200 3 0.0002485 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The modification worked in principle, that is variance components are now quite similar, dBA = 3.16 anddCA = 3.22 in zcpLMM m5_2. However, intLMM m9_2 fits better than zcpLMM m5_2, despite the extra degree of freedom for the latter. This is a problem, because in this case intLMM m9_2 cannot be nested under zcpLMM m5_2. Again, this is evidence agaist intLMM being nested under zcpLMM.
Is the contrast the problem? In other words, do the fit statistics align with the degrees of freedom for the sum contrast?
contrasts(d$Machine) <- contr.sum(3)
mm3 <- model.matrix(~ 1 + Machine, d)
dA_GM <- mm3[, 2]
dB_GM <- mm3[, 3]
print(summary(m9_2b <- lmer(score_2 ~ Machine + (1 | Worker) + (1 | Machine:Worker), d, REML=FALSE)))## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: score_2 ~ Machine + (1 | Worker) + (1 | Machine:Worker)
## Data: d
##
## AIC BIC logLik deviance df.resid
## 220.5 232.5 -104.3 208.5 48
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.34708 -0.52782 0.05395 0.50259 2.41783
##
## Random effects:
## Groups Name Variance Std.Dev.
## Machine:Worker (Intercept) 3.6593 1.9129
## Worker (Intercept) 11.2694 3.3570
## Residual 0.9246 0.9616
## Number of obs: 54, groups: Machine:Worker, 18; Worker, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 60.706 1.449 41.904
## Machine1 -8.350 0.664 -12.576
## Machine2 2.783 0.664 4.192
##
## Correlation of Fixed Effects:
## (Intr) Machn1
## Machine1 0.000
## Machine2 0.000 -0.500print(summary(m5_2b <- lmer(score_2 ~ 1 + dA_GM + dB_GM + (+ dA_GM + dB_GM || Worker), d, REML=FALSE))) ## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula:
## score_2 ~ 1 + dA_GM + dB_GM + ((1 | Worker) + (0 + dA_GM | Worker) +
## (0 + dB_GM | Worker))
## Data: d
##
## AIC BIC logLik deviance df.resid
## 223.6 237.6 -104.8 209.6 47
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.44418 -0.56151 0.08045 0.48211 2.39276
##
## Random effects:
## Groups Name Variance Std.Dev.
## Worker (Intercept) 12.488 3.5338
## Worker.1 dA_GM 4.504 2.1223
## Worker.2 dB_GM 1.023 1.0116
## Residual 0.938 0.9685
## Number of obs: 54, groups: Worker, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 60.7056 1.4487 41.904
## dA_GM -8.3500 0.8863 -9.422
## dB_GM 2.7833 0.4531 6.143
##
## Correlation of Fixed Effects:
## (Intr) dA_GM
## dA_GM 0.000
## dB_GM 0.000 -0.043print(summary(m10_2b <- lmer(score_2 ~ Machine + (1 | Worker), d, REML=FALSE)))## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: score_2 ~ Machine + (1 | Worker)
## Data: d
##
## AIC BIC logLik deviance df.resid
## 254 264 -122 244 49
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.4522 -0.5530 0.0685 0.5359 1.7610
##
## Random effects:
## Groups Name Variance Std.Dev.
## Worker (Intercept) 12.184 3.491
## Residual 3.669 1.915
## Number of obs: 54, groups: Worker, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 60.7056 1.4487 41.90
## Machine1 -8.3500 0.3686 -22.65
## Machine2 2.7833 0.3686 7.55
##
## Correlation of Fixed Effects:
## (Intr) Machn1
## Machine1 0.000
## Machine2 0.000 -0.500print(summary(m1_2b <- lmer(score_2 ~ 1 + Machine + (1 + Machine | Worker), d, REML=FALSE)))## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: score_2 ~ 1 + Machine + (1 + Machine | Worker)
## Data: d
##
## AIC BIC logLik deviance df.resid
## 215.8 235.7 -97.9 195.8 44
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.2654 -0.6364 0.1528 0.5789 2.4308
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Worker (Intercept) 12.4917 3.5344
## Machine1 4.9143 2.2168 -0.23
## Machine2 1.1598 1.0769 0.18 -1.00
## Residual 0.9016 0.9495
## Number of obs: 54, groups: Worker, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 60.7056 1.4487 41.904
## Machine1 -8.3500 0.9233 -9.044
## Machine2 2.7833 0.4761 5.846
##
## Correlation of Fixed Effects:
## (Intr) Machn1
## Machine1 -0.223
## Machine2 0.163 -0.942anova(m10_2b, m9_2b, m5_2b, m1_2b)## Data: d
## Models:
## m10_2b: score_2 ~ Machine + (1 | Worker)
## m9_2b: score_2 ~ Machine + (1 | Worker) + (1 | Machine:Worker)
## m5_2b: score_2 ~ 1 + dA_GM + dB_GM + ((1 | Worker) + (0 + dA_GM | Worker) +
## m5_2b: (0 + dB_GM | Worker))
## m1_2b: score_2 ~ 1 + Machine + (1 + Machine | Worker)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m10_2b 5 254.02 263.97 -122.012 244.02
## m9_2b 6 220.53 232.46 -104.263 208.53 35.497 1 2.554e-09 ***
## m5_2b 7 223.64 237.56 -104.819 209.64 0.000 1 1.000000
## m1_2b 10 215.79 235.68 -97.894 195.79 13.850 3 0.003116 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(m10_2b, m9_2b, m1_2b)## Data: d
## Models:
## m10_2b: score_2 ~ Machine + (1 | Worker)
## m9_2b: score_2 ~ Machine + (1 | Worker) + (1 | Machine:Worker)
## m1_2b: score_2 ~ 1 + Machine + (1 + Machine | Worker)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m10_2b 5 254.02 263.97 -122.012 244.02
## m9_2b 6 220.53 232.46 -104.263 208.53 35.497 1 2.554e-09 ***
## m1_2b 10 215.79 235.68 -97.894 195.79 12.740 4 0.01262 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(m10_2b, m5_2b, m1_2b)## Data: d
## Models:
## m10_2b: score_2 ~ Machine + (1 | Worker)
## m5_2b: score_2 ~ 1 + dA_GM + dB_GM + ((1 | Worker) + (0 + dA_GM | Worker) +
## m5_2b: (0 + dB_GM | Worker))
## m1_2b: score_2 ~ 1 + Machine + (1 + Machine | Worker)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m10_2b 5 254.02 263.97 -122.012 244.02
## m5_2b 7 223.64 237.56 -104.819 209.64 34.387 2 3.412e-08 ***
## m1_2b 10 215.79 235.68 -97.894 195.79 13.850 3 0.003116 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The results are clear. IntLMM and zcpLMM were not in the expected order with respect to goodness of fit; intLMM fits the modified data better than the zcpLMM. For a two-level factor mathematical equivalence was possible for sum contrast, not treatment contrast. For the three-level factor the outcome did not depend on the choice of contrast for factor. Also when using sum contrast for the original data, a better goodness-of-fit statistic was obtained for intLMM than zcpLMM.
8.5 Conclusion
For the LMMs considered in this script, there are three hiearchical structures that can be compared with LRTs:
maxLMM_FE1_RE1 -> zcpLMM_FE1_RE1 -> minLMM_FE1_RE1
maxLMM_FE1_RE1 -> intLMM_FE1 -> minLMM_FE1_RE1
maxLMM_FE1_RE0 -> zcpLMM_FE1_RE0 where “->” means reduction of LMM complexity.
Empirically (from what I have seen so far), goodness of fit is better for zcpLMM_RE1 than zcpLMM_RE0. Note that these LMMs use the same number of model parameters (degrees of freedom).
In conclusion, intLMM and zcpLMM are not in a hierarchal order despite the extra degree of freedom for the latter. Rune Haubo’s argument that there are two hierarchical orders with respect to these models desribed in the R-Sig-ME thread is valid.
9 Parsimonious LMMs: Treatment contrasts
Parsimonious LMMs still need to be sorted into the hierarchical model sequences. Here I explore whether treatment contrasts in prsmLMM_FE1_RE1 lead to mathematically equivalent goodness of fit in prsmLMM_FE1_RE0.
9.1 dBA and dCA vs. B and C
print(summary(prsmLMM_FE1_RE1 <- lmer(score ~ 1 + Machine + (1 | Worker) + (0 + dBA + dCA | Worker), d, REML=FALSE)))## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: score ~ 1 + Machine + (1 | Worker) + (0 + dBA + dCA | Worker)
## Data: d
##
## AIC BIC logLik deviance df.resid
## 236.6 252.5 -110.3 220.6 46
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.29836 -0.50665 -0.01621 0.45950 2.42303
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Worker (Intercept) 13.7151 3.7034
## Worker.1 dBA 29.1133 5.3957
## dCA 11.0432 3.3231 0.30
## Residual 0.9257 0.9621
## Number of obs: 54, groups: Worker, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 59.6500 1.8021 33.101
## Machine1 -7.2944 0.9893 -7.373
## Machine2 0.6722 1.4119 0.476
##
## Correlation of Fixed Effects:
## (Intr) Machn1
## Machine1 -0.530
## Machine2 0.383 -0.710print(summary(prsmLMM_FE1_RE0 <- lmer(score ~ 1 + Machine + (1 | Worker) + (0 + B + C | Worker), d, REML=FALSE)))## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: score ~ 1 + Machine + (1 | Worker) + (0 + B + C | Worker)
## Data: d
##
## AIC BIC logLik deviance df.resid
## 236.6 252.5 -110.3 220.6 46
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.29836 -0.50665 -0.01621 0.45950 2.42303
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Worker (Intercept) 13.7151 3.7034
## Worker.1 B 29.1133 5.3957
## C 11.0432 3.3231 0.30
## Residual 0.9257 0.9621
## Number of obs: 54, groups: Worker, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 59.6500 1.8021 33.101
## Machine1 -7.2944 0.9893 -7.373
## Machine2 0.6722 1.4119 0.476
##
## Correlation of Fixed Effects:
## (Intr) Machn1
## Machine1 -0.530
## Machine2 0.383 -0.710anova(prsmLMM_FE1_RE1, prsmLMM_FE1_RE0, maxLMM_FE1_RE1, maxLMM_FE1_RE0)## Data: d
## Models:
## prsmLMM_FE1_RE1: score ~ 1 + Machine + (1 | Worker) + (0 + dBA + dCA | Worker)
## prsmLMM_FE1_RE0: score ~ 1 + Machine + (1 | Worker) + (0 + B + C | Worker)
## maxLMM_FE1_RE1: score ~ 1 + dBA + dCA + (1 + dBA + dCA | Worker)
## maxLMM_FE1_RE0: score ~ 1 + dBA + dCA + (0 + A + B + C | Worker)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## prsmLMM_FE1_RE1 8 236.58 252.49 -110.29 220.58
## prsmLMM_FE1_RE0 8 236.58 252.49 -110.29 220.58 0.0000 0 1.0000
## maxLMM_FE1_RE1 10 236.42 256.31 -108.21 216.42 4.1578 2 0.1251
## maxLMM_FE1_RE0 10 236.42 256.31 -108.21 216.42 0.0000 0 <2e-16
##
## prsmLMM_FE1_RE1
## prsmLMM_FE1_RE0
## maxLMM_FE1_RE1
## maxLMM_FE1_RE0 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1prsmLMM_FE1_RE1 and prsmLMM_FE1_RE0 are mathematically equivalent.
9.2 dBA and dCA vs. A and B
print(summary(prsmLMM_FE1_RE1 <- lmer(score ~ 1 + Machine + (1 | Worker) + (0 + dBA + dCA | Worker), d, REML=FALSE)))## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: score ~ 1 + Machine + (1 | Worker) + (0 + dBA + dCA | Worker)
## Data: d
##
## AIC BIC logLik deviance df.resid
## 236.6 252.5 -110.3 220.6 46
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.29836 -0.50665 -0.01621 0.45950 2.42303
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Worker (Intercept) 13.7151 3.7034
## Worker.1 dBA 29.1133 5.3957
## dCA 11.0432 3.3231 0.30
## Residual 0.9257 0.9621
## Number of obs: 54, groups: Worker, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 59.6500 1.8021 33.101
## Machine1 -7.2944 0.9893 -7.373
## Machine2 0.6722 1.4119 0.476
##
## Correlation of Fixed Effects:
## (Intr) Machn1
## Machine1 -0.530
## Machine2 0.383 -0.710print(summary(prsmLMM_FE1_RE0 <- lmer(score ~ 1 + Machine + (1 | Worker) + (0 + A + B | Worker), d, REML=FALSE)))## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: score ~ 1 + Machine + (1 | Worker) + (0 + A + B | Worker)
## Data: d
##
## AIC BIC logLik deviance df.resid
## 237.4 253.3 -110.7 221.4 46
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.29838 -0.55106 0.01563 0.46128 2.54804
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Worker (Intercept) 15.7729 3.9715
## Worker.1 A 10.9837 3.3142
## B 29.6279 5.4432 0.33
## Residual 0.9269 0.9627
## Number of obs: 54, groups: Worker, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 59.6500 1.9023 31.356
## Machine1 -7.2944 0.9772 -7.465
## Machine2 0.6722 1.4109 0.476
##
## Correlation of Fixed Effects:
## (Intr) Machn1
## Machine1 -0.017
## Machine2 0.374 -0.703anova(prsmLMM_FE1_RE1, prsmLMM_FE1_RE0, maxLMM_FE1_RE1, maxLMM_FE1_RE0)## Data: d
## Models:
## prsmLMM_FE1_RE1: score ~ 1 + Machine + (1 | Worker) + (0 + dBA + dCA | Worker)
## prsmLMM_FE1_RE0: score ~ 1 + Machine + (1 | Worker) + (0 + A + B | Worker)
## maxLMM_FE1_RE1: score ~ 1 + dBA + dCA + (1 + dBA + dCA | Worker)
## maxLMM_FE1_RE0: score ~ 1 + dBA + dCA + (0 + A + B + C | Worker)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## prsmLMM_FE1_RE1 8 236.58 252.49 -110.29 220.58
## prsmLMM_FE1_RE0 8 237.40 253.31 -110.70 221.40 0.0000 0 1.0000
## maxLMM_FE1_RE1 10 236.42 256.31 -108.21 216.42 4.9802 2 0.0829
## maxLMM_FE1_RE0 10 236.42 256.31 -108.21 216.42 0.0000 0 <2e-16
##
## prsmLMM_FE1_RE1
## prsmLMM_FE1_RE0
## maxLMM_FE1_RE1 .
## maxLMM_FE1_RE0 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1prsmLMM_FE1_RE1 and prsmLMM_FE1_RE0 are mathematically not equivalent.
9.3 dBA and dCA vs. A and B
print(summary(prsmLMM_FE1_RE1 <- lmer(score ~ 1 + Machine + (1 | Worker) + (0 + dBA + dCA | Worker), d, REML=FALSE)))## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: score ~ 1 + Machine + (1 | Worker) + (0 + dBA + dCA | Worker)
## Data: d
##
## AIC BIC logLik deviance df.resid
## 236.6 252.5 -110.3 220.6 46
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.29836 -0.50665 -0.01621 0.45950 2.42303
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Worker (Intercept) 13.7151 3.7034
## Worker.1 dBA 29.1133 5.3957
## dCA 11.0432 3.3231 0.30
## Residual 0.9257 0.9621
## Number of obs: 54, groups: Worker, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 59.6500 1.8021 33.101
## Machine1 -7.2944 0.9893 -7.373
## Machine2 0.6722 1.4119 0.476
##
## Correlation of Fixed Effects:
## (Intr) Machn1
## Machine1 -0.530
## Machine2 0.383 -0.710print(summary(prsmLMM_FE1_RE0 <- lmer(score ~ 1 + Machine + (1 | Worker) + (0 + A + C | Worker), d, REML=FALSE)))## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: score ~ 1 + Machine + (1 | Worker) + (0 + A + C | Worker)
## Data: d
##
## AIC BIC logLik deviance df.resid
## 245.4 261.3 -114.7 229.4 46
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.2582 -0.4981 -0.0030 0.4498 2.5500
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Worker (Intercept) 61.0492 7.8134
## Worker.1 A 28.2979 5.3196
## C 28.9294 5.3786 0.80
## Residual 0.9269 0.9627
## Number of obs: 54, groups: Worker, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 59.6500 3.4790 17.146
## Machine1 -7.2944 0.9813 -7.433
## Machine2 0.6722 1.3949 0.482
##
## Correlation of Fixed Effects:
## (Intr) Machn1
## Machine1 0.275
## Machine2 -0.394 -0.698anova(prsmLMM_FE1_RE1, prsmLMM_FE1_RE0, maxLMM_FE1_RE1, maxLMM_FE1_RE0)## Data: d
## Models:
## prsmLMM_FE1_RE1: score ~ 1 + Machine + (1 | Worker) + (0 + dBA + dCA | Worker)
## prsmLMM_FE1_RE0: score ~ 1 + Machine + (1 | Worker) + (0 + A + C | Worker)
## maxLMM_FE1_RE1: score ~ 1 + dBA + dCA + (1 + dBA + dCA | Worker)
## maxLMM_FE1_RE0: score ~ 1 + dBA + dCA + (0 + A + B + C | Worker)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## prsmLMM_FE1_RE1 8 236.58 252.49 -110.29 220.58
## prsmLMM_FE1_RE0 8 245.43 261.35 -114.72 229.43 0.000 0 1.000000
## maxLMM_FE1_RE1 10 236.42 256.31 -108.21 216.42 13.016 2 0.001492
## maxLMM_FE1_RE0 10 236.42 256.31 -108.21 216.42 0.000 0 < 2.2e-16
##
## prsmLMM_FE1_RE1
## prsmLMM_FE1_RE0
## maxLMM_FE1_RE1 **
## maxLMM_FE1_RE0 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1prsmLMM_FE1_RE1 and prsmLMM_FE1_RE0 are mathematically not equivalent.
9.4 Summary
Using the reference level (A) in prsmLMM_FE1_RE0 destroys mathematical equivalence. This is a good example that contrast specifications matter for mathematical equivalence.
10 Parsimonious LMMs: Sum contrasts
Here I explore whether sum contrasts in prsmLMM_FE1_RE1 lead to mathematically equivalent goodness of fit in prsmLMM_FE1_RE0.
10.1 dA_GM and dB_GM vs. A and B / A and C / B and C
contrasts(d$Machine) <- contr.sum(3)
print(summary(prsmLMM_FE1_RE1 <- lmer(score ~ 1 + dA_GM + dB_GM + (1 | Worker) + (0 + dA_GM + dB_GM | Worker), d, REML=FALSE)))## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: score ~ 1 + dA_GM + dB_GM + (1 | Worker) + (0 + dA_GM + dB_GM |
## Worker)
## Data: d
##
## AIC BIC logLik deviance df.resid
## 239.5 255.5 -111.8 223.5 46
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.27099 -0.50449 -0.00538 0.45602 2.51023
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Worker (Intercept) 22.8953 4.7849
## Worker.1 dA_GM 5.6165 2.3699
## dB_GM 11.6388 3.4116 -0.71
## Residual 0.9246 0.9616
## Number of obs: 54, groups: Worker, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 59.6500 1.9578 30.468
## dA_GM -7.2944 0.9851 -7.405
## dB_GM 0.6722 1.4050 0.478
##
## Correlation of Fixed Effects:
## (Intr) dA_GM
## dA_GM 0.000
## dB_GM 0.000 -0.701print(summary(prsmLMM_FE1_RE0 <- lmer(score ~ 1 + Machine + (1 | Worker) + (0 + B + C | Worker), d, REML=FALSE)))## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: score ~ 1 + Machine + (1 | Worker) + (0 + B + C | Worker)
## Data: d
##
## AIC BIC logLik deviance df.resid
## 236.6 252.5 -110.3 220.6 46
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.29836 -0.50665 -0.01621 0.45950 2.42303
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Worker (Intercept) 13.7151 3.7034
## Worker.1 B 29.1133 5.3957
## C 11.0432 3.3231 0.30
## Residual 0.9257 0.9621
## Number of obs: 54, groups: Worker, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 59.6500 1.8021 33.101
## Machine1 -7.2944 0.9893 -7.373
## Machine2 0.6722 1.4119 0.476
##
## Correlation of Fixed Effects:
## (Intr) Machn1
## Machine1 -0.530
## Machine2 0.383 -0.710anova(prsmLMM_FE1_RE1, prsmLMM_FE1_RE0, maxLMM_FE1_RE1, maxLMM_FE1_RE0)## Data: d
## Models:
## prsmLMM_FE1_RE1: score ~ 1 + dA_GM + dB_GM + (1 | Worker) + (0 + dA_GM + dB_GM |
## prsmLMM_FE1_RE1: Worker)
## prsmLMM_FE1_RE0: score ~ 1 + Machine + (1 | Worker) + (0 + B + C | Worker)
## maxLMM_FE1_RE1: score ~ 1 + dBA + dCA + (1 + dBA + dCA | Worker)
## maxLMM_FE1_RE0: score ~ 1 + dBA + dCA + (0 + A + B + C | Worker)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## prsmLMM_FE1_RE1 8 239.54 255.46 -111.77 223.54
## prsmLMM_FE1_RE0 8 236.58 252.49 -110.29 220.58 2.9698 0 <2e-16
## maxLMM_FE1_RE1 10 236.42 256.31 -108.21 216.42 4.1578 2 0.1251
## maxLMM_FE1_RE0 10 236.42 256.31 -108.21 216.42 0.0000 0 <2e-16
##
## prsmLMM_FE1_RE1
## prsmLMM_FE1_RE0 ***
## maxLMM_FE1_RE1
## maxLMM_FE1_RE0 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 110.2 Summary
For sum contrasts, prsmLMM_FE1_RE1 and prsmLMM_FE1_RE0 are never mathematically equivalent (3 tests).
11 Bates (2008) on interaction LMM (intLMM)
Bates (2008, Slide 84-95) elaborated on “[i]nteractions of grouping factors and other covariates”. He distinguihes between two ways to define the interaction between a fixed factor and a random (grouping) factor:
- Define simple, scalar random effects for
Workerand for theWorker:Machineinteraction, yielding random effects forWorkerandWorker/Machine.
intLMM <- fm1 <- lmer(score ~ 1 + Machine + (1 | Worker) + (1 | Worker:Machine), Machines)- Define vector-valued random effects for
Worker
max_LMM_RE0 <- fm2 <- lmer(score ~ 1 + Machine + (0 + Machine | Worker), Machines)
Bates (2008, Slide 91) writes: “Although not obvious from the specifications, the model fits are nested. If the variance-covariance matrix for the vector-valued random effects [fm2, max_LMM_RE0] has a special form, called compound symmetry, the model reduces to model [fm1, int_LMM_RE0].” He goes on to show that fm2 fits marginally better than fm1, but it is the data from Worker 6 on Machine B that cause the necessity of the more complex model.
Bates’s (2008, Slide 95) main concern are “[T]radeoffs when defining interactions”. Specifically, taken verbatim from the slide:
- It is important to realize that estimating scale parameters (i.e. variances and covariances) is considerably more difficult than estimating location parameters (i.e. means or fixed-effects coefficients).
- A vector-valued random effect term having q_i random effects per level of the grouping factor requires q_i(q_i + 1)/2 variance-covariance parameters to be estimated. A simple, scalar random effect for the interaction of a “random-effects” factor and a “fixed-effects” factor requires only 1 additional variance-covariance parameter.
- Especially when the “fixed-effects” factor has a moderate to large number of levels, the trade-off in model complexity argues against the vector-valued approach.
- One of the major sources of difficulty in using the
lme4package is the tendency to overspecify the number of random effects per level of a grouping factor.
12 Appendix
sessionInfo()## R version 3.5.0 (2018-04-23)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS High Sierra 10.13.4
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] bindrcpp_0.2.2 forcats_0.3.0 stringr_1.3.1
## [4] dplyr_0.7.5 purrr_0.2.4 readr_1.1.1
## [7] tidyr_0.8.1 tibble_1.4.2 ggplot2_2.2.1
## [10] tidyverse_1.2.1 RePsychLing_0.0.4 lme4_1.1-18
## [13] Matrix_1.2-14
##
## loaded via a namespace (and not attached):
## [1] httr_1.3.1 jsonlite_1.5 splines_3.5.0
## [4] carData_3.0-1 modelr_0.1.2 assertthat_0.2.0
## [7] stats4_3.5.0 coin_1.2-2 cellranger_1.1.0
## [10] yaml_2.1.19 numDeriv_2016.8-1 pillar_1.2.3
## [13] backports_1.1.2 lattice_0.20-35 glue_1.2.0
## [16] digest_0.6.15 rvest_0.3.2 minqa_1.2.4
## [19] sandwich_2.4-0 colorspace_1.3-2 htmltools_0.3.6
## [22] plyr_1.8.4 psych_1.8.4 pkgconfig_2.0.1
## [25] broom_0.4.4 haven_1.1.1 xtable_1.8-2
## [28] mvtnorm_1.0-7 scales_0.5.0 openxlsx_4.0.17
## [31] rio_0.5.10 emmeans_1.2.1 car_3.0-0
## [34] TH.data_1.0-8 lazyeval_0.2.1 cli_1.0.0
## [37] mnormt_1.5-5 survival_2.42-3 magrittr_1.5
## [40] crayon_1.3.4 readxl_1.1.0 estimability_1.3
## [43] evaluate_0.10.1 nlme_3.1-137 MASS_7.3-50
## [46] xml2_1.2.0 foreign_0.8-70 tools_3.5.0
## [49] data.table_1.11.2 hms_0.4.2 multcomp_1.4-8
## [52] munsell_0.4.3 compiler_3.5.0 rlang_0.2.0
## [55] grid_3.5.0 nloptr_1.0.4 rstudioapi_0.7
## [58] labeling_0.3 rmarkdown_1.9 lmerTest_3.0-1.9001
## [61] codetools_0.2-15 gtable_0.2.0 abind_1.4-5
## [64] curl_3.2 reshape2_1.4.3 R6_2.2.2
## [67] zoo_1.8-1 lubridate_1.7.4 knitr_1.20
## [70] afex_0.20-2 utf8_1.1.4 bindr_0.1.1
## [73] rprojroot_1.3-2 modeltools_0.2-21 stringi_1.2.2
## [76] parallel_3.5.0 Rcpp_0.12.17 coda_0.19-1
## [79] tidyselect_0.2.4