constructing the components of a mixed model

back to the general definition

\[ \begin{split} \underbrace{Y_i}_{\text{response}} & \sim \overbrace{\text{Distr}}^{\substack{\text{conditional} \\ \text{distribution}}}(\underbrace{g^{-1}(\eta_i)}_{\substack{\text{inverse} \\ \text{link} \\ \text{function}}},\underbrace{\phi}_{\substack{\text{scale} \\ \text{parameter}}}) \\ \underbrace{\boldsymbol \eta}_{\substack{\text{linear} \\ \text{predictor}}} & = \underbrace{\boldsymbol X\boldsymbol \beta}_{\substack{\text{fixed} \\ \text{effects}}} + \underbrace{\boldsymbol Z\boldsymbol b}_{\substack{\text{random} \\ \text{effects}}} \\ \underbrace{\boldsymbol b}_{\substack{\text{conditional} \\ \text{modes}}} & \sim \text{MVN}(\boldsymbol 0, \underbrace{\Sigma(\boldsymbol \theta)}_{\substack{\text{variance-} \\ \text{covariance} \\ \text{matrix}}}) \end{split} \]

fixed-effect model matrix

• \(\boldsymbol X\) (\(n \times p_\textrm{fix}\))
• terms, input variables, predictor variables
• Wilkinson-Rogers notation¹
• defined by model.matrix(), terms()²
• details hard to find!
• default intercept term (1); -1 or +0 to suppress
• numeric input variables → single predictor variable
• variable transformations on the fly
• contrasts automatically constructed for factors (unless +0)
• interaction : = pairwise products of predictor variables
• data-dependent bases: poly(), splines::ns(), etc.

packages

library(Matrix)
library(reformulas)

example

terms(~ mpg * disp, data = mtcars)
model.matrix(hp~ mpg * disp, data = mtcars)

random-effect formulas

• formula: typically ( f | g ) (f = varying term, g = cluster/grouping variable)
• g can be expanded into separate terms ((f | g1/g2) → (f | g1 + g1:g2)); *, + may be expanded (reformulas package)
• f = any set of terms that vary within levels of g
• dialect differences
  • random component may be specified separately from fixed (e.g. nlme::lme, MCMCglmm)
  • … or as part of the same formula (e.g. lme4, glmmTMB)
  • covariance structures may be specified (brms, MCMCglmm, glmmTMB) as (e.g.) cs( f | g )
  • AS-REML-like syntax (e.g. MCMCglmm): e.g. ~us(f):g
  • extensions for multi-membership models, shared variables between terms, …
• each RE term (after separation) is independent

from RE formula to model objects: RE model matrix \(\boldsymbol Z\)

• details in lmer vignette³
• \(\boldsymbol Z\)
  • Khatri-Rao product of \(\boldsymbol J\) (indicator matrix for cluster variable g) with model matrix for f
  • e.g.:

f <- factor(rep(1:3, each = 2))
Ji <- t(as(f, "sparseMatrix"))
Xi <- cbind(1, rep.int(c(-1, 1), 3L))
Zi <- t(Matrix::KhatriRao(t(Ji), t(Xi)))

• make \(\boldsymbol J\) a sparse matrix, \(\boldsymbol Z\) inherits structure
• concatenate \(\boldsymbol Z\) from multiple terms columnwise

from RE formula to model objects: RE covariance matrix \(\Sigma\)

• \(\Sigma_i\) (for a single term) is block-diagonal (one block per cluster); terms are blocks too
• we may be interested in the lower-triangular Cholesky factor (\(\Sigma = \Lambda \Lambda^\top\))
  • or the precision matrix (for a Gaussian model, “\(x_{i}\) is conditionally independent of \(x_{j}\) \(\to\) \(\sigma^{-1}_{ij} = 0\); enables sparse computations

lme4 example

dd <- expand.grid(f1 = factor(1:3), f2 = factor(1:2), g1 = factor(1:5), g2 = factor(1:6))
dd$y <- 1
form <- y ~ 1 + (f1 | g1) + (f2 | g2)
rt <- reformulas::mkReTrms(
                      bar = findbars(form),
                      fr = model.frame(subbars(form), data = dd))
Sigma <- rt$Lambdat %*% t(rt$Lambdat)
image(Sigma)

references

1.

Wilkinson, G. N. & Rogers, C. E. Symbolic description of factorial models for analysis of variance. Applied Statistics 22, 392–399 (1973).

2.

Chambers, J. M. & Hastie, T. J. Statistical Models. in Statistical Models in S (Routledge, 1992).

3.

Bates, D., Mächler, M., Bolker, B. M. & Walker, S. C. Fitting linear mixed-effects models using lme4. Journal of Statistical Software 67, 1–48 (2015).