Mind the Gap: Mixed Models in R

Problem: dispersed functionality

Table of a some example functionality and the packages where it is enabled. This is not a comprehensive list of possible functionality in mixed models or what it implemented in R.

The most popular mixed model package is likely lme4, perhaps followed by nlme (it’s difficult to quantify this because downloads don’t always reflect usage, and since nlme is now part of base R, download stats for that are really wonky).

These packages by themselves can good aims, but the problems arise in complex situations that require combining functionality. While lme4 is very powerful (and fast!), it lacks essential functionality for including spatial covariates. nlme does have that functionality (although not for Matérn, arguable one of the most popular spatial covariance structures), but it is difficult (but not impossible) to specify crossed random effect (they are nested by default in lme() calls). There are several packages that extend lme4 but to permit inclusion of a kinship based covariance matrix, but combining that with spatial covariates is not possible. And of course, lme4 famously does not compute p-values, a limitation that can be overcome by computing them manually (fun!!!) or using lmerTest. But, often users can only use one lme4 package extender at a time, leaving them to choose which functionality they need most. And while I agree with the spirit of why p-values are not included in the default ANOVA output for lme4, this is rarely a persuasive argument for my clients, who want p-values come hell or high water.

Problem: inconsistent syntax

Inconsistent syntax is a long standing problem with R and mixed models or any modelling package in R has this problem (and for all those who say “use Tidymodels!” - agreed, and this is addressed later in the post).

Imagine we want to construct a linear mixed model for one fixed effect, x, and one random effect, rep, for a dependent variable, y. These variables are in a data.frame called “mydata” where y and x are numeric, while rep is an unordered factor.

The model:

\[y_{ij} = x_i + rep_j + \epsilon_{ij}\]

where the error terms, \(\epsilon_{ij}\) are distributed as thus:

\[\epsilon_{ij} \sim N(0, \sigma_{error})\] and \(rep_j\) is distributed as:

\[rep_j \sim N(0, \sigma_{rep})\]

Below are syntax examples for this model in 4 major mixed modelling packages in R.

lme4::lmer(y ~ x + (1|rep), data = mydata,  
        na.action = na.exclude, REML = TRUE)

nlme::lme(y ~ x, random = ~ 1|rep, data = mydata, 
       na.action = na.exclude, method = “REML”)

asreml::asreml(y ~ x, random = ~ rep, data = mydata, 
              na.method.X = 'exclude', 
             na.method.Y = ‘exclude’)

rrBLUP::mixed.solve(y = mydata$y, x = {fixed_design_matrix}, 
                    z = {random_design_matrix}, SE = TRUE)

Problem: varied output

There are so many special rules to follow.

1. Use lmerTest to output p-values for lme4 (but this doesn’t work with lme4GS or other add-on packages).
   
2. emmeans can extract fixed effects, but not BLUPs! (But, marginaleffects can calculate BLUPs! 🤩🤩🤩) Instead, check the package for package functions (fingers crossed that BLUP calculation is enabled!) and deal with the wildly varying syntax.
   
3. Not all these packages have methods written to extract output, so you have to hunt through the model object or the summary of the model object. This is annoying (and for many, it can be very mystifying) and results in very brittle code.
   
4. The AR1xAR1 covariance structure (kronecker product of AR1 in two dimensions, each with its own \(\rho\)) is implemented in sommer, but it does not estimate \(\rho\); breedR can do that (but it’s not open source); INLA can estimate \(\rho\), but the package uses a Bayesian approach and takes awhile to install and learn.
   
5. …and so much more.

Unsurprisingly, people have a hard time navigating this!

So, users end up implementing sub-optimal solutions (e.g. running aov(), equivalent to PROC ANOVA in SAS, when a mixed model is needed). I can (usually) help them, but is this the right long-term solution? Once these clients leave my office, are they ready to implement different mixed models on their own? And what about the people who never consult with a statistician? Mixed models are a cornerstone of agricultural research (and likely academic research overall), yet the tools available to do this correctly are insufficient.

multilevelmod package

This is part of the Tidymodels family, whose goal (among many) is to standardize the syntax for mixed modelling in the R-verse. I am so glad this package exists; it has been sorely needed for years and does simplify the syntax. However, it addresses 4 packages and an extremely limited number of scenarios:

Again, these package authors are busy, and this task of unifying the mixed model ecosystem is difficult, so no shade intended towards these folks! My point is that implementing a mixed model is still rather challenging from from the user standpoint.

broom.mixed package

This is a fantastic package that supports a large number of packages (see below). The highlighted packages are those also supported by multilevelmod.

• brms
  
• gamlss
  
• GLMMadaptive
  
• glmmADMB
  
• lme4
  
• mcmc
  
• MCMCglmm
  
• mediation
  
• nlme
  
• pscl
  
• rstanarm
  
• TMB

multilevelmod and broom.mixed work like their counterparts in the Tidymodels world. I don’t want to go into the syntax, but you can find examples here and here.

Many accessory packages:

There are loads of accessory packages written to ease the burden of extracting output from a mixed model object and/or extending the functionality. These add to the general complexity of specifying a mixed model in R and hence are both a solution and part of the problem.

It’s a big task to name all helper packages for mixed models that I have no interest in doing comprehensively, so here is a short list:

• afex: conduct ANOVA
• emmeans and marginaleffects: 2 helpful packages for extracted BLUE’s and in in the case of marginaleffects, BLUPs.
• lmerTest: get those p-values from lme4!

There are more (see the google sheets links listed above), but wow, keeping track of this is a major undertaking.

Tidymodels style guide

• Naming conventions for arguments
  
• How information should be returned in ‘predict’ calls, ‘summary’ statements
  
• what the user-facing function should look like
  
• usage of S3 method conventions
  
• how parallel processing occurs
  
• Tidyverse style guide (opinionated)

Missing from this guide are conventions on the standardized output for (mixed) linear models (log likelihood, residuals, variance components, etc). These things are “missing” because the Tidymodels guide is broadly written for a diverse group of packages/modelling applications in which the standard output from a linear model is not possible or desired. However, in the mixed model landscape, there is common output that users expect; hence it would be useful to codify standards for these types of models.