Linear mixed-effect models (LMMs)

Recall, blocking helps control variability by making treatment groups more alike. Experimental units are divided into subsets (called blocks) so that units within the same block are more similar than units from different subsets or blocks. Blocking is a technique for dealing with nuisance factors. A nuisance factor is a factor that has some effect on the response, but is of no interest (e.g., age class).

Fixed effects are terms (parameters) in a statistical model which are fixed, or non-random, quantities (e.g., treatment group’s mean response). For the same treatment, we expect this quantity to be the same from experiment to experiment.

Random effects are terms (parameters) in a statistical model which are considered as random quantities or variables (e.g., block id). Specifically, terms whose levels are a representative sample from a population, and where the variance of the population is of interest should be allocated as random. Setting a block as a random effect allows us to infer variation between blocks as well as the variation between experimental units within blocks.

Key idea: Partition known sources of variation which are unimportant to key scientific question(s) to improve precision of comparisons between treatment means.

A Randomised Controlled Block Design (RCBD)

The rcbd data:

To use predictmeans later on we have to ensure that the relevant variables are coded as factors:

Run as a random effect

Note Both the lmerTest and lme4 packages export a function called lmer, which fits linear mixed effect models. For all estimation purposes etc. they essentially do the same think. However, the structure of the output differs and this MATTERS when we later come to use other funtions on the outputs. Therefore, I take care below to specify which version of the lmer function I use by explicitly calling it from either lmerTest or lme4 using the :: syntax.

There are, confusingly, two ways of fitting the same model. For inference we require both!

Option 1 uses the lmer function from the lme4 package:

Option 2 uses the lmer function from the lmerTest package:

As you can see they give the same output! Why bother, you might ask?! This will become apparent later on.

Inference about the random effects

We have two variance components

  • Between Groups (Runs) \(\hat{\sigma^2}_{\text{Run}}\) = 1.479
    • Within Runs (between observations) \(\hat{\sigma_2}\) = 1.447

Note that aov() presents the same information, but in a different way:

  • Within Runs (Residuals) \(\hat{\sigma}_2\) = 1.447 (same as lmer)
  • Between Run variance = \(\hat{\sigma}^2\) + \(3\:\hat{\sigma}^2_{\text{Run}}\) so \(\hat{\sigma}^2_{\text{Run}} = \frac{5.883 - \hat{\sigma}^2 }{3} = \frac{5.883 - 1.447}{3} = 1.479\)

Inference about the fixed effects

Specifying Run as random effect changes our estimated baseline (i.e., Intercept coefficient) as now and effect due to Run is attributed to the structural component of the model.

We can interpret the fixed effects of a LMM as we might for a linear model (now the Intercept estimate changes depending on Run).

What about an ANOVA table? NOTE that specifying the type of \(SS\) (e.g., Type I, II, or III) in an anova call only works with a model fitted using lmerTest::lmer:

Now, as we only have a single fixed effect in our model (Surgery) the ANOVA Type I and Type II tables above are equivalent!

Inference using predictmeans() (Note: output will be the same for the lmerTest model.)

Using the elements of the predictmeans object (as in the last section) we can extract the pairwise comparison information:

Gives us 1) pairwise differences (upper diagonal) and the LSD values (lower diagonal), and 2) the pairwise comparison statistic (the \(t\)-statistic in this case) on the upper diagonal and the associated p-value’s on the lower diagonal.

So is there a pairwise difference in means? Let’s organise the information in a table. Below we read a specifically designed function to do this from the given URL.

Have a look at the CIs and p-values. What do you conclude?

We could also plot the pairwise comparisons using emmeans

A Split-plot design

To use predictmeans later on we have to ensure that the relevant variables are coded as factors:

Animal as a random effect

Using aov()

Using lmer() (from lmeTest and lmer4)

Recall that specifying the type of \(SS\) (e.g., Type I, II, or III) in an anova call only works with a model fitted using lmerTest::lmer. As we now have specified an interaction model then the type of \(SS\) calculated will have an effect on inference for unbalanced designs.

Note here, though, our design is balanced hence the ordering of terms in our model does not make a difference (no need to specify Type II \(SS\))

Pairwise comparison of time means

When a design has blocking, to get summary stats using predictmeans you should fit the model using lme4::lmer():

Calling predictmeans alone produces plots:

## 
## $ciPlot

Recall from previous sections that the LSD vale is essentially the buffer around the point estimate (the radius of the CI if you like), beyond the limit of which we might believe there to be a significant difference (a bit lax with terminology there!).

How would be interpret the above plot? I would conclude that there is a big difference between Diabetic and Healthy for InnerLV wall, but no difference for the outerLV wall. We can construct the pairwise comparison table:

What do the CI coverages indicate?

We can also use emmeans to plot the pairwise comparisons

From all the output above, what do you conclude about the interaction effect, if any!

A repeated measures design

Below these data are plotted.

Animal as a random effect

Using aov()

Using lmer() (from lmerTest and lme4)

Recall that specifying the type of \(SS\) (e.g., Type I, II, or III) in an anova call only works with a model fitted using lmerTest::lmer. As we now have specified an interaction model then the type of \(SS\) calculated will have an effect on inference for unbalanced designs.

Pairwise comparison of time means

As above, when a design has blocking, to get summary stats using predictmeans you should fit the model using lme4::lmer():

TASK What inference do you draw? You may also find the outputs from the code below useful.