In courtiol/LM2GLMM: Data science for biologists: Generalized linear modelling with R
library(LM2GLMM)
library(doSNOW)
spaMM::spaMM.options(nb_cores = parallel::detectCores() - 1)
options(width = 120)
knitr::opts_chunk$set(cache = TRUE, cache.path = "./cache_knitr/MM_intro/", fig.path = "./fig_knitr/MM_intro/", fig.width = 5, fig.height = 5, fig.align = "center")
Mixed-effects models
- 4.0 Introduction to LMM & GLMM
- 4.1 Solving LM problems using LMM
- 4.2 A showcase of some useful applications
[Back to main menu](./Title.html#2)
You will learn in this session r .emo("goal")
- what mixed models (LMM & GLMM) are
- how to use {lme4}, {spaMM} and {glmmTMB} to fit mixed models
- how to make predictions, tests, CI...
- how to check assumptions
- when to consider effects as fixed or random
- how to implement different random effects structure
Linear Mixed-effects Models
Why mixed models? r .emo("info")
Goal:
To study, or to account for, unobservable sources of heterogeneity between observations.
Mixed-effects models allow for:
- the study of other questions than LM/GLM (e.g. heritability)
- the fixing of assumption violations in LM/GLM (lack of independence, some cases of overdispersion)
- the reduction of the uncertainty in estimates and predictions in cases where many parameters would have to be estimated at the cost of an additional hypothesis (the distribution of the random effects)
The main sources of heterogeneity considered by mixed-effects models are:
- origin (in its widest sense)
- time
- space
The linear mixed-effects model r .emo("info")
Definition
- a LMM is a specific linear model for which, given the design matrix $\mathbf{X}$, the responses ($\mathbf{Y}$) are no longer independent, but where the correlations can be described in terms of a random effect, i.e. a random variable that is not included in the predictor variables
Mathematical notation of LM r .emo("info")
LM, which we have seen before:
$\mathbf{Y} = \mathbf{X} \beta + \epsilon$
$$
\begin{bmatrix} y_1 \ y_2 \ y_3 \ \vdots \ y_n \end{bmatrix} =
\begin{bmatrix}
1 & x_{1,1} & x_{1,2} & \dots & x_{1,p} \
1 & x_{2,1} & x_{2,2} & \dots & x_{2,p} \
1 & x_{3,1} & x_{3,2} & \dots & x_{3,p} \
\vdots & \vdots & \vdots & \ddots & \vdots \
1 & x_{n,1} & x_{n,2} & \dots & x_{n,p}
\end{bmatrix}
\begin{bmatrix}
\beta_0 \ \beta_1 \ \beta_2 \ \vdots \ \beta_p
\end{bmatrix}+
\begin{bmatrix}
\epsilon_1 \ \epsilon_2 \ \epsilon_3 \ \vdots \ \epsilon_n
\end{bmatrix}
$$ with,
$$
\epsilon \sim
\mathcal{N}\left(0,
\begin{bmatrix}
\sigma^2 & 0 & 0 & \dots & 0 \
0 & \sigma^2 & 0 & \dots & 0 \
0 & 0 & \sigma^2 & \dots & 0 \
\vdots & \vdots & \vdots & \ddots & \vdots \
0 & 0 & 0 & \dots & \sigma^2
\end{bmatrix}\right)
$$
Mathematical notation of LMM r .emo("info")
LMM with one random factor with $q$ levels:
$\mathbf{Y} = \mathbf{X} \beta + \mathbf{Z}b + \epsilon$
$$
\begin{bmatrix} y_1 \ y_2 \ y_3 \ \vdots \ y_n \end{bmatrix} =
\begin{bmatrix}
1 & x_{1,1} & x_{1,2} & \dots & x_{1,p} \
1 & x_{2,1} & x_{2,2} & \dots & x_{2,p} \
1 & x_{3,1} & x_{3,2} & \dots & x_{3,p} \
\vdots & \vdots & \vdots & \ddots & \vdots \
1 & x_{n,1} & x_{n,2} & \dots & x_{n,p}
\end{bmatrix}
\begin{bmatrix}
\beta_0 \ \beta_1 \ \beta_2 \ \vdots \ \beta_p
\end{bmatrix}+
\begin{bmatrix}
z_{1,1} & z_{1,2} & z_{1,3} & \dots & z_{1,q} \
z_{2,1} & z_{2,2} & z_{2,3} & \dots & z_{2,q} \
z_{3,1} & z_{3,2} & z_{3,3} & \dots & z_{3,q} \
\vdots & \vdots & \vdots & \ddots & \vdots \
z_{n,1} & z_{n,2} & z_{n,3} & \dots & z_{n,q} \
\end{bmatrix}
\begin{bmatrix}
b_1 \ b_2 \ b_3 \ \vdots \ b_q
\end{bmatrix}+
\begin{bmatrix}
\epsilon_1 \ \epsilon_2 \ \epsilon_3 \ \vdots \ \epsilon_n
\end{bmatrix}
$$ with,
$$
b \sim
\mathcal{N}\left(0,
\begin{bmatrix}
c_{1,1} & c_{1,2} & c_{1,3} & \dots & c_{1,q} \
c_{2,1} & c_{2,2} & c_{2,3} & \dots & c_{2,q} \
c_{3,1} & c_{3,2} & c_{3,3} & \dots & c_{3,q} \
\vdots & \vdots & \vdots & \ddots & \vdots \
c_{q,1} & c_{q,2} & c_{q,3} & \dots & c_{q,q} \
\end{bmatrix}\right)
\text{& }
\epsilon \sim
\mathcal{N}\left(0,
\begin{bmatrix}
\phi & 0 & 0 & \dots & 0 \
0 & \phi & 0 & \dots & 0 \
0 & 0 & \phi & \dots & 0 \
\vdots & \vdots & \vdots & \ddots & \vdots \
0 & 0 & 0 & \dots & \phi
\end{bmatrix}\right)
$$
with $\text{E}(b) = 0$, $\text{Cov}(b) = \mathbf{C}$, which is symmetrical ($c_{i, j} = c_{j, i}$). Also, $\text{Cov}(b, \epsilon) = 0$.
Mathematical notation of LMM r .emo("info")
LMM with one random factor with $q$ levels:
$\mathbf{Y} = \mathbf{X} \beta + \mathbf{Z}b + \epsilon$
$$
\begin{bmatrix} y_1 \ y_2 \ y_3 \ \vdots \ y_n \end{bmatrix} =
\begin{bmatrix}
1 & x_{1,1} & x_{1,2} & \dots & x_{1,p} \
1 & x_{2,1} & x_{2,2} & \dots & x_{2,p} \
1 & x_{3,1} & x_{3,2} & \dots & x_{3,p} \
\vdots & \vdots & \vdots & \ddots & \vdots \
1 & x_{n,1} & x_{n,2} & \dots & x_{n,p}
\end{bmatrix}
\begin{bmatrix}
\beta_0 \ \beta_1 \ \beta_2 \ \vdots \ \beta_p
\end{bmatrix}+
\begin{bmatrix}
z_{1,1} & z_{1,2} & z_{1,3} & \dots & z_{1,q} \
z_{2,1} & z_{2,2} & z_{2,3} & \dots & z_{2,q} \
z_{3,1} & z_{3,2} & z_{3,3} & \dots & z_{3,q} \
\vdots & \vdots & \vdots & \ddots & \vdots \
z_{n,1} & z_{n,2} & z_{n,3} & \dots & z_{n,q} \
\end{bmatrix}
\begin{bmatrix}
b_1 \ b_2 \ b_3 \ \vdots \ b_q
\end{bmatrix}+
\begin{bmatrix}
\epsilon_1 \ \epsilon_2 \ \epsilon_3 \ \vdots \ \epsilon_n
\end{bmatrix}
$$ and often,
$$
b \sim
\mathcal{N}\left(0,
\begin{bmatrix}
\lambda & 0 & 0 & \dots & 0 \
0 & \lambda & 0 & \dots & 0 \
0 & 0 & \lambda & \dots & 0 \
\vdots & \vdots & \vdots & \ddots & \vdots \
0 & 0 & 0 & \dots & \lambda \
\end{bmatrix}\right)
\text{& }
\epsilon \sim
\mathcal{N}\left(0,
\begin{bmatrix}
\phi & 0 & 0 & \dots & 0 \
0 & \phi & 0 & \dots & 0 \
0 & 0 & \phi & \dots & 0 \
\vdots & \vdots & \vdots & \ddots & \vdots \
0 & 0 & 0 & \dots & \phi
\end{bmatrix}\right)
$$
Note: the dimensions of the matrices differ: $n \times n$ for $\mathbf{X}$ vs $q \times q$ for $\mathbf{Z}$
Fitting procedure
A simple simulation function r .emo("alien")
simulate_Mix <- function(intercept, slope, n, group_nb, var.rand, var.error){
data <- data.frame(intercept = intercept, slope = slope, x = runif(n))
group_compo <- rmultinom(n = 1, size = n, prob = c(rep(1/group_nb, group_nb)))
data$group <- factor(rep(paste("group", 1:group_nb, sep = "_"), group_compo))
data$b <- rep(rnorm(group_nb, mean = 0, sd = sqrt(var.rand)), group_compo)
data$error <- rnorm(n, mean = 0, sd = sqrt(var.error))
data$y <- data$intercept + data$slope*data$x + data$b + data$error
return(data)
}
set.seed(1)
Aliens <- simulate_Mix(intercept = 50, slope = 1.5, n = 30, group_nb = 10, var.rand = 2, var.error = 0.5)
Note: look carefully at the next slide to better understand.
Our toy dataset r .emo("alien")
Aliens
Fitting the model with {lme4} r .emo("practice")
library(lme4)
(fit_lme4 <- lmer(y ~ x + (1|group), data = Aliens)) ## REML fit by default!
print(VarCorr(fit_lme4), comp = "Variance") ## extract REML variance estimates
Fitting the model with {spaMM} r .emo("practice")
library(spaMM)
(fit_spaMM <- fitme(y ~ x + (1|group), data = Aliens, method = "REML"))
Functions for fitting the mixed model numerically r .emo("nerd")
lik_b <- function(b.vec, level, intercept, slope, var.rand, var.error, data, scale = 1) {
lik <- sapply(b.vec, function(b){
sub_data <- data[which(data$group == level), ]
sub_data$pred <- intercept + slope*sub_data$x + b
sub_data$conditional.density <- dnorm(sub_data$y, mean = sub_data$pred, sd = sqrt(var.error))
return(dnorm(b, mean = 0, sd = sqrt(var.rand)) * prod(sub_data$conditional.density))
})
return(scale * lik)
}
log_lik_b_prod <- function(param, data, scale = 1){
log_lik_vec <- sapply(levels(data$group), function(level) {
log(integrate(lik_b, -Inf, Inf, level = level, intercept = param[1], slope = param[2],
var.rand = param[3], var.error = param[4], data = data)$value)})
return(scale * sum(log_lik_vec))
}
Note: this is a ML fit for simplicity.
Testing the functions r .emo("nerd")
Let's test the function under fixed parameter values:
log_lik_b_prod(param = c(50, 1.5, 2, 0.5), data = Aliens) ## test functions above
fit_constr <- fitme(y ~ 0 + (1|group) + offset(50 + 1.5*x), data = Aliens,
fixed = list(lambda = 2, phi = 0.5))
logLik(fit_constr)
Fitting the mixed model numerically r .emo("nerd")
bad_mod <- lm(y ~ x + group, data = Aliens)
(init_values <- c(bad_mod$coefficients[1], bad_mod$coefficients[2],
var.group = var(bad_mod$coefficients[-c(1:2)]),
var.error = deviance(bad_mod) / bad_mod$df.residual))
system.time(
opt <- nloptr::nloptr(x0 = init_values, eval_f = log_lik_b_prod, data = Aliens, scale = -1,
lb = 0.5*init_values, ub = 2*init_values,
opts = list(algorithm = "NLOPT_LN_BOBYQA", xtol_rel = 1.0e-4, maxeval = -1))
)
Fitting the mixed model numerically r .emo("proof")
fit_spaMM_ML <- fitme(y ~ x + (1|group), data = Aliens, method = "ML")
estimates <- rbind(opt$solution, as.numeric(c(fit_spaMM_ML$fixef, fit_spaMM_ML$lambda, fit_spaMM_ML$phi)))
colnames(estimates) <- c("intercept", "slope", "var.group", "var.error")
rownames(estimates) <- c("numeric", "spaMM")
estimates
c(logLik.num = -1 * opt$objective, logLik.spaMM = logLik(fit_spaMM_ML)[[1]])
The estimation of random effects r .emo("nerd")
We estimate the realized values of the random variable:
ranef_num <- sapply(levels(Aliens$group), function(group) {
nloptr::nloptr(x0 = 0, lik_b, level = group, intercept = opt$solution[1],
slope = opt$solution[2], var.rand = opt$solution[3],
var.error = opt$solution[4], data = Aliens, scale = -1,
lb = -4*sqrt(opt$solution[3]), ub = 4*sqrt(opt$solution[3]),
opts = list(algorithm = "NLOPT_LN_BOBYQA", xtol_rel = 1.0e-4, maxeval = -1))$solution})
rbind(ranef_num,
ranef_spaMM = ranef(fit_spaMM_ML)[1][[1]])
Best Linear Unbiased Predictions (BLUPs) r .emo("info")
- The numbers we obtained with
ranef() are BLUPs of the random effects.
- Even if BLUPs can be estimated, they are not parameters of the statistical model!
- BLUPs can be computed using ML or REML (REML is best when BLUPs are not added to fitted values, unclear otherwise)
curve(dnorm(x, mean = 0, sd = sqrt(estimates[1, "var.group"])), -5, 5, ylab = "density", xlab = "b")
points(dnorm(ranef_num, mean = 0, sd = sqrt(estimates[1, "var.group"])) ~ ranef_num, col = "blue", type = "h")
points(dnorm(ranef_num, mean = 0, sd = sqrt(estimates[1, "var.group"])) ~ ranef_num, col = "blue")
A simple example of LMM
The oatsyield dataset r .emo("alien")
head(oatsyield)
str(oatsyield)
The oatsyield dataset r .emo("alien")
coplot(yield ~ nitro | Variety + Block, data = oatsyield, type = "b")
LMM using {lme4} r .emo("practice")
(fit_lmm_lme4 <- lmer(yield ~ nitro + Variety + (1|Block), data = oatsyield))
LMM using {lme4} r .emo("practice")
head(model.matrix(fit_lmm_lme4)) ## X
fixef(fit_lmm_lme4) ## beta estimates
print(VarCorr(fit_lmm_lme4), comp = "Variance") ## extract REML variance estimates
LMM using {lme4} r .emo("practice")
head(getME(fit_lmm_lme4, "Z")) ## Z
getME(fit_lmm_lme4, "Zt") ## Z transposed
LMM using {lme4} r .emo("practice")
ranef(fit_lmm_lme4) ## b estimates
LMM using {spaMM} r .emo("practice")
(fit_lmm_spaMM <- fitme(yield ~ nitro + Variety + (1|Block), data = oatsyield, method = "REML"))
LMM using {spaMM} r .emo("practice")
head(model.matrix(fit_lmm_spaMM)) ## X
fixef(fit_lmm_spaMM) ## beta estimates
VarCorr(fit_lmm_spaMM) ## lambda and phi estimates
LMM using {spaMM} r .emo("practice")
head(get_ZALMatrix(fit_lmm_spaMM)) ## Z
LMM using {spaMM} r .emo("practice")
ranef(fit_lmm_spaMM) ## b estimates
Predictions
Predictions with LM r .emo("recap")
fit_lm <- lm(yield ~ nitro + Variety + Block, data = oatsyield)
data.for.pred <- data.frame(nitro = 0.3, Variety = "Victory", Block = levels(oatsyield$Block))
(p <- predict(fit_lm, newdata = data.for.pred, interval = "confidence", se.fit = TRUE))
Predictions with LM r .emo("recap")
library(ggplot2)
ggplot() +
geom_jitter(aes(y = yield, x = Block), data = oatsyield, shape = 1, width = 0.05, color = "blue") +
geom_point(aes(y = fit, x = Block), data = cbind(data.for.pred, p$fit)) +
geom_errorbar(aes(ymin = lwr, ymax = upr, x = Block), data = cbind(data.for.pred, p$fit),
width = 0.1) + theme_classic()
Marginal predictions for LMM using {lme4} r .emo("practice")
In LMM (but not in GLMM), you can obtain predictions averaged over the random variable (i.e. marginal prediction), by considering a realisation of the random effects equals to 0:
data.for.pred <- data.frame(nitro = 0.3, Variety = "Victory", Block = "new")
p <- predict(fit_lmm_lme4, newdata = data.for.pred, re.form = NA) ## NA to ignore ranef
X <- matrix(c(1, 0.3, 0, 1), nrow = 1) ## column of X corresponding to newdata
se.fit <- sqrt(X %*% vcov(fit_lmm_lme4) %*% t(X))
se.rand <- attr(VarCorr(fit_lmm_lme4)$Block, "stddev")
se.predVar <- as.numeric(sqrt(se.fit^2 + se.rand^2))
lwr <- as.numeric(p + qnorm(0.025) * se.predVar)
upr <- as.numeric(p + qnorm(0.975) * se.predVar)
c(fit = as.numeric(p), lwr = lwr, upr = upr)
Notes:
- {lme4} does not compute CI (because method shown here is not perfect).
- I don't know why people don't use the Student's here, but it would not make a big difference.
Marginal predictions for LMM using {spaMM} r .emo("practice")
In LMM (but not in GLMM), you can obtain predictions averaged over the random variable (i.e. marginal prediction), by considering a realisation of the random effects equals to 0:
data.for.pred <- data.frame(nitro = 0.3, Variety = "Victory", Block = "new")
predict(fit_lmm_spaMM, newdata = data.for.pred)
get_intervals(fit_lmm_spaMM, newdata = data.for.pred, intervals = "predVar")
Notes:
- {spaMM} directly produces CI as we computed them by hand for {lme4}.
- we use
intervals = "predVar" and not intervals = "fixefVar" to account for the uncertainty in both fixed and random effects.
Conditional predictions using {lme4} r .emo("practice")
Predictions conditional on the random variable ($\mathbf{X}\widehat{\beta}+\mathbf{Z}\widehat{b}$):
data.for.pred.with.b <- data.frame(nitro = 0.3, Variety = "Victory", Block = levels(oatsyield$Block))
library(merTools)
predictInterval(fit_lmm_lme4, newdata = data.for.pred.with.b, include.resid.var = FALSE, seed = 1)
Notes:
- using the companion package {merTools} makes things easier (but it cannot be used in the marginal case as it neglects random effect variance).
merTools::predictInterval() uses simulations to build CI, so results differ slightly each time you run the function (or set seed).
Conditional predictions using {spaMM} r .emo("practice")
Predictions conditional on the random variable ($\mathbf{X}\widehat{\beta}+\mathbf{Z}\widehat{b}$):
data.for.pred.with.b <- data.frame(nitro = 0.3, Variety = "Victory", Block = levels(oatsyield$Block))
predict(fit_lmm_spaMM, newdata = data.for.pred.with.b)
get_intervals(fit_lmm_spaMM, newdata = data.for.pred.with.b, intervals = "predVar")
Conditional predictions from LM and LMM r .emo("info")
p <- predict(fit_lm, newdata = data.for.pred.with.b, interval = "confidence", se.fit = TRUE)
p3 <- predict(fit_lmm_spaMM, newdata = data.for.pred.with.b, intervals = "predVar")
p4 <- predictInterval(fit_lmm_lme4, newdata = data.for.pred.with.b, include.resid.var = FALSE)
a <- cbind(data.for.pred.with.b, p$fit, method = "LM")
b <- cbind(data.for.pred.with.b, fit = p3[ ,1],
get_intervals(fit_lmm_spaMM, newdata = data.for.pred.with.b, intervals = "predVar"),
method = "spaMM")
colnames(b)[colnames(b) == "predVar_0.025"] <- "lwr"
colnames(b)[colnames(b) == "predVar_0.975"] <- "upr"
c <- cbind(data.for.pred.with.b, p4[, c(1, 3, 2)], method = "merTools")
res <- rbind(a, b, c)
ggplot() +
geom_jitter(aes(y = yield, x = Block), data = oatsyield, shape = 1, width = 0.05, color = "blue") +
geom_point(aes(y = fit, x = Block, colour = method), data = res,
position = position_dodge(width = 0.2)) +
geom_errorbar(aes(ymin = lwr, ymax = upr, x = Block, colour = method), data = res,
width = 0.1, position = position_dodge(width = 0.2)) +
geom_hline(yintercept = mean(tapply(oatsyield$yield, oatsyield$Block, mean)),
linetype = "dashed") +
scale_colour_manual(values = c("red", "purple", "orange")) +
theme_classic() +
theme(legend.position = "top")
Note: in LMM predictions (with Gaussian random effects) are always attracted toward the mean.
Tests
Testing the effect of nitro r .emo("practice")
With {lme4} using an asymptotic LRT:
fit_lmm_lme4_nonitro <- lmer(yield ~ Variety + (1|Block), data = oatsyield)
anova(fit_lmm_lme4, fit_lmm_lme4_nonitro)
Note: always use ML fits for testing fixed effects (but anova() from {lme4} does that automatically!)
Testing the effect of the effect of nitro r .emo("practice")
With {spaMM} using an asymptotic LRT:
fit_lmm_spaMM_ML <- fitme(yield ~ nitro + Variety + (1|Block), data = oatsyield, method = "ML")
fit_lmm_spaMM_nonitro <- fitme(yield ~ Variety + (1|Block), data = oatsyield, method = "ML")
anova(fit_lmm_spaMM_ML, fit_lmm_spaMM_nonitro)
Note: always use ML fits for testing fixed effects (anova() from {spaMM} does NOT do that automatically)
Reliability of the test of the asymptotic LRT r .emo("proof")
pvalues <- replicate(1000, {
oatsyield$yield <- simulate(fit_lmm_lme4_nonitro)[, 1]
fit_lmm_lme4_new <- lmer(yield ~ nitro + Variety + (1|Block), data = oatsyield, REML = FALSE)
## We use ML to save time because anova won't trigger refit
fit_lmm_lme4_new_nonitro <- lmer(yield ~ Variety + (1|Block), data = oatsyield, REML = FALSE)
anova(fit_lmm_lme4_new, fit_lmm_lme4_new_nonitro)$"Pr(>Chisq)"[2]})
plot(ecdf(pvalues), xlim = c(0, 1), ylim = c(0, 1)) ## looks good here, but design is ideal
abline(0, 1, col = 2, lwd = 2, lty = 2)
Testing the effect of the effect of nitro r .emo("practice")
With {lme4} + {pbkrtest} using parametric bootstrap:
pbkrtest::PBmodcomp(fit_lmm_lme4, fit_lmm_lme4_nonitro, nsim = 999)
Testing the effect of the effect of nitro r .emo("practice")
With {spaMM} using parametric bootstrap:
anova(fit_lmm_spaMM_ML, fit_lmm_spaMM_nonitro, boot.repl = 999)
res <- anova(fit_lmm_spaMM_ML, fit_lmm_spaMM_nonitro, boot.repl = 999)
print(res)
Never use the outdated aov()! r .emo("warn")
pvalues <- replicate(1000, {
oatsyield$nitro <- runif(1:nrow(oatsyield)) ## simulate H0 for nitro effect
fit_aov_sim <- aov(yield ~ nitro + Variety + Error(Block), data = oatsyield)
summary(fit_aov_sim)[[2]][[1]][2, "Pr(>F)"]})
plot(ecdf(pvalues), xlim = c(0, 1), ylim = c(0, 1))
abline(0, 1, col = 2, lwd = 2, lty = 2)
Note:
- even the help file tells you that tests using
aov() are "not particularly sensible statistically"...
A summary table for mixed models? r .emo("practice")
{lme4} or {spaMM} do not report p-values in the summary table because t-tests and Wald tests have bad properties for MMs. Some alternatives have been proposed:
library(lmerTest) ## overwrite lmer!
fit_lmm_lme4_nitro_bis <- lmer(yield ~ nitro + Variety + (1|Block), data = oatsyield, REML = FALSE)
summary(fit_lmm_lme4_nitro_bis)$coefficients
detach(package:lmerTest) ## to restore original lmer
Note: here you must use a ML fit (not automatic here, and the results would change if you don't!).
The method relies on the so-called Satterthwaite's degrees of freedom method, which corrects the degrees of freedom to account for the facts that the design is not always balanced.
A summary table for mixed models? r .emo("practice")
You can modify options in {spaMM} to give you Wald p-values, but those are unlikely to be reliable:
summary(fit_lmm_spaMM_ML, details = c(p_value = "Wald"), verbose = FALSE)$beta
Instead, it should be better to compute confidence intervals by parametric bootstrap:
boot_res <- confint(fit_lmm_spaMM_ML, parm = "nitro", boot_args = list(nsim = 1000, seed = 123), verbose = FALSE)
boot_res$normal
Testing random effects r .emo("practice")
It is possible, at least for relatively simple LMM to test if the variance of the random effect is null:
fit1 <- lmer(yield ~ nitro + Variety + (1|Block), data = oatsyield, REML = FALSE)
fit2 <- lm(yield ~ nitro + Variety, data = oatsyield)
RLRsim::exactLRT(fit1, fit2)
Testing random effects r .emo("practice")
The same package can be used with model fitted with {spaMM} with some more code:
fit1 <- fitme(yield ~ nitro + Variety + (1|Block), data = oatsyield)
fit2 <- fitme(yield ~ nitro + Variety, data = oatsyield)
(obs.LRT <- 2*(logLik(fit1) - logLik(fit2)))
args <- get_RLRsim_args(fit1, fit2)
sim.LRT <- do.call(RLRsim::LRTSim, args)
(RLRpval <- (sum(sim.LRT >= obs.LRT) + 1) / (length(sim.LRT) + 1))
Note: you can alternatively use anova(fit1, fit2, boot.repl = 999) but it is much slower.
Information Criteria
Information Criteria r .emo("info")
As for predictions, you can compute a marginal or a conditional AIC depending on whether you are interested in the fit under the average realisations of the random effect or conditional on the estimates for such random values under the observed levels in particular.
AIC(fit_lmm_spaMM)
print(AIC(fit_lmm_spaMM))
- rely on the marginal AIC if you are interested in marginal predictions.
- rely on the conditional AIC if you are interested in conditional predictions.
{lme4} only returns the marginal AIC, but (again) you can use another companion package for computing the conditional AIC (note: methods implemented in {spaMM} and {cAIC4} are not the same, so results may differ):
cAIC4::cAIC(fit_lmm_lme4)
Assumptions
What are the assumptions in LMM? r .emo("info")
Same as LM (except for independence) +
- given the fixed effects and the realized values of the random effects, the errors must be independent
- the random effects must follow the assumed distribution
- the levels of the factorial variable for which random effects are estimated must be representative from the whole population
Plotting residuals with {lme4} r .emo("practice")
plot(fit_lmm_lme4, type = c("p", "smooth")) ## see ?lme4:::plot.merMod for details
Plotting residuals with {lme4} r .emo("practice")
plot(fit_lmm_lme4, resid(., scaled=TRUE) ~ fitted(.) | Block, abline = 0)
Plotting residuals with {lme4} r .emo("practice")
lattice::qqmath(fit_lmm_lme4, id = 0.05) ## id allows to see outliers
Plotting BLUPs with {lme4} r .emo("practice")
lattice::qqmath(lme4::ranef(fit_lmm_lme4))
Plotting residuals with {spaMM} r .emo("practice")
plot(fit_lmm_spaMM)
Note: it does not work on the slide, but should work in your RStudio...
Using simulated residuals r .emo("practice")
You can also rely on {DHARMa} as we have seen it for GLM:
library(DHARMa)
plot(simulateResiduals(fit_lmm_spaMM, n = 1000))
Note: {DHARMa} works for both {spaMM} and {lme4}!
Random or fixed?
Choosing between fixed and random effects r .emo("info")
The choice between considering a predictor as a fixed or random effect can be difficult; there is a trade-off between both approaches.
Fixed
- no assumption about the distribution of the parameter values associated with each level of a predictor
- requires many datapoints for each additional level to get reliable results
- allows for the prediction of the effect of observed levels only (for factors)
- simple to handle
Random
- the values associated with the levels of a predictor follow a probability distribution (usually Gaussian)
- requires at least one datapoint for each additional level to get reliable results (more is better)
- requires many levels for the variance estimates to be reliable (5 being strict minimum)
- allow for the prediction of the effect of both observed and unobserved levels (for factors)
- more difficult to handle
Specifying multiple random effects
The lme4::Penicillin dataset r .emo("alien")
str(Penicillin)
table(Penicillin$sample, Penicillin$plate)
The lme4::Penicillin dataset r .emo("info")
The random effects are crossed, i.e. we have now 2 independent random effects:
fit_peni <- fitme(diameter ~ 1 + (1|plate) + (1|sample), data = Penicillin)
fit_peni$lambda
The lme4::cake dataset r .emo("alien")
head(cake)
str(cake)
The lme4::cake dataset r .emo("alien")
table(cake$recipe, cake$replicate, cake$temperature)
The lme4::cake dataset r .emo("info")
The random effect is nested within a fixed effect, so we must account for that:
fit_cake <- fitme(angle ~ recipe + temperature + (1|recipe:replicate), data = cake)
fit_cake$lambda
Note:
- considering random effects as crossed here would be wrong. For example, it would imply that the random realization of the effect associated with replicate 1 from recipe A would be the same as that for the replicate 1 from recipe B.
The lme4::cake dataset r .emo("practice")
The random effect is nested within a fixed effect, so we must account for that, (alternative specification):
cake$replicate_tot <- factor(paste(cake$recipe, cake$replicate, sep = "_"))
levels(cake$replicate_tot)
fit_cake <- fitme(angle ~ recipe + temperature + (1|replicate_tot), data = cake)
fit_cake$lambda
Note:
- whatever you do, the important thing is to make sure that replicates across recipes do not share the same names.
The carnivora dataset r .emo("alien")
carnivora$log_brain <- log(carnivora$Brain)
carnivora$log_body <- log(carnivora$Weight)
str(carnivora)
The carnivora dataset r .emo("info")
Genus names are unique across families:
length(unique(paste(carnivora$Genus, carnivora$Family))) == length(unique(carnivora$Genus))
coplot(log_brain ~ log_body | Family, data = carnivora)
The carnivora dataset r .emo("info")
Here we could consider one random effect nested into another one:
fit_carnivora <- fitme(log_brain ~ log_body + (1|Family/Genus), data = carnivora, method = "REML")
fit_carnivora
The carnivora dataset r .emo("info")
Here we could consider one random effect nested into another one (alternative specification):
fit_carnivora_bis <- fitme(log_brain ~ log_body + (1|Family) + (1|Family:Genus), data = carnivora, method = "REML")
fit_carnivora_bis
The carnivora dataset r .emo("info")
But since genus names are not (here) replicated across families, nested random effects are equivalent to crossed random effects:
fit_carnivora_ter <- fitme(log_brain ~ log_body + (1|Family) + (1|Genus), data = carnivora, method = "REML")
fit_carnivora_ter
Checking the random structure r .emo("practice")
Checking the names of the BLUPs structure is an easy solution to check that you did specify the random structure correctly:
lapply(ranef(fit_carnivora), function(x) head(names(x))) ## or simply ranef(fit_carnivora)
Random slopes
Fitting a random slope model r .emo("practice")
(fit_rand_slope_spaMM <- fitme(log_brain ~ log_body + (log_body|Family) + (1|Genus),
data = carnivora, method = "REML"))
The BLUPs for the slopes r .emo("info")
A random slope model is a model in which one slope is a random variable.
Here the effect of log_body is considered to vary across families:
ranef(fit_rand_slope_spaMM)$`( log_body | Family )`
Predictions r .emo("practice")
data_pred <- predict(fit_rand_slope_spaMM,
newdata = expand.grid(log_body = range(carnivora$log_body),
Family = levels(carnivora$Family), Genus = "new"),
binding = "log_brain_pred") ## Tip: binding adds predictions to newdata
ggplot() +
geom_line(aes(y = log_brain_pred, x = log_body, colour = Family), data = data_pred) +
geom_point(aes(y = log_brain, x = log_body, colour = Family), shape = 1, data = carnivora) +
theme_minimal()
Testing the random slope using {lme4} r .emo("practice")
Asymptotic LRT are really poor for that, so better use parametric bootstrap:
fit_rand_slope_lme4 <- lmer(log_brain ~ log_body + (log_body|Family) + (1|Genus), data = carnivora)
fit_no_rand_slope_lme4 <- lmer(log_brain ~ log_body + (1|Family) + (1|Genus), data = carnivora)
pbkrtest::PBmodcomp(fit_rand_slope_lme4, fit_no_rand_slope_lme4)
Notes:
- fits differ by 2 degrees of freedom (1 for the fit of the variance of the random slope + 1 for its covariance with the random intercept).
pbkrtest::PBmodcomp() uses ML logLik.
Testing the random slope using {spaMM} r .emo("practice")
Asymptotic LRT are really poor for that, so better use parametric bootstrap r .emo("slow"):
fit_no_rand_slope_spaMM <- fitme(log_brain ~ log_body + (1|Family) + (1|Genus),
data = carnivora, method = "REML")
anova(fit_rand_slope_spaMM, fit_no_rand_slope_spaMM, boot.repl = 999)
Note:
- I used REML fit unlike
pbkrtest::PBmodcomp() (but it is not clear what the best choice is).
Generalized Linear Mixed-effects Models
GLM + LMM = GLMM r .emo("info")
$$\begin{array}{lcl}
\mu &=& g^{-1}(\eta)\
\mu &=& g^{-1}(\mathbf{X}\beta + \mathbf{Z}b)\
\end{array}
$$
with (as for GLM):
- $\text{E}(\text{Y}) = \mu = g^{-1}(\eta)$
- $\text{Var}(\text{Y}) = \phi\text{V}(\mu)$
Notes:
- If $g^{-1}$ is the identity function, $\phi = \sigma^2$ and $\text{V}(\mu) = 1$, we have the LMM.
- If $\mathbf{Z}b = 0$, we have the GLM.
- If $g^{-1}$ is the identity function, $\phi = \sigma^2$, $\text{V}(\mu) = 1$, and $\mathbf{Z}b = 0$, we have the LM.
If you combine everything we said about GLMs and LMMs, then you should know how to handle GLMMs.
Let illustrate this with an example!
The Flatwork dataset r .emo("alien")
Flatwork
The Flatwork dataset r .emo("alien")
str(Flatwork)
Fitting GLMM with {lme4} r .emo("practice")
(fit_fw_lme4_pois <- glmer(shopping ~ gender + (1|individual) + (1|month), family = poisson(),
data = Flatwork))
Note:
- with {lme4} you cannot choose between ML and REML for GLMMs (it does a ML fit).
Fitting GLMM with {spaMM} r .emo("practice")
With {spaMM} you can still choose between ML and REML for GLMMs (REML sensu stricto does not exist for GLMMs, but {spaMM} has implemented a generalisation of the concept proposed in the statistical literature).
By default, as for LMMs, spaMM::fitme() produces a ML fit.
(fit_fw_spaMM_pois <- fitme(shopping ~ gender + (1|individual) + (1|month), family = poisson(),
data = Flatwork))
Checking residuals r .emo("practice")
You can still test the assumptions of GLMMs with {DHARMa}:
fw_resid_spaMM_pois <- simulateResiduals(fit_fw_spaMM_pois)
plot(fw_resid_spaMM_pois)
Overdispersion? r .emo("practice")
testDispersion(fw_resid_spaMM_pois)
Extra 0s? r .emo("practice")
testZeroInflation(fw_resid_spaMM_pois)
Alternative 1: negative binomial r .emo("practice")
(fit_fw_spaMM_nb <- fitme(shopping ~ gender + (1|individual) + (1|month), family = spaMM::negbin(),
data = Flatwork))
Alternative 1: negative binomial r .emo("practice")
fw_resid_spaMM_nb <- simulateResiduals(fit_fw_spaMM_nb)
plot(fw_resid_spaMM_nb)
Alternative 1: negative binomial r .emo("practice")
testDispersion(fw_resid_spaMM_nb)
Alternative 1: negative binomial r .emo("practice")
testZeroInflation(fw_resid_spaMM_nb)
Alternative 2: zero-augmented Poisson r .emo("practice")
library(glmmTMB)
fit_fw_glmmTMB_zipois <- glmmTMB(shopping ~ gender + (1|individual) + (1|month), family = poisson(),
data = Flatwork, ziformula = ~ 1)
Like {spaMM}, {glmmTMB} is extremely flexible package (it can fit models {spaMM} cannot fit, and the reverse is also true).
{glmmTMB} allows for both zero-inflation and hurdle models (using truncated distributions), as well as many covariance functions, modelling the residual variance...
There is no full-fledged parametric bootstrap methods implemented as far as I know, nor conditional AIC, but this could change in the near future.
The package is a wrapper around general numerical methods, so it may not be as robust as packages build with mixed models in mind (e.g. {spaMM} and {lme4}) in some cases; so convergence issues could happen more often.
Alternative 2: zero-inflated Poisson r .emo("practice")
{glmmTMB} is again compatible with {DHARMa}:
fw_resid_glmmTMB_zipois <- simulateResiduals(fit_fw_glmmTMB_zipois)
plot(fw_resid_glmmTMB_zipois)
Alternative 2: zero-inflated Poisson r .emo("practice")
{glmmTMB} is also compatible with {DHARMa}:
testDispersion(fw_resid_glmmTMB_zipois)
Alternative 2: zero-inflated Poisson r .emo("practice")
{glmmTMB} is also compatible with {DHARMa}:
testZeroInflation(fw_resid_glmmTMB_zipois)
Alternative 3: zero-inflated negative binomial r .emo("practice")
(fit_fw_glmmTMB_zinb <- glmmTMB(shopping ~ gender + (1|individual) + (1|month), family = nbinom1(),
data = Flatwork, ziformula = ~ 1))
Alternative 3: zero-inflated negative binomial r .emo("practice")
fw_resid_glmmTMB_zinb <- simulateResiduals(fit_fw_glmmTMB_zinb)
plot(fw_resid_glmmTMB_zinb)
Alternative 3: zero-inflated negative binomial r .emo("practice")
testDispersion(fw_resid_glmmTMB_zinb)
Alternative 3: zero-inflated negative binomial r .emo("practice")
testZeroInflation(fw_resid_glmmTMB_zinb)
What is the best alternative for our GLMM? r .emo("practice")
AIC(fit_fw_spaMM_pois)
print(AIC(fit_fw_spaMM_pois))
AIC(fit_fw_spaMM_nb)
print(AIC(fit_fw_spaMM_nb))
AIC(fit_fw_glmmTMB_zipois, fit_fw_glmmTMB_zinb) #Marginal AIC
Note:
- the zero-inflated negative binomial model fits the data much better than the alternative models we tried.
Predictions for GLMM? r .emo("info")
-
building conditional predictions is similar to what we have seen in GLM and LMM: just use predict() (or fitme::pdep_effects()) for that.
-
building marginal predictions is more complex because you need to integrate the predictions over the distribution of the random effects... There is no function doing that automatically and this is a little too advanced for this course.
What you need to remember r .emo("goal")
- what mixed models (LMM & GLMM) are
- how to use {lme4}, {spaMM} and {glmmTMB} to fit mixed models
- how to make predictions, tests, CI...
- how to check assumptions
- when to consider effects as fixed or random
- how to implement different random effects structure
Table of contents
Mixed-effects models
- 4.0 Introduction to LMM & GLMM
- 4.1 Solving LM problems using LMM
- 4.2 A showcase of some useful applications
[Back to main menu](./Title.html#2)
courtiol/LM2GLMM documentation built on July 3, 2022, 7:42 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
library(LM2GLMM) library(doSNOW) spaMM::spaMM.options(nb_cores = parallel::detectCores() - 1) options(width = 120) knitr::opts_chunk$set(cache = TRUE, cache.path = "./cache_knitr/MM_intro/", fig.path = "./fig_knitr/MM_intro/", fig.width = 5, fig.height = 5, fig.align = "center")
Mixed-effects models
- 4.0 Introduction to LMM & GLMM
- 4.1 Solving LM problems using LMM
- 4.2 A showcase of some useful applications
[Back to main menu](./Title.html#2)
You will learn in this session r .emo("goal")
- what mixed models (LMM & GLMM) are
- how to use {lme4}, {spaMM} and {glmmTMB} to fit mixed models
- how to make predictions, tests, CI...
- how to check assumptions
- when to consider effects as fixed or random
- how to implement different random effects structure
Linear Mixed-effects Models
Why mixed models? r .emo("info")
Goal:
To study, or to account for, unobservable sources of heterogeneity between observations.
Mixed-effects models allow for:
- the study of other questions than LM/GLM (e.g. heritability)
- the fixing of assumption violations in LM/GLM (lack of independence, some cases of overdispersion)
- the reduction of the uncertainty in estimates and predictions in cases where many parameters would have to be estimated at the cost of an additional hypothesis (the distribution of the random effects)
The main sources of heterogeneity considered by mixed-effects models are:
- origin (in its widest sense)
- time
- space
The linear mixed-effects model r .emo("info")
Definition
- a LMM is a specific linear model for which, given the design matrix $\mathbf{X}$, the responses ($\mathbf{Y}$) are no longer independent, but where the correlations can be described in terms of a random effect, i.e. a random variable that is not included in the predictor variables
Mathematical notation of LM r .emo("info")
LM, which we have seen before:
Mathematical notation of LMM r .emo("info")
LMM with one random factor with $q$ levels:
with $\text{E}(b) = 0$, $\text{Cov}(b) = \mathbf{C}$, which is symmetrical ($c_{i, j} = c_{j, i}$). Also, $\text{Cov}(b, \epsilon) = 0$.
Mathematical notation of LMM r .emo("info")
LMM with one random factor with $q$ levels:
Note: the dimensions of the matrices differ: $n \times n$ for $\mathbf{X}$ vs $q \times q$ for $\mathbf{Z}$
Fitting procedure
A simple simulation function r .emo("alien")
simulate_Mix <- function(intercept, slope, n, group_nb, var.rand, var.error){ data <- data.frame(intercept = intercept, slope = slope, x = runif(n)) group_compo <- rmultinom(n = 1, size = n, prob = c(rep(1/group_nb, group_nb))) data$group <- factor(rep(paste("group", 1:group_nb, sep = "_"), group_compo)) data$b <- rep(rnorm(group_nb, mean = 0, sd = sqrt(var.rand)), group_compo) data$error <- rnorm(n, mean = 0, sd = sqrt(var.error)) data$y <- data$intercept + data$slope*data$x + data$b + data$error return(data) } set.seed(1) Aliens <- simulate_Mix(intercept = 50, slope = 1.5, n = 30, group_nb = 10, var.rand = 2, var.error = 0.5)
Note: look carefully at the next slide to better understand.
Our toy dataset r .emo("alien")
Aliens
Fitting the model with {lme4} r .emo("practice")
library(lme4) (fit_lme4 <- lmer(y ~ x + (1|group), data = Aliens)) ## REML fit by default! print(VarCorr(fit_lme4), comp = "Variance") ## extract REML variance estimates
Fitting the model with {spaMM} r .emo("practice")
library(spaMM) (fit_spaMM <- fitme(y ~ x + (1|group), data = Aliens, method = "REML"))
Functions for fitting the mixed model numerically r .emo("nerd")
lik_b <- function(b.vec, level, intercept, slope, var.rand, var.error, data, scale = 1) { lik <- sapply(b.vec, function(b){ sub_data <- data[which(data$group == level), ] sub_data$pred <- intercept + slope*sub_data$x + b sub_data$conditional.density <- dnorm(sub_data$y, mean = sub_data$pred, sd = sqrt(var.error)) return(dnorm(b, mean = 0, sd = sqrt(var.rand)) * prod(sub_data$conditional.density)) }) return(scale * lik) }
log_lik_b_prod <- function(param, data, scale = 1){ log_lik_vec <- sapply(levels(data$group), function(level) { log(integrate(lik_b, -Inf, Inf, level = level, intercept = param[1], slope = param[2], var.rand = param[3], var.error = param[4], data = data)$value)}) return(scale * sum(log_lik_vec)) }
Note: this is a ML fit for simplicity.
Testing the functions r .emo("nerd")
Let's test the function under fixed parameter values:
log_lik_b_prod(param = c(50, 1.5, 2, 0.5), data = Aliens) ## test functions above fit_constr <- fitme(y ~ 0 + (1|group) + offset(50 + 1.5*x), data = Aliens, fixed = list(lambda = 2, phi = 0.5)) logLik(fit_constr)
Fitting the mixed model numerically r .emo("nerd")
bad_mod <- lm(y ~ x + group, data = Aliens)
(init_values <- c(bad_mod$coefficients[1], bad_mod$coefficients[2], var.group = var(bad_mod$coefficients[-c(1:2)]), var.error = deviance(bad_mod) / bad_mod$df.residual))
system.time( opt <- nloptr::nloptr(x0 = init_values, eval_f = log_lik_b_prod, data = Aliens, scale = -1, lb = 0.5*init_values, ub = 2*init_values, opts = list(algorithm = "NLOPT_LN_BOBYQA", xtol_rel = 1.0e-4, maxeval = -1)) )
Fitting the mixed model numerically r .emo("proof")
fit_spaMM_ML <- fitme(y ~ x + (1|group), data = Aliens, method = "ML") estimates <- rbind(opt$solution, as.numeric(c(fit_spaMM_ML$fixef, fit_spaMM_ML$lambda, fit_spaMM_ML$phi))) colnames(estimates) <- c("intercept", "slope", "var.group", "var.error") rownames(estimates) <- c("numeric", "spaMM") estimates c(logLik.num = -1 * opt$objective, logLik.spaMM = logLik(fit_spaMM_ML)[[1]])
The estimation of random effects r .emo("nerd")
We estimate the realized values of the random variable:
ranef_num <- sapply(levels(Aliens$group), function(group) { nloptr::nloptr(x0 = 0, lik_b, level = group, intercept = opt$solution[1], slope = opt$solution[2], var.rand = opt$solution[3], var.error = opt$solution[4], data = Aliens, scale = -1, lb = -4*sqrt(opt$solution[3]), ub = 4*sqrt(opt$solution[3]), opts = list(algorithm = "NLOPT_LN_BOBYQA", xtol_rel = 1.0e-4, maxeval = -1))$solution}) rbind(ranef_num, ranef_spaMM = ranef(fit_spaMM_ML)[1][[1]])
Best Linear Unbiased Predictions (BLUPs) r .emo("info")
- The numbers we obtained with
ranef()are BLUPs of the random effects. - Even if BLUPs can be estimated, they are not parameters of the statistical model!
- BLUPs can be computed using ML or REML (REML is best when BLUPs are not added to fitted values, unclear otherwise)
curve(dnorm(x, mean = 0, sd = sqrt(estimates[1, "var.group"])), -5, 5, ylab = "density", xlab = "b") points(dnorm(ranef_num, mean = 0, sd = sqrt(estimates[1, "var.group"])) ~ ranef_num, col = "blue", type = "h") points(dnorm(ranef_num, mean = 0, sd = sqrt(estimates[1, "var.group"])) ~ ranef_num, col = "blue")
A simple example of LMM
The oatsyield dataset r .emo("alien")
head(oatsyield) str(oatsyield)
The oatsyield dataset r .emo("alien")
coplot(yield ~ nitro | Variety + Block, data = oatsyield, type = "b")
LMM using {lme4} r .emo("practice")
(fit_lmm_lme4 <- lmer(yield ~ nitro + Variety + (1|Block), data = oatsyield))
LMM using {lme4} r .emo("practice")
head(model.matrix(fit_lmm_lme4)) ## X fixef(fit_lmm_lme4) ## beta estimates print(VarCorr(fit_lmm_lme4), comp = "Variance") ## extract REML variance estimates
LMM using {lme4} r .emo("practice")
head(getME(fit_lmm_lme4, "Z")) ## Z getME(fit_lmm_lme4, "Zt") ## Z transposed
LMM using {lme4} r .emo("practice")
ranef(fit_lmm_lme4) ## b estimates
LMM using {spaMM} r .emo("practice")
(fit_lmm_spaMM <- fitme(yield ~ nitro + Variety + (1|Block), data = oatsyield, method = "REML"))
LMM using {spaMM} r .emo("practice")
head(model.matrix(fit_lmm_spaMM)) ## X fixef(fit_lmm_spaMM) ## beta estimates VarCorr(fit_lmm_spaMM) ## lambda and phi estimates
LMM using {spaMM} r .emo("practice")
head(get_ZALMatrix(fit_lmm_spaMM)) ## Z
LMM using {spaMM} r .emo("practice")
ranef(fit_lmm_spaMM) ## b estimates
Predictions
Predictions with LM r .emo("recap")
fit_lm <- lm(yield ~ nitro + Variety + Block, data = oatsyield) data.for.pred <- data.frame(nitro = 0.3, Variety = "Victory", Block = levels(oatsyield$Block)) (p <- predict(fit_lm, newdata = data.for.pred, interval = "confidence", se.fit = TRUE))
Predictions with LM r .emo("recap")
library(ggplot2) ggplot() + geom_jitter(aes(y = yield, x = Block), data = oatsyield, shape = 1, width = 0.05, color = "blue") + geom_point(aes(y = fit, x = Block), data = cbind(data.for.pred, p$fit)) + geom_errorbar(aes(ymin = lwr, ymax = upr, x = Block), data = cbind(data.for.pred, p$fit), width = 0.1) + theme_classic()
Marginal predictions for LMM using {lme4} r .emo("practice")
In LMM (but not in GLMM), you can obtain predictions averaged over the random variable (i.e. marginal prediction), by considering a realisation of the random effects equals to 0:
data.for.pred <- data.frame(nitro = 0.3, Variety = "Victory", Block = "new") p <- predict(fit_lmm_lme4, newdata = data.for.pred, re.form = NA) ## NA to ignore ranef X <- matrix(c(1, 0.3, 0, 1), nrow = 1) ## column of X corresponding to newdata se.fit <- sqrt(X %*% vcov(fit_lmm_lme4) %*% t(X)) se.rand <- attr(VarCorr(fit_lmm_lme4)$Block, "stddev") se.predVar <- as.numeric(sqrt(se.fit^2 + se.rand^2)) lwr <- as.numeric(p + qnorm(0.025) * se.predVar) upr <- as.numeric(p + qnorm(0.975) * se.predVar) c(fit = as.numeric(p), lwr = lwr, upr = upr)
Notes:
- {lme4} does not compute CI (because method shown here is not perfect).
- I don't know why people don't use the Student's here, but it would not make a big difference.
Marginal predictions for LMM using {spaMM} r .emo("practice")
In LMM (but not in GLMM), you can obtain predictions averaged over the random variable (i.e. marginal prediction), by considering a realisation of the random effects equals to 0:
data.for.pred <- data.frame(nitro = 0.3, Variety = "Victory", Block = "new") predict(fit_lmm_spaMM, newdata = data.for.pred) get_intervals(fit_lmm_spaMM, newdata = data.for.pred, intervals = "predVar")
Notes:
- {spaMM} directly produces CI as we computed them by hand for {lme4}.
- we use
intervals = "predVar"and notintervals = "fixefVar"to account for the uncertainty in both fixed and random effects.
Conditional predictions using {lme4} r .emo("practice")
Predictions conditional on the random variable ($\mathbf{X}\widehat{\beta}+\mathbf{Z}\widehat{b}$):
data.for.pred.with.b <- data.frame(nitro = 0.3, Variety = "Victory", Block = levels(oatsyield$Block)) library(merTools) predictInterval(fit_lmm_lme4, newdata = data.for.pred.with.b, include.resid.var = FALSE, seed = 1)
Notes:
- using the companion package {merTools} makes things easier (but it cannot be used in the marginal case as it neglects random effect variance).
merTools::predictInterval()uses simulations to build CI, so results differ slightly each time you run the function (or setseed).
Conditional predictions using {spaMM} r .emo("practice")
Predictions conditional on the random variable ($\mathbf{X}\widehat{\beta}+\mathbf{Z}\widehat{b}$):
data.for.pred.with.b <- data.frame(nitro = 0.3, Variety = "Victory", Block = levels(oatsyield$Block)) predict(fit_lmm_spaMM, newdata = data.for.pred.with.b) get_intervals(fit_lmm_spaMM, newdata = data.for.pred.with.b, intervals = "predVar")
Conditional predictions from LM and LMM r .emo("info")
p <- predict(fit_lm, newdata = data.for.pred.with.b, interval = "confidence", se.fit = TRUE) p3 <- predict(fit_lmm_spaMM, newdata = data.for.pred.with.b, intervals = "predVar") p4 <- predictInterval(fit_lmm_lme4, newdata = data.for.pred.with.b, include.resid.var = FALSE) a <- cbind(data.for.pred.with.b, p$fit, method = "LM") b <- cbind(data.for.pred.with.b, fit = p3[ ,1], get_intervals(fit_lmm_spaMM, newdata = data.for.pred.with.b, intervals = "predVar"), method = "spaMM") colnames(b)[colnames(b) == "predVar_0.025"] <- "lwr" colnames(b)[colnames(b) == "predVar_0.975"] <- "upr" c <- cbind(data.for.pred.with.b, p4[, c(1, 3, 2)], method = "merTools") res <- rbind(a, b, c) ggplot() + geom_jitter(aes(y = yield, x = Block), data = oatsyield, shape = 1, width = 0.05, color = "blue") + geom_point(aes(y = fit, x = Block, colour = method), data = res, position = position_dodge(width = 0.2)) + geom_errorbar(aes(ymin = lwr, ymax = upr, x = Block, colour = method), data = res, width = 0.1, position = position_dodge(width = 0.2)) + geom_hline(yintercept = mean(tapply(oatsyield$yield, oatsyield$Block, mean)), linetype = "dashed") + scale_colour_manual(values = c("red", "purple", "orange")) + theme_classic() + theme(legend.position = "top")
Note: in LMM predictions (with Gaussian random effects) are always attracted toward the mean.
Tests
Testing the effect of nitro r .emo("practice")
With {lme4} using an asymptotic LRT:
fit_lmm_lme4_nonitro <- lmer(yield ~ Variety + (1|Block), data = oatsyield) anova(fit_lmm_lme4, fit_lmm_lme4_nonitro)
Note: always use ML fits for testing fixed effects (but anova() from {lme4} does that automatically!)
Testing the effect of the effect of nitro r .emo("practice")
With {spaMM} using an asymptotic LRT:
fit_lmm_spaMM_ML <- fitme(yield ~ nitro + Variety + (1|Block), data = oatsyield, method = "ML") fit_lmm_spaMM_nonitro <- fitme(yield ~ Variety + (1|Block), data = oatsyield, method = "ML") anova(fit_lmm_spaMM_ML, fit_lmm_spaMM_nonitro)
Note: always use ML fits for testing fixed effects (anova() from {spaMM} does NOT do that automatically)
Reliability of the test of the asymptotic LRT r .emo("proof")
pvalues <- replicate(1000, { oatsyield$yield <- simulate(fit_lmm_lme4_nonitro)[, 1] fit_lmm_lme4_new <- lmer(yield ~ nitro + Variety + (1|Block), data = oatsyield, REML = FALSE) ## We use ML to save time because anova won't trigger refit fit_lmm_lme4_new_nonitro <- lmer(yield ~ Variety + (1|Block), data = oatsyield, REML = FALSE) anova(fit_lmm_lme4_new, fit_lmm_lme4_new_nonitro)$"Pr(>Chisq)"[2]}) plot(ecdf(pvalues), xlim = c(0, 1), ylim = c(0, 1)) ## looks good here, but design is ideal abline(0, 1, col = 2, lwd = 2, lty = 2)
Testing the effect of the effect of nitro r .emo("practice")
With {lme4} + {pbkrtest} using parametric bootstrap:
pbkrtest::PBmodcomp(fit_lmm_lme4, fit_lmm_lme4_nonitro, nsim = 999)
Testing the effect of the effect of nitro r .emo("practice")
With {spaMM} using parametric bootstrap:
anova(fit_lmm_spaMM_ML, fit_lmm_spaMM_nonitro, boot.repl = 999)
res <- anova(fit_lmm_spaMM_ML, fit_lmm_spaMM_nonitro, boot.repl = 999)
print(res)
Never use the outdated aov()! r .emo("warn")
pvalues <- replicate(1000, { oatsyield$nitro <- runif(1:nrow(oatsyield)) ## simulate H0 for nitro effect fit_aov_sim <- aov(yield ~ nitro + Variety + Error(Block), data = oatsyield) summary(fit_aov_sim)[[2]][[1]][2, "Pr(>F)"]}) plot(ecdf(pvalues), xlim = c(0, 1), ylim = c(0, 1)) abline(0, 1, col = 2, lwd = 2, lty = 2)
Note:
- even the help file tells you that tests using
aov()are "not particularly sensible statistically"...
A summary table for mixed models? r .emo("practice")
{lme4} or {spaMM} do not report p-values in the summary table because t-tests and Wald tests have bad properties for MMs. Some alternatives have been proposed:
library(lmerTest) ## overwrite lmer! fit_lmm_lme4_nitro_bis <- lmer(yield ~ nitro + Variety + (1|Block), data = oatsyield, REML = FALSE) summary(fit_lmm_lme4_nitro_bis)$coefficients detach(package:lmerTest) ## to restore original lmer
Note: here you must use a ML fit (not automatic here, and the results would change if you don't!).
The method relies on the so-called Satterthwaite's degrees of freedom method, which corrects the degrees of freedom to account for the facts that the design is not always balanced.
A summary table for mixed models? r .emo("practice")
You can modify options in {spaMM} to give you Wald p-values, but those are unlikely to be reliable:
summary(fit_lmm_spaMM_ML, details = c(p_value = "Wald"), verbose = FALSE)$beta
Instead, it should be better to compute confidence intervals by parametric bootstrap:
boot_res <- confint(fit_lmm_spaMM_ML, parm = "nitro", boot_args = list(nsim = 1000, seed = 123), verbose = FALSE) boot_res$normal
Testing random effects r .emo("practice")
It is possible, at least for relatively simple LMM to test if the variance of the random effect is null:
fit1 <- lmer(yield ~ nitro + Variety + (1|Block), data = oatsyield, REML = FALSE) fit2 <- lm(yield ~ nitro + Variety, data = oatsyield) RLRsim::exactLRT(fit1, fit2)
Testing random effects r .emo("practice")
The same package can be used with model fitted with {spaMM} with some more code:
fit1 <- fitme(yield ~ nitro + Variety + (1|Block), data = oatsyield) fit2 <- fitme(yield ~ nitro + Variety, data = oatsyield) (obs.LRT <- 2*(logLik(fit1) - logLik(fit2))) args <- get_RLRsim_args(fit1, fit2) sim.LRT <- do.call(RLRsim::LRTSim, args) (RLRpval <- (sum(sim.LRT >= obs.LRT) + 1) / (length(sim.LRT) + 1))
Note: you can alternatively use anova(fit1, fit2, boot.repl = 999) but it is much slower.
Information Criteria
Information Criteria r .emo("info")
As for predictions, you can compute a marginal or a conditional AIC depending on whether you are interested in the fit under the average realisations of the random effect or conditional on the estimates for such random values under the observed levels in particular.
AIC(fit_lmm_spaMM)
print(AIC(fit_lmm_spaMM))
- rely on the marginal AIC if you are interested in marginal predictions.
- rely on the conditional AIC if you are interested in conditional predictions.
{lme4} only returns the marginal AIC, but (again) you can use another companion package for computing the conditional AIC (note: methods implemented in {spaMM} and {cAIC4} are not the same, so results may differ):
cAIC4::cAIC(fit_lmm_lme4)
Assumptions
What are the assumptions in LMM? r .emo("info")
Same as LM (except for independence) +
- given the fixed effects and the realized values of the random effects, the errors must be independent
- the random effects must follow the assumed distribution
- the levels of the factorial variable for which random effects are estimated must be representative from the whole population
Plotting residuals with {lme4} r .emo("practice")
plot(fit_lmm_lme4, type = c("p", "smooth")) ## see ?lme4:::plot.merMod for details
Plotting residuals with {lme4} r .emo("practice")
plot(fit_lmm_lme4, resid(., scaled=TRUE) ~ fitted(.) | Block, abline = 0)
Plotting residuals with {lme4} r .emo("practice")
lattice::qqmath(fit_lmm_lme4, id = 0.05) ## id allows to see outliers
Plotting BLUPs with {lme4} r .emo("practice")
lattice::qqmath(lme4::ranef(fit_lmm_lme4))
Plotting residuals with {spaMM} r .emo("practice")
plot(fit_lmm_spaMM)
Note: it does not work on the slide, but should work in your RStudio...
Using simulated residuals r .emo("practice")
You can also rely on {DHARMa} as we have seen it for GLM:
library(DHARMa) plot(simulateResiduals(fit_lmm_spaMM, n = 1000))
Note: {DHARMa} works for both {spaMM} and {lme4}!
Random or fixed?
Choosing between fixed and random effects r .emo("info")
The choice between considering a predictor as a fixed or random effect can be difficult; there is a trade-off between both approaches.
Fixed
- no assumption about the distribution of the parameter values associated with each level of a predictor
- requires many datapoints for each additional level to get reliable results
- allows for the prediction of the effect of observed levels only (for factors)
- simple to handle
Random
- the values associated with the levels of a predictor follow a probability distribution (usually Gaussian)
- requires at least one datapoint for each additional level to get reliable results (more is better)
- requires many levels for the variance estimates to be reliable (5 being strict minimum)
- allow for the prediction of the effect of both observed and unobserved levels (for factors)
- more difficult to handle
Specifying multiple random effects
The lme4::Penicillin dataset r .emo("alien")
str(Penicillin) table(Penicillin$sample, Penicillin$plate)
The lme4::Penicillin dataset r .emo("info")
The random effects are crossed, i.e. we have now 2 independent random effects:
fit_peni <- fitme(diameter ~ 1 + (1|plate) + (1|sample), data = Penicillin)
fit_peni$lambda
The lme4::cake dataset r .emo("alien")
head(cake) str(cake)
The lme4::cake dataset r .emo("alien")
table(cake$recipe, cake$replicate, cake$temperature)
The lme4::cake dataset r .emo("info")
The random effect is nested within a fixed effect, so we must account for that:
fit_cake <- fitme(angle ~ recipe + temperature + (1|recipe:replicate), data = cake) fit_cake$lambda
Note:
- considering random effects as crossed here would be wrong. For example, it would imply that the random realization of the effect associated with replicate 1 from recipe A would be the same as that for the replicate 1 from recipe B.
The lme4::cake dataset r .emo("practice")
The random effect is nested within a fixed effect, so we must account for that, (alternative specification):
cake$replicate_tot <- factor(paste(cake$recipe, cake$replicate, sep = "_")) levels(cake$replicate_tot) fit_cake <- fitme(angle ~ recipe + temperature + (1|replicate_tot), data = cake) fit_cake$lambda
Note:
- whatever you do, the important thing is to make sure that replicates across recipes do not share the same names.
The carnivora dataset r .emo("alien")
carnivora$log_brain <- log(carnivora$Brain) carnivora$log_body <- log(carnivora$Weight) str(carnivora)
The carnivora dataset r .emo("info")
Genus names are unique across families:
length(unique(paste(carnivora$Genus, carnivora$Family))) == length(unique(carnivora$Genus))
coplot(log_brain ~ log_body | Family, data = carnivora)
The carnivora dataset r .emo("info")
Here we could consider one random effect nested into another one:
fit_carnivora <- fitme(log_brain ~ log_body + (1|Family/Genus), data = carnivora, method = "REML") fit_carnivora
The carnivora dataset r .emo("info")
Here we could consider one random effect nested into another one (alternative specification):
fit_carnivora_bis <- fitme(log_brain ~ log_body + (1|Family) + (1|Family:Genus), data = carnivora, method = "REML") fit_carnivora_bis
The carnivora dataset r .emo("info")
But since genus names are not (here) replicated across families, nested random effects are equivalent to crossed random effects:
fit_carnivora_ter <- fitme(log_brain ~ log_body + (1|Family) + (1|Genus), data = carnivora, method = "REML") fit_carnivora_ter
Checking the random structure r .emo("practice")
Checking the names of the BLUPs structure is an easy solution to check that you did specify the random structure correctly:
lapply(ranef(fit_carnivora), function(x) head(names(x))) ## or simply ranef(fit_carnivora)
Random slopes
Fitting a random slope model r .emo("practice")
(fit_rand_slope_spaMM <- fitme(log_brain ~ log_body + (log_body|Family) + (1|Genus), data = carnivora, method = "REML"))
The BLUPs for the slopes r .emo("info")
A random slope model is a model in which one slope is a random variable.
Here the effect of log_body is considered to vary across families:
ranef(fit_rand_slope_spaMM)$`( log_body | Family )`
Predictions r .emo("practice")
data_pred <- predict(fit_rand_slope_spaMM, newdata = expand.grid(log_body = range(carnivora$log_body), Family = levels(carnivora$Family), Genus = "new"), binding = "log_brain_pred") ## Tip: binding adds predictions to newdata ggplot() + geom_line(aes(y = log_brain_pred, x = log_body, colour = Family), data = data_pred) + geom_point(aes(y = log_brain, x = log_body, colour = Family), shape = 1, data = carnivora) + theme_minimal()
Testing the random slope using {lme4} r .emo("practice")
Asymptotic LRT are really poor for that, so better use parametric bootstrap:
fit_rand_slope_lme4 <- lmer(log_brain ~ log_body + (log_body|Family) + (1|Genus), data = carnivora) fit_no_rand_slope_lme4 <- lmer(log_brain ~ log_body + (1|Family) + (1|Genus), data = carnivora) pbkrtest::PBmodcomp(fit_rand_slope_lme4, fit_no_rand_slope_lme4)
Notes:
- fits differ by 2 degrees of freedom (1 for the fit of the variance of the random slope + 1 for its covariance with the random intercept).
pbkrtest::PBmodcomp()uses ML logLik.
Testing the random slope using {spaMM} r .emo("practice")
Asymptotic LRT are really poor for that, so better use parametric bootstrap r .emo("slow"):
fit_no_rand_slope_spaMM <- fitme(log_brain ~ log_body + (1|Family) + (1|Genus), data = carnivora, method = "REML") anova(fit_rand_slope_spaMM, fit_no_rand_slope_spaMM, boot.repl = 999)
Note:
- I used REML fit unlike
pbkrtest::PBmodcomp()(but it is not clear what the best choice is).
Generalized Linear Mixed-effects Models
GLM + LMM = GLMM r .emo("info")
$$\begin{array}{lcl} \mu &=& g^{-1}(\eta)\ \mu &=& g^{-1}(\mathbf{X}\beta + \mathbf{Z}b)\ \end{array} $$
with (as for GLM):
- $\text{E}(\text{Y}) = \mu = g^{-1}(\eta)$
- $\text{Var}(\text{Y}) = \phi\text{V}(\mu)$
Notes:
- If $g^{-1}$ is the identity function, $\phi = \sigma^2$ and $\text{V}(\mu) = 1$, we have the LMM.
- If $\mathbf{Z}b = 0$, we have the GLM.
- If $g^{-1}$ is the identity function, $\phi = \sigma^2$, $\text{V}(\mu) = 1$, and $\mathbf{Z}b = 0$, we have the LM.
If you combine everything we said about GLMs and LMMs, then you should know how to handle GLMMs.
Let illustrate this with an example!
The Flatwork dataset r .emo("alien")
Flatwork
The Flatwork dataset r .emo("alien")
str(Flatwork)
Fitting GLMM with {lme4} r .emo("practice")
(fit_fw_lme4_pois <- glmer(shopping ~ gender + (1|individual) + (1|month), family = poisson(), data = Flatwork))
Note:
- with {lme4} you cannot choose between ML and REML for GLMMs (it does a ML fit).
Fitting GLMM with {spaMM} r .emo("practice")
With {spaMM} you can still choose between ML and REML for GLMMs (REML sensu stricto does not exist for GLMMs, but {spaMM} has implemented a generalisation of the concept proposed in the statistical literature).
By default, as for LMMs, spaMM::fitme() produces a ML fit.
(fit_fw_spaMM_pois <- fitme(shopping ~ gender + (1|individual) + (1|month), family = poisson(), data = Flatwork))
Checking residuals r .emo("practice")
You can still test the assumptions of GLMMs with {DHARMa}:
fw_resid_spaMM_pois <- simulateResiduals(fit_fw_spaMM_pois) plot(fw_resid_spaMM_pois)
Overdispersion? r .emo("practice")
testDispersion(fw_resid_spaMM_pois)
Extra 0s? r .emo("practice")
testZeroInflation(fw_resid_spaMM_pois)
Alternative 1: negative binomial r .emo("practice")
(fit_fw_spaMM_nb <- fitme(shopping ~ gender + (1|individual) + (1|month), family = spaMM::negbin(), data = Flatwork))
Alternative 1: negative binomial r .emo("practice")
fw_resid_spaMM_nb <- simulateResiduals(fit_fw_spaMM_nb) plot(fw_resid_spaMM_nb)
Alternative 1: negative binomial r .emo("practice")
testDispersion(fw_resid_spaMM_nb)
Alternative 1: negative binomial r .emo("practice")
testZeroInflation(fw_resid_spaMM_nb)
Alternative 2: zero-augmented Poisson r .emo("practice")
library(glmmTMB) fit_fw_glmmTMB_zipois <- glmmTMB(shopping ~ gender + (1|individual) + (1|month), family = poisson(), data = Flatwork, ziformula = ~ 1)
Like {spaMM}, {glmmTMB} is extremely flexible package (it can fit models {spaMM} cannot fit, and the reverse is also true).
{glmmTMB} allows for both zero-inflation and hurdle models (using truncated distributions), as well as many covariance functions, modelling the residual variance...
There is no full-fledged parametric bootstrap methods implemented as far as I know, nor conditional AIC, but this could change in the near future.
The package is a wrapper around general numerical methods, so it may not be as robust as packages build with mixed models in mind (e.g. {spaMM} and {lme4}) in some cases; so convergence issues could happen more often.
Alternative 2: zero-inflated Poisson r .emo("practice")
{glmmTMB} is again compatible with {DHARMa}:
fw_resid_glmmTMB_zipois <- simulateResiduals(fit_fw_glmmTMB_zipois) plot(fw_resid_glmmTMB_zipois)
Alternative 2: zero-inflated Poisson r .emo("practice")
{glmmTMB} is also compatible with {DHARMa}:
testDispersion(fw_resid_glmmTMB_zipois)
Alternative 2: zero-inflated Poisson r .emo("practice")
{glmmTMB} is also compatible with {DHARMa}:
testZeroInflation(fw_resid_glmmTMB_zipois)
Alternative 3: zero-inflated negative binomial r .emo("practice")
(fit_fw_glmmTMB_zinb <- glmmTMB(shopping ~ gender + (1|individual) + (1|month), family = nbinom1(), data = Flatwork, ziformula = ~ 1))
Alternative 3: zero-inflated negative binomial r .emo("practice")
fw_resid_glmmTMB_zinb <- simulateResiduals(fit_fw_glmmTMB_zinb) plot(fw_resid_glmmTMB_zinb)
Alternative 3: zero-inflated negative binomial r .emo("practice")
testDispersion(fw_resid_glmmTMB_zinb)
Alternative 3: zero-inflated negative binomial r .emo("practice")
testZeroInflation(fw_resid_glmmTMB_zinb)
What is the best alternative for our GLMM? r .emo("practice")
AIC(fit_fw_spaMM_pois)
print(AIC(fit_fw_spaMM_pois))
AIC(fit_fw_spaMM_nb)
print(AIC(fit_fw_spaMM_nb))
AIC(fit_fw_glmmTMB_zipois, fit_fw_glmmTMB_zinb) #Marginal AIC
Note:
- the zero-inflated negative binomial model fits the data much better than the alternative models we tried.
Predictions for GLMM? r .emo("info")
-
building conditional predictions is similar to what we have seen in GLM and LMM: just use
predict()(orfitme::pdep_effects()) for that. -
building marginal predictions is more complex because you need to integrate the predictions over the distribution of the random effects... There is no function doing that automatically and this is a little too advanced for this course.
What you need to remember r .emo("goal")
- what mixed models (LMM & GLMM) are
- how to use {lme4}, {spaMM} and {glmmTMB} to fit mixed models
- how to make predictions, tests, CI...
- how to check assumptions
- when to consider effects as fixed or random
- how to implement different random effects structure
Table of contents
Mixed-effects models
- 4.0 Introduction to LMM & GLMM
- 4.1 Solving LM problems using LMM
- 4.2 A showcase of some useful applications
[Back to main menu](./Title.html#2)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.