• Understand types of random effects structures

• Understand how random effects are estimated

• Understand restricted maximum likelihood

• Understand approaches to make inference from mixed models

This online tool is quite helpful

We have seen how to write a general linear model as

\[ \mathbf{y} = \mathbf{X} \boldsymbol{\beta} + \mathbf{e} \\ \mathbf{e} \sim \text{MVN}(\mathbf{0}, \sigma^2 \mathbf{I}) \]

where

\(\mathbf{X}\) is the design matrix

\(\boldsymbol{\beta}\) contains the fixed effects of \(\mathbf{X}\) on \(\mathbf{y}\)

\(\mathbf{e}\) are the model errors and \(\mathbf{I}\) is the identity matrix

We can extend the general linear model to include both of fixed and random effects (a mixed effects model)

\[ \mathbf{y} = \mathbf{X} \boldsymbol{\beta} + \mathbf{Z} \boldsymbol{\alpha} + \mathbf{e} \\ \boldsymbol{\alpha} \sim \text{MVN}(\mathbf{0}, \sigma^2 \mathbf{D}) \\ \mathbf{e} \sim \text{MVN}(\mathbf{0}, \sigma^2 \mathbf{I}) \]

where

\(\mathbf{Z}\) is a design matrix containing a mix of factors and covariates

\(\boldsymbol{\alpha}\) are the random effects and \(\mathbf{D}\) is a correlation matrix

\(\mathbf{e}\) are the model errors and \(\mathbf{I}\) is the identity matrix

\[ \text{Var}(DATA) = \text{Var}(MODEL) + \text{Var}(ERRORS) \]

\[ \text{Var}(DATA) = \text{Var}(MODEL) + \text{Var}(RE's) + \text{Var}(ERRORS) \]

\[ \mathbf{y} = \mathbf{X} \boldsymbol{\beta} + \mathbf{Z} \boldsymbol{\alpha} + \mathbf{e} \\ \Downarrow \\ \text{Var}(\mathbf{y}) = \text{Var}(\mathbf{\mathbf{X} \boldsymbol{\beta}}) + \text{Var}(\mathbf{Z} \boldsymbol{\alpha}) + \text{Var}(\mathbf{e}) \]

\[ \text{Var}(\mathbf{y} | \mathbf{\mathbf{X} \boldsymbol{\beta}}) = \text{Var}(\mathbf{Z} \boldsymbol{\alpha}) + \text{Var}(\mathbf{e}) \\ \Downarrow \\ \begin{aligned} \mathbf{V} &= \mathbf{Z} \text{Var} (\boldsymbol{\alpha}) \mathbf{Z}^\top + \text{Var}(\mathbf{e}) \\ &= \mathbf{Z} (\sigma^2 \mathbf{D}) \mathbf{Z}^\top + \sigma^2 \mathbf{I} \\ &= \sigma^2 (\mathbf{Z} \mathbf{D} \mathbf{Z}^\top + \mathbf{I}) \end{aligned} \]

Recall that we think of likelihoods in terms of the observed data

Recall that we think of likelihoods in terms of the observed data

But the random effects in our model are unobserved random variables, so we need to integrate them out of the likelihood

The log-likelihood for the fixed effects \(\boldsymbol{\beta}\)

\[ \log \mathcal{L}(\mathbf{y}; \boldsymbol{\beta}, \sigma^2) = - \frac{1}{2} \log \left| \mathbf{V} \right| - \frac{1}{2}( \mathbf{y} - \mathbf{X} \boldsymbol{\beta})^\top \mathbf{V}^{-1} (\mathbf{y} - \mathbf{X} \boldsymbol{\beta}) \]

This leads us to our familiar statement for the weighted least squares estimate for \(\boldsymbol{\beta}\)

\[ \begin{aligned} \hat{\boldsymbol{\beta}} &= \min ~ (\mathbf{y} - \mathbf{X} \boldsymbol{\beta})^{\top} \mathbf{V}^{-1} (\mathbf{y} - \mathbf{X} \boldsymbol{\beta}) \\ &= (\mathbf{X}^{\top} \mathbf{V}^{-1} \mathbf{X}) \mathbf{X}^{\top} \mathbf{V}^{-1} \mathbf{y} \end{aligned} \]

Our variance estimate for \(\boldsymbol{\beta}\) is then

\[ \text{Var}(\hat{\boldsymbol{\beta}}) = (\mathbf{X}^{\top} \mathbf{V}^{-1} \mathbf{X})^{-1} \]

The log-likelihood for the random effects is given by

\[ \begin{aligned} \log \mathcal{L}(\mathbf{y}; \boldsymbol{\beta}, \sigma^2) = - \frac{\sigma^2}{2} &- \frac{1}{2 \sigma^2}( \mathbf{y} - \mathbf{X} \boldsymbol{\beta} - \mathbf{Z} \boldsymbol{\alpha})^\top (\mathbf{y} - \mathbf{X} \boldsymbol{\beta} - \mathbf{Z} \boldsymbol{\alpha}) \\ &- \frac{1}{2} \left| \mathbf{Z} \mathbf{D} \mathbf{Z}^\top\right| - \frac{1}{2} \boldsymbol{\alpha}^\top (\mathbf{Z} \mathbf{D} \mathbf{Z}^\top)^{-1} \boldsymbol{\alpha} \end{aligned} \]

This leads to the best linear unbiased predictor for \(\boldsymbol{\alpha}\)

\[ \hat{\boldsymbol{\alpha}} = \sigma^2 (\mathbf{Z} \mathbf{D} \mathbf{Z}^\top) \mathbf{Z}^\top \mathbf{V}^{-1} (\mathbf{y} - \mathbf{X} \boldsymbol{\beta}) \]

Recall the difference between the sample variance and MLE estimate

\[ \hat{\sigma}^2 = \frac{1}{n-1} \sum_{i = 1}^n (y_i - \bar{y}) \\ ~ \\ \hat{\sigma}^2_{MLE} = \frac{1}{n} \sum_{i = 1}^n (y_i - \bar{y}) \]

such that

\[ \hat{\sigma}^2_{MLE} = \frac{n - 1}{n} \hat{\sigma}^2 \]

Maximum likelihood estimates of the variances in mixed models are also biased

To correct for this, we use restricted maximum likelihood

REML works by

1. estimating the fixed effects \((\hat{\boldsymbol{\beta}})\) via ML

2. using the \(\hat{\boldsymbol{\beta}}\) to estimate the \(\hat{\boldsymbol{\alpha}}\)

lme4 makes this easy for us

library(lme4)

lmer()

Benefits:

• better performance on large problems

• better for crossed random effects

• profile likelihood confidence intervals on RE’s

library(nlme)

lme()

Benefits:

• Fitting autocorrelation structure (spatial and temporal)

• Fitting heteroscedasticity

With random effects models, we can’t use our standard inference tools because we don’t know the distributions for our test statistic

({lme4} won’t give \(p\)-values)

We can use a likelihood ratio test for nested models, but the assumption of a \(\chi^2\) distribution can be poor

We can also use \(F\)-tests to evaluate a single fixed effect, but again the assumption of a \(F\) distribution can be poor

We can use bootstrapping to conduct likelihood ratio tests

1. simulate data from the simple model

2. fit simple & full model and calculate likelihood ratio

3. see where test statistic falls within estimated distribution from (2)

We can report parameter estimates and CI’s via bootstrapping

We can generate predictions given fixed and random effects and estimate their uncertainty via bootstrapping

Recall that \(AIC = 2 k - 2 \log \mathcal{L}\)

The problem with mixed effects models is that it’s not clear what \(k\) equals

It works well to select among fixed effects if random effects are held constant

To use AIC, we can follow these steps

1. Fit a model with all of the possible fixed-effects included

2. Keep the fixed effects constant and search for random effects

3. Keep random effects as is and fit different fixed effects

General advice for random effects

• Assume we know what random effects to include and move forward to evaluating fixed effects

• Work out what RE’s might be needed from understanding of the sampling situation

Other options include

• BIC

• cross-validation

• Once random effects are chosen, select fixed effects

• Be cautious with AIC (use ML to get AIC estimates)

• AIC can’t help us choose our random effects

• Inference will generally require bootstrapping