What a GLMM can (and can’t!) do

The utility of GLMMs covers a wide swath of problems across disciplines, and includes two very common types of response distributions: (1) count data, and (2) binary data. Alongside these two response variables, we can handle repeated measures, nested designs, and longitudinal data with this approach.¹

It is, however, important to point out up front what a GLMM cannot do. The “generalized” in GLMM refers specifically to the ability to accommodate response distributions in the exponential family with a link function² — not to flexibility in any broader sense. The model is still fundamentally linear on the link-transformed scale, meaning predictors enter additively with constant coefficients, and the framework has no inherent ability to capture truly nonlinear mean structures. Common structures of this kind that are better modeled with other tools are asymptotic growth curves, dose-response relationships with inflection points, or other phenomena where the functional form is inherently nonlinear. For those situations, nonlinear mixed effects models are more appropriate. The model is also not inherently flexible in a machine learning type of way.

Unhelpfully, the reader might feel as though the scope of the GLMM, from this short introduction, is both incredibly wide and also quite narrow. This is in many ways true, and requires some delicacy in both (a) our selection or rejection of a GLMM as a suitable tool for our problem, and (b) the inference using a model like this. GLMMs are actually relatively new to the statistical universe, their formalism was only brought together by N. E. Breslow and Clayton (1993), building on the GLM framework set out (in an excellent text worth a read) by McCullagh and Nelder (1989). Given that they are both new on the scene and so flexible, their implementation is an active area of statistical research.

In the rest of this guide we will try to summarize some of the important theoretical underpinnings of a GLMM, and then discuss some important considerations in the implementation of these models. While we will try to compare and contrast various estimation procedures, as noted above, accurately doing estimation for GLMMs is an active area of research, so we err on the side of more established methods and implementation tools, though the adventurous modeler may check out the FAQs for pointers to some newer approaches.

flowchart LR
    Q1{"Non-normal response?<br/>binary, count, proportion"}
    Q2{"Observations grouped?<br/>grouped or repeated"}
    Q3{"Linear on link scale?<br/>no inherent curve"}
    E1["Linear models<br/>OLS or LMM"]
    E2["GLM<br/>no clustering"]
    E3["NLMM<br/>nonlinear mean"]
    GLMM["GLMM<br/>continue with this guide"]

    Q1 -->|Yes| Q2
    Q1 -->|No| E1
    Q2 -->|Yes| Q3
    Q2 -->|No| E2
    Q3 -->|Yes| GLMM
    Q3 -->|No| E3

    classDef question fill:#F1EFE8,stroke:#888780,color:#2C2C2A
    classDef exit    fill:#EEEDFE,stroke:#7F77DD,color:#3C3489
    classDef outcome fill:#E1F5EE,stroke:#1D9E75,color:#085041

    class Q1,Q2,Q3 question
    class E1,E2,E3 exit
    class GLMM outcome
    linkStyle default stroke:#2C2C2A,stroke-width:1.5px,fill:none

Theoretical Intuition

Ultimately, we are best served, when thinking about a new model or approach, by starting with a motivating example. Suppose in some midwestern American state, each farmer plants some number of plants on their farm. They meticulously keep track of the survival of a given plant through the course of the growing season, and at the end of the growing season they have survival data for each plant. These are obviously very small farms or very industrious farmers…

Since these data are considered to be Bernoulli (i.e. composed of 0s and 1s), it is likely that our primary assumption of standard linear regression (that of normality) is immediately violated. If you need a refresher on what that assumption actually entails, and on the other core assumptions of the linear model, we recommend reviewing the OLS regression guide before continuing. In brief, the four assumptions of the linear model encoded in \(\mathbf{Y}|\mathbf{X} \sim \mathcal{N}(\mathbf{X}\boldsymbol{\beta}, \sigma^2\mathbf{I})\) are: normality of errors, linearity of the mean in the predictors, homoscedasticity, and independence of observations. Bernoulli data violate the first three immediately because the support is \(\{0,1\}\) rather than \(\mathbb{R}\), and the variance \(p(1-p)\) depends on the mean rather than being constant.

Formalizing our Model

Alright, so we know at this point, we need two things to deal with our farm data. We need to account for the bernoulli response distribution, which we’ll do with a link function, and we need to account for our error non-independence with a random effect. We now have the tools to write down our model. Let’s say that our survival data, \(\mathbf{Y}\) still has the linear predictor \(\eta = \mathbf{X} \boldsymbol{\beta}+ \mathbf{Z}\mathbf{b}\) which we just derived in the previous section, and that our expected value \(E(\mathbf{Y}|\mathbf{b}) = \mu = g^{-1}(\eta)\).

In general, we can say that the probability model for a GLMM is

\[\begin{equation} \mathbf{Y}|\mathbf{b} \sim \mathcal{D}(g^{-1}(\mathbf{X}\boldsymbol{\beta}+ \mathbf{Z}\mathbf{b}), \phi) \end{equation}\]

where \(\mathcal{D}\) is the distribution family (e.g. Bernoulli), which is characterized by the mean (where the mean, as we’ve seen, is \(\mu = g^{-1}(\mathbf{X}\boldsymbol{\beta}+ \mathbf{Z}\mathbf{b})\)), and a scale parameter, \(\phi\).

In our case, we don’t need a scale parameter since the Bernoulli distribution takes only one parameter. So, for our model,

\[\begin{equation} \mathbf{Y}|\mathbf{b} \sim \text{Bernoulli}(g^{-1}(\mathbf{X}\boldsymbol{\beta}+ \mathbf{Z}\mathbf{b})) \end{equation}\]

Link Functions

We have already seen, in the motivating example above, why we cannot simply apply a standard linear model to bernoulli data: the linear predictor \(\eta = \mathbf{X}\boldsymbol{\beta}\) is unbounded, ranging over \((-\infty, +\infty)\), while a probability (since our survival data is technically a probability) is constrained to \([0, 1]\). A link function is the tool that resolves this incompatibility. Formally, a link function \(g\) is a monotone⁵, differentiable function that maps the conditional mean \(\mu = \text{E}(Y \mid \mathbf{x})\) to the linear predictor \(\eta\), such that \(g(\mu) = \eta = \mathbf{x}^\top\boldsymbol{\beta}\). Equivalently, and as we saw in the probability model above, we can express the conditional mean as \(\mu = g^{-1}(\eta)\), where \(g^{-1}\) is the inverse link function. Importantly, we are not transforming the data — we are transforming the mean of the response distribution, allowing the linear predictor to operate on an unconstrained scale while the response itself remains on its natural scale.

The choice of link function is not arbitrary. Each distribution in the exponential family has a canonical link function, which arises naturally from the mathematical form of the distribution itself and has convenient theoretical properties for estimation. For the bernoulli distribution, the canonical link is the logit function, \(g(\mu) = \log\left(\frac{\mu}{1-\mu}\right)\), which maps probabilities in \([0,1]\) to the real line \((-\infty, +\infty)\). Its inverse, \(g^{-1}(\eta) = \frac{1}{1 + e^{-\eta}}\), is the familiar sigmoid function, which squashes the linear predictor back into a valid probability. For count data following a Poisson distribution, the canonical link is the log function, \(g(\mu) = \log(\mu)\), which maps positive means to the real line, with inverse \(g^{-1}(\eta) = e^{\eta}\).

For the normal distribution, the canonical link is simply the identity, \(g(\mu) = \mu\), which recovers the standard linear model as a special case — confirming that everything we discussed in the conditional normality section above is just a GLM with an identity link.

It is worth being precise about what the link function is not doing. A common misconception is that the link function is a data transformation — that we are log-transforming our response variable and then fitting a linear model, which is not the case. When we fit a GLM with a log link, we are saying that \(\log(\mu) = \mathbf{x}^\top\boldsymbol{\beta}\), not that \(\log(y_i) \sim \mathcal{N}(\mathbf{x}_i^\top\boldsymbol{\beta}, \sigma^2)\). The response \(y_i\) retains its original distribution — Poisson, bernoulli, or otherwise — and the likelihood is evaluated accordingly. Transforming the data and then fitting a linear model changes the likelihood entirely, which has real consequences for inference and is not equivalent to fitting a GLM. The link function operates on the mean structure; the distributional assumption governs how the data scatter around that mean.

Returning to our farm example, we will use the logit link, giving us the model \(\log\left(\frac{\mu_i}{1-\mu_i}\right) = \mathbf{x}_i^\top\boldsymbol{\beta}\), or equivalently \(\mu_i = \frac{1}{1 + e^{-\mathbf{x}_i^\top\boldsymbol{\beta}}}\). The coefficients \(\boldsymbol{\beta}\) are then interpreted on the log-odds scale: a unit increase in a predictor shifts the log-odds of survival by \(\beta_j\), or equivalently multiplies the odds of survival by \(e^{\beta_j}\). This is the price of using a non-identity link — the linear predictor is no longer on the same scale as the response, and coefficients require careful interpretation. We will return to this point when we discuss the output of our fitted model below.

We now have formulated our full model! Unfortunately, this model is not quite as simple to fit as one that might be amenable to OLS (ordinary least squares), or even a direct usage of maximum likelihood. The source of this difficulty is actually the likelihood function itself,

The Likelihood

Most of the challenges of fitting a GLMM come from the fact that the response is non-normal, which means we cannot use standard maximum likelihood directly. To understand why, it helps to first write down what the likelihood actually is for our model.

Recall that in our farm example, we have \(i = 1, \ldots, I\) farms, each with \(j = 1, \ldots, n_i\) plants. Each plant has a binary survival outcome \(y_{ij}\), a vector of predictors \(\mathbf{x}_{ij}\), and shares a farm-level random effect \(\mathbf{b}_i\). If we knew the value of the random effects \(\mathbf{b}\), the likelihood would be straightforward — conditional on \(\mathbf{b}_i\), the observations within each farm are independent Bernoulli trials, and we could write the conditional likelihood as:

\[\begin{equation} L(\boldsymbol{\beta} \mid \mathbf{b}) = \prod_i \prod_j f(y_{ij} \mid \mathbf{b}_i) \end{equation}\]

The left-hand side, \(L(\boldsymbol{\beta} \mid \mathbf{b})\), is the probability of observing all of our data simultaneously, given the fixed-effect parameters \(\boldsymbol{\beta}\) and the random effects \(\mathbf{b}\). The double product \(\prod_i \prod_j\) is just a compact way of saying “multiply together one term for every plant on every farm”. The term \(f(y_{ij} \mid \mathbf{b}_i)\) is the Bernoulli probability mass function evaluated at plant \(j\)’s observed outcome:

\[\begin{equation} f(y_{ij} \mid \mathbf{b}_i) = \mu_{ij}^{y_{ij}}(1 - \mu_{ij})^{1 - y_{ij}} \end{equation}\]

with \(\mu_{ij} = g^{-1}(\mathbf{x}_{ij}^\top\boldsymbol{\beta} + \mathbf{Z}_i\mathbf{b}_i)\). This is just saying: for each plant, the probability of the observed outcome is \(\mu_{ij}\) if the plant survived and \(1 - \mu_{ij}\) if it did not. The product over all plants then gives the joint probability of the entire dataset, assuming observations are independent once we condition on the random effects.

But we do not know \(\mathbf{b}\) — the random effects are unobservable. We only know that they come from the distribution \(\mathbf{b}_i \sim \mathcal{N}(\mathbf{0}, G)\). To obtain the marginal likelihood (the likelihood in terms of only the parameters we actually want to estimate, \(\boldsymbol{\beta}\) and \(G\)) we have to integrate the random effects out. This means we can ask, across all possible values the farm random effect \(\mathbf{b}_i\) could take, what is the average probability of the data, weighted by how plausible each value of \(\mathbf{b}_i\) is?

\[\begin{equation} L(\boldsymbol{\beta}, G) = \prod_i \int \left[ \prod_j f(y_{ij} \mid \mathbf{b}_i) \right] f(\mathbf{b}_i \mid G) \, d\mathbf{b}_i \end{equation}\]

Reading this expression from the inside out we see that the term in brackets, \(\prod_j f(y_{ij} \mid \mathbf{b}_i)\), is exactly the conditional likelihood for farm \(i\) that we wrote down above. We then multiply it by \(f(\mathbf{b}_i \mid G)\), the Gaussian density of the random effect, which acts as a weight telling us how probable each value of \(\mathbf{b}_i\) is. The integral \(\int \ldots \, d\mathbf{b}_i\) averages over the entire range of possible random effect values, weighting each by that Gaussian density. Finally, the outer \(\prod_i\) multiplies this averaged quantity across all farms to give the full likelihood. Intuitively, we are saying: since we can’t observe the farm random effects, we account for all of them at once by averaging over their distribution.

This means that to fit our model, we cannot evaluate \(L(\boldsymbol{\beta}, G)\) exactly, and therefore cannot maximise it directly. Instead, we must either approximate the integral or find an alternative approach to estimation altogether. The following section covers the main strategies for doing so.

Estimation Procedures

Since we cannot evaluate the marginal likelihood \(L(\boldsymbol{\beta}, G)\) directly, we need alternative strategies. There are (arguably) three that dominate applied practice. The Laplace approximation replaces the integrand with a Gaussian centred at its mode and integrates that instead. Adaptive Gauss-Hermite quadrature (AGQ) evaluates the integral numerically at a set of points chosen around that same mode, with the Laplace approximation recovered as the one-point special case. MCMC sidesteps approximation altogether by sampling. A fourth method, penalized quasi-likelihood (PQL), has some historical importance for those interested in that, but is generally considered a last resort in modern practice and is treated in the FAQs, alongside other approaches such as variational inference and Monte Carlo EM.

For most applied work, we default to recommending the simplest option of going with the default (Laplace), which is is the default in both glmmTMB and lme4 and is accurate enough for the overwhelming majority of problems. The situation where it is worth reaching past it is a binary response with few observations per group, where the integrand is furthest from Gaussian and AGQ does provide some accuracy at modest computational cost. MCMC is the most flexible and the most expensive, and is rarely the right first move in a frequentist framework.

The derivations behind these methods, what the optimizer is actually doing when you call glmmTMB(), and why the resulting likelihood surface is such an awkward object to optimize over, are all treated in the appendix on estimation. None of it is required to fit a model, but we think it is worth reading eventually, because most of the warnings you will encounter in practice are the surface geometry being difficult and naughty.

The Data

Since we are working with a synthetic dataset, we have the unusual luxury of knowing the true parameter values.

The data-generating process mirrors the probability model we derived earlier in that each plant’s survival is a Bernoulli draw whose probability is determined by the inverse logit of a linear predictor combining an intercept, a genotype effect, an irrigation effect, and a farm-level random draw. The specific parameter values are visible in the code in the appendix.

plants <- readr::read_csv(
    here::here("./research-guides/glmms/assets/data/plants.csv"),
    show_col_types = FALSE
) |>
    dplyr::mutate(
        genotype = factor(genotype, levels = c("G1", "G2", "G3", "G4")),
        irrigation = factor(irrigation, levels = c("rainfed", "overhead", "drip")),
        farm = factor(farm)
    )

One structural feature of these data is worth flagging before we fit anything, because it has real consequences for inference. Irrigation method is assigned at the farm level — every plant on a given farm shares the same irrigation method — which means that the irrigation effect and the farm random effect are not fully separable without sufficient between-farm replication. This is not a quirk of the simulation; it reflects a genuine challenge in applied hierarchical data where a key predictor of interest lives at the group level rather than the observation level. We return to this point in the sample size discussion below.

ggplot(data = df) +
    geom_point(aes(x = x, y = y)) +
    # ... whatever else I need to plot
    theme_better() # and I add my theme to the end!

theme_foundation <- function(base_size = 20, base_family = "") {
    thm <- theme_grey(base_size = base_size, base_family = base_family)
    for (i in names(thm)) {
        if ("colour" %in% names(thm[[i]])) {
            thm[[i]]["colour"] <- list(NULL)
        }
        if ("fill" %in% names(thm[[i]])) {
            thm[[i]]["fill"] <- list(NULL)
        }
    }
    thm + theme(
        panel.border = element_rect(fill = NA),
        legend.background = element_rect(colour = NA),
        line = element_line(colour = "black"),
        rect = element_rect(fill = "white", colour = "black"),
        text = element_text(colour = "black")
    )
}

theme_better <- function(base_size = 20, base_family = "") {
    theme_foundation() +
        theme(
            line = element_line(
                colour = "black",
                lineend = "round",
                linetype = "solid"
            ),
            rect = element_rect(
                fill = "white",
                colour = "black",
                linetype = "solid"
            ),
            text = element_text(
                colour = "black",
                face = "plain",
                family = "Helvetica",
                size = base_size,
                vjust = 0.5,
                hjust = 0.5,
                lineheight = 1
            ),
            panel.grid = element_blank(),
            strip.background = element_rect(colour = NA),
            legend.key = element_rect(colour = NA),
            title = element_text(size = rel(1)),
            plot.title = element_text(
                size = rel(1.4), face = "bold",
                hjust = 0.5
            ),
            axis.title = element_text(
                size = rel(1.15),
                face = "bold"
            ),
            axis.text = element_text(
                size = rel(1)
            ),
            strip.text = element_text(),
            axis.ticks.length = unit(0.5, "lines"),
            plot.background = element_rect(fill = "transparent", colour = NA),
            panel.background = element_rect(fill = "transparent", colour = NA)
        )
}

library(ggplot2)
library(patchwork)

# a small helper to compute Wilson binomial confidence intervals —
# these behave better than Wald intervals near 0 and 1
wilson_ci <- function(x) {
    n <- length(x)
    phat <- mean(x)
    z <- qnorm(0.975)
    denom <- 1 + z^2 / n
    centre <- (phat + z^2 / (2 * n)) / denom
    half <- (z / denom) * sqrt(phat * (1 - phat) / n + z^2 / (4 * n^2))
    c(estimate = centre, lower = centre - half, upper = centre + half)
}

# --- plot 1: survival by genotype ---

geno_summary <- do.call(rbind, lapply(levels(plants$genotype), function(g) {
    ci <- wilson_ci(plants$survival[plants$genotype == g])
    data.frame(
        genotype = g,
        estimate = ci["estimate"],
        lower = ci["lower"],
        upper = ci["upper"]
    )
}))

p1 <- ggplot(geno_summary, aes(x = genotype, y = estimate)) +
    geom_errorbar(aes(ymin = lower, ymax = upper),
        width = 0, colour = "#2c7bb6", linewidth = 1.2
    ) +
    geom_point(size = 3.5, colour = "#2c7bb6") +
    scale_y_continuous(
        limits = c(0.45, 0.85),
        labels = function(x) paste0(round(x * 100), "%")
    ) +
    labs(
        x = "Genotype", y = "Observed survival rate",
        title = "Survival by genotype"
    ) +
    theme_better(base_size = 13) +
    theme(panel.grid.minor = element_blank())

# --- plot 2: survival by irrigation ---

irr_summary <- do.call(rbind, lapply(levels(plants$irrigation), function(m) {
    ci <- wilson_ci(plants$survival[plants$irrigation == m])
    data.frame(
        irrigation = m,
        estimate = ci["estimate"],
        lower = ci["lower"],
        upper = ci["upper"]
    )
}))
irr_summary$irrigation <- factor(irr_summary$irrigation,
    levels = c("rainfed", "overhead", "drip")
)

p2 <- ggplot(irr_summary, aes(x = irrigation, y = estimate)) +
    geom_errorbar(aes(ymin = lower, ymax = upper),
        width = 0, colour = "#d7191c", linewidth = 1.2
    ) +
    geom_point(size = 3.5, colour = "#d7191c") +
    scale_y_continuous(
        limits = c(0.45, 0.85),
        labels = function(x) paste0(round(x * 100), "%")
    ) +
    labs(
        x = "Irrigation method", y = "Observed survival rate",
        title = "Survival by irrigation method"
    ) +
    theme_better(base_size = 13) +
    theme(panel.grid.minor = element_blank())

# --- plot 3: farm-level survival rates coloured by irrigation ---

farm_summary <- do.call(rbind, lapply(levels(plants$farm), function(f) {
    sub <- plants[plants$farm == f, ]
    ci <- wilson_ci(sub$survival)
    data.frame(
        farm = f,
        irrigation = sub$irrigation[1],
        estimate = ci["estimate"],
        lower = ci["lower"],
        upper = ci["upper"]
    )
}))

# sort farms within each irrigation group by observed survival rate
farm_summary <- farm_summary[order(
    farm_summary$irrigation,
    farm_summary$estimate
), ]
farm_summary$farm_order <- seq_len(nrow(farm_summary))

irr_colours <- c(rainfed = "#d95f02", overhead = "#7570b3", drip = "#1b9e77")

p3 <- ggplot(
    farm_summary,
    aes(x = farm_order, y = estimate, colour = irrigation)
) +
    geom_errorbar(aes(ymin = lower, ymax = upper), width = 0, linewidth = 1.2) +
    geom_point(size = 2.5) +
    scale_colour_manual(values = irr_colours, name = "Irrigation") +
    scale_y_continuous(labels = function(x) paste0(round(x * 100), "%")) +
    labs(
        x = "Farm (sorted within irrigation group)",
        y = "Observed survival rate",
        title = "Farm-level survival rates"
    ) +
    theme_better(base_size = 13) +
    theme(
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        panel.grid.minor = element_blank()
    )

p1 / p2 / p3

Specifying the Model

With the data in hand, we can now write down our model in R. Recall what we are specifying: plant survival as a binary response, genotype and irrigation as fixed effects, and a random intercept for each farm to absorb the between-farm variation we saw in the plots above. In mathematical notation this is the model we already derived, \[\text{logit}(\mu_{ij}) = \beta_0 + \beta_{\text{geno}[ij]} + \beta_{\text{irr}[i]} + b_i,\] where \(b_i \sim \mathcal{N}(0, \sigma^2_{\text{farm}})\). Translating this into R syntax is, for the most part, pleasingly direct.

Before diving into the code, it is worth briefly introducing the two R packages we will be using throughout this section, and connecting each one back to the estimation procedures we discussed above. The R ecosystem for GLMMs is rich and occasionally bewildering, so having a clear sense of what each tool is doing under the hood — and why we have chosen it — will serve us well as we move through the implementation.

The workhorse of this guide is glmmTMB, developed by Brooks and colleagues. It uses the Template Model Builder (TMB) framework to fit models via Laplace approximation, and is our recommended default for most applied problems. If you would like to know what that framework is actually doing when you call it, we cover it in the appendix on estimation. It has a really nice combination of flexibility, speed, and reliability, and supports an exceptionally wide range of response distributions and random effects structures, and its formula interface is intuitive for anyone already familiar with base R modelling functions. Though we don’t talk about this much (though see the FAQs), it also provides (in our opinion) best-in-class support for zero-inflated models. For a thorough treatment of the package and its capabilities, we recommend the paper by Brooks et al. (2017), and a thorough perusal of the docs.

We also demonstrate lme4, the long-standing standard for mixed-effects models in R, developed by D. Bates et al. (2014). Like glmmTMB, it fits GLMMs via Laplace approximation (with the option of adaptive Gauss-Hermite quadrature for models with a single scalar random effect), and its formula syntax is essentially identical. For many straightforward problems lme4 and glmmTMB will give you nearly identical results, and lme4 remains extremely well-documented and widely understood in the applied literature. We include it here both for completeness and because it remains the most commonly encountered package for this class of models in published work.

Both packages share the same general convention for the fixed-effects part of the formula, which follows the standard R model formula syntax: response ~ predictor1 + predictor2. The key syntactic question is how each package handles the random effects, and this is where the packages diverge in ways that are worth understanding rather than just memorizing.

The (1 | farm) notation, shared by both glmmTMB and lme4, is read as “a random intercept (1) grouped by farm.” The 1 here is the same intercept term you would find in any R model formula, and the | operator separates what is varying (the left side) from what it is varying over (the right side, the grouping factor). A random slope, for comparison, would look like (1 + genotype | farm), specifying both a random intercept and a random slope for genotype within each farm. More complex structures — crossed random effects, nested groupings, and so on — are handled by extending this notation further, and we point the reader to the FAQs for a treatment of those cases.

fit_tmb <- glmmTMB::glmmTMB(
    survival ~ genotype + irrigation + (1 | farm),
    data = plants,
    family = binomial(link = "logit")
)

fit_lme4 <- lme4::glmer(
    survival ~ genotype + irrigation + (1 | farm),
    data = plants,
    family = binomial(link = "logit")
)

# MCMCglmm uses a different family argument and requires explicit priors —
# what follows is a minimal weakly informative specification
prior <- list(
    R = list(V = 1, fix = 1),
    G = list(G1 = list(V = 1, nu = 1, alpha.mu = 0, alpha.V = 1000))
)

fit_mcmc <- MCMCglmm::MCMCglmm(
    survival ~ genotype + irrigation,
    random = ~farm,
    family = "categorical",
    prior = prior,
    data = plants,
    verbose = FALSE
)

Fitting the Model

With the model specified, we can now actually run it! The code chunks below are identical in structure to what we saw in the specification section above, but this time we actually evaluate the code, so these are the fitted objects we will carry forward through the rest of the implementation.

fit_tmb <- glmmTMB::glmmTMB(
    survival ~ genotype + irrigation + (1 | farm),
    data = plants,
    family = binomial(link = "logit")
)

fit_lme4 <- lme4::glmer(
    survival ~ genotype + irrigation + (1 | farm),
    data = plants,
    family = binomial(link = "logit")
)

# running the chain — nitt controls total iterations, burnin discards
# early samples before the chain has converged, thin retains every nth sample
fit_mcmc <- MCMCglmm::MCMCglmm(
    survival ~ genotype + irrigation,
    random = ~farm,
    family = "categorical",
    prior = prior,
    nitt = 13000,
    burnin = 3000,
    thin = 10,
    data = plants,
    verbose = FALSE
)

summary(fit_tmb)

 Family: binomial  ( logit )
Formula:          survival ~ genotype + irrigation + (1 | farm)
Data: plants

      AIC       BIC    logLik -2*log(L)  df.resid 
   2548.1    2588.2   -1267.1    2534.1      2260 

Random effects:

Conditional model:
 Groups Name        Variance Std.Dev.
 farm   (Intercept) 0.3328   0.5769  
Number of obs: 2267, groups:  farm, 25

Conditional model:
                   Estimate Std. Error z value Pr(>|z|)    
(Intercept)          0.3333     0.1980   1.683 0.092281 .  
genotypeG2           0.7238     0.1487   4.866 1.14e-06 ***
genotypeG3          -0.6711     0.1317  -5.094 3.51e-07 ***
genotypeG4           0.3141     0.1386   2.266 0.023458 *  
irrigationoverhead   0.8697     0.3356   2.592 0.009553 ** 
irrigationdrip       0.9807     0.2890   3.394 0.000689 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(fit_lme4)

Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: binomial  ( logit )
Formula: survival ~ genotype + irrigation + (1 | farm)
   Data: plants

      AIC       BIC    logLik -2*log(L)  df.resid 
   2548.1    2588.2   -1267.1    2534.1      2260 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.0747 -0.8397  0.4187  0.6644  1.9894 

Random effects:
 Groups Name        Variance Std.Dev.
 farm   (Intercept) 0.3328   0.5768  
Number of obs: 2267, groups:  farm, 25

Fixed effects:
                   Estimate Std. Error z value Pr(>|z|)    
(Intercept)          0.3333     0.1979   1.684 0.092126 .  
genotypeG2           0.7238     0.1484   4.876 1.08e-06 ***
genotypeG3          -0.6711     0.1315  -5.101 3.37e-07 ***
genotypeG4           0.3141     0.1384   2.270 0.023213 *  
irrigationoverhead   0.8696     0.3354   2.593 0.009528 ** 
irrigationdrip       0.9807     0.2888   3.396 0.000684 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) gntyG2 gntyG3 gntyG4 irrgtnv
genotypeG2  -0.326                             
genotypeG3  -0.357  0.479                      
genotypeG4  -0.351  0.460  0.515               
irrigtnvrhd -0.490  0.014  0.003  0.024        
irrigatndrp -0.559  0.007 -0.018  0.003  0.331 

# simulate residuals from the fitted model — 250 simulations is usually
# sufficient for a dataset of this size
sim_res_tmb <- DHARMa::simulateResiduals(fit_tmb)
plot(sim_res_tmb)

DHARMa::testDispersion(sim_res_tmb)

    DHARMa nonparametric dispersion test via sd of residuals fitted vs.
    simulated

data:  simulationOutput
dispersion = 1.0058, p-value = 0.696
alternative hypothesis: two.sided

sim_res_lme4 <- DHARMa::simulateResiduals(fit_lme4)
plot(sim_res_lme4)

DHARMa::testDispersion(sim_res_lme4)

    DHARMa nonparametric dispersion test via sd of residuals fitted vs.
    simulated

data:  simulationOutput
dispersion = 1.0058, p-value = 0.696
alternative hypothesis: two.sided

Testing Fixed Effects

The p-values in the summary() output come from Wald z-tests, which are fast and usually fine as a first pass. However, for hypothesis testing — especially for categorical predictors with multiple levels like genotype, where you want to ask “does genotype matter at all?” rather than “is this particular contrast significant?” — likelihood ratio tests (LRTs) are generally preferred. An LRT compares the full model against a reduced model with the predictor of interest removed, and the test statistic is twice the difference in log-likelihoods, which follows a \(\chi^2\) distribution under the null. This is more principled than the Wald approach because it doesn’t rely on the asymptotic normality of the coefficient estimates in the same way, and it gives you a single omnibus test for each term rather than one test per level of a factor.

The easiest way to run LRTs for all fixed effects simultaneously is drop1() with test = "Chisq". This refits the model once for each term, each time dropping that term and comparing the fit to the full model. It’s a little slower than a Wald test since it requires actually refitting, but for a model this size it’s instantaneous.

drop1(fit_tmb, test = "Chisq")

Single term deletions

Model:
survival ~ genotype + irrigation + (1 | farm)
           Df    AIC     LRT  Pr(>Chi)    
<none>        2548.1                      
genotype    3 2654.2 112.134 < 2.2e-16 ***
irrigation  2 2555.3  11.153  0.003786 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

drop1(fit_lme4, test = "Chisq")

Single term deletions

Model:
survival ~ genotype + irrigation + (1 | farm)
           npar    AIC     LRT   Pr(Chi)    
<none>          2548.1                      
genotype      3 2654.2 112.134 < 2.2e-16 ***
irrigation    2 2555.3  11.153  0.003786 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

It is also worth knowing how to test whether the random intercept itself is earning its place. The natural approach is to fit a no-random-effect equivalent and compare the two models. The key here is to keep both models in the same package — mixing a base R glm() with a glmmTMB object will produce unreliable output, so we fit the reduced model using glmmTMB as well, simply omitting the (1 | farm) term.

# fit an equivalent model with no random effect — same package, same family
fit_tmb_no_re <- glmmTMB::glmmTMB(
    survival ~ genotype + irrigation,
    data = plants,
    family = binomial(link = "logit")
)

# AIC comparison — lower is better
AIC(fit_tmb_no_re, fit_tmb)

              df      AIC
fit_tmb_no_re  6 2633.904
fit_tmb        7 2548.106

# LRT via anova()
anova(fit_tmb_no_re, fit_tmb)

Data: plants
Models:
fit_tmb_no_re: survival ~ genotype + irrigation, zi=~0, disp=~1
fit_tmb: survival ~ genotype + irrigation + (1 | farm), zi=~0, disp=~1
              Df    AIC    BIC logLik deviance  Chisq Chi Df Pr(>Chisq)    
fit_tmb_no_re  6 2633.9 2668.3  -1311   2621.9                             
fit_tmb        7 2548.1 2588.2  -1267   2534.1 87.798      1  < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Both comparisons tell the same story, but they are answering slightly different questions. The AIC difference quantifies how much predictive information is lost by dropping the random effect, with a rule of thumb that differences greater than ~2 favour the more complex model. The LRT p-value comes from a \(\chi^2\) distribution, but with a caveat worth knowing about: testing whether a variance component equals zero is a boundary hypothesis, because variance can’t go negative — the null sits right on the edge of the allowable parameter space. The theoretical consequence is that the \(\chi^2\) null distribution is anti-conservative, and the true p-value is approximately twice the reported value, making the test conservative for this use case. In practice, if the LRT comes back significant you can trust it, but don’t read the precise p-value too literally. For most applied purposes, a large AIC improvement combined with a visual inspection of the farm-level spread (as in our earlier plot) is sufficient justification for keeping the random effect.

Visualizing Model Results

Numbers in a summary table are all well and good, but a model is rarely truly understood until its results have been turned into something you can actually look at. In this section we’ll work through a set of visualizations that collectively tell us what our fitted model is saying: where the coefficients land relative to the truth we built in, what the model predicts on the probability scale, and how much the farms vary around the population average. Before we get there, though, we mustn’t skip over some pieces of interpretation that often trip up folks when working on Bernoulli data for the first time.

# what are our raw estimates?
glmmTMB::fixef(fit_tmb)$cond

       (Intercept)         genotypeG2         genotypeG3         genotypeG4 
         0.3333242          0.7237880         -0.6710730          0.3141378 
irrigationoverhead     irrigationdrip 
         0.8697466          0.9807044 

# exponentiate to get odds ratios
exp(glmmTMB::fixef(fit_tmb)$cond)

       (Intercept)         genotypeG2         genotypeG3         genotypeG4 
         1.3955996          2.0622302          0.5111598          1.3690784 
irrigationoverhead     irrigationdrip 
         2.3863061          2.6663336 

Odds ratios are better than log-odds, but they’re still not probabilities, and probabilities are what most people want. The inverse logit function plogis() converts a value on the log-odds scale directly to a probability. So if we want to know the predicted survival probability for a G1 plant on a rainfed farm with a zero random effect (i.e., an “average” farm), we can simply do:

# extract the intercept and back-transform it to the probability scale
intercept_est <- glmmTMB::fixef(fit_tmb)$cond["(Intercept)"]
plogis(intercept_est)

(Intercept) 
   0.582568 

And to get the predicted probability for, say, a G2 plant under the same conditions, we add the G2 coefficient before back-transforming — because the addition happens on the log-odds scale, not the probability scale. This is the key thing to remember: you must not add effects on the probability scale. Always add on the log-odds scale first, then back-transform.

# add the G2 effect on the log-odds scale, then back-transform
beta <- glmmTMB::fixef(fit_tmb)$cond
plogis(beta["(Intercept)"] + beta["genotypeG2"])

(Intercept) 
  0.7421383 

We will use exactly this logic — construct a linear predictor on the log-odds scale and then apply plogis() — to generate all of the predicted probability plots below.

Here are the true values (found in Section 7.5) that we’ll want for our results plots:

# true vals for plotting
intercept <- 0.5
beta_overhead <- 0.20
beta_drip <- 0.80
beta_G2 <- 0.50
beta_G3 <- -0.75
beta_G4 <- 0.20
farm_re_sd <- 0.8

irrigation_types <- c("rainfed", "overhead", "drip")

# true vals for the count model, same idea as the plant model above
insect_intercept <- 1.2
insect_beta_G2 <- 0.40
insect_beta_G3 <- -0.30
insect_beta_G4 <- 0.15
insect_beta_overhead <- 0.25
insect_beta_drip <- -0.50
insect_farm_re_sd <- 0.5
insect_theta <- 2
insect_zi_intercept <- -2.2
insect_zi_drip <- 1.50

# extract the fixed effects summary and build a data frame for plotting
cs_tmb <- coef(summary(fit_tmb))$cond
coef_tmb <- data.frame(
    param = rownames(cs_tmb),
    estimate = cs_tmb[, "Estimate"],
    lower = cs_tmb[, "Estimate"] - 1.96 * cs_tmb[, "Std. Error"],
    upper = cs_tmb[, "Estimate"] + 1.96 * cs_tmb[, "Std. Error"]
)

# build a matching data frame of true values to overlay
true_tmb <- data.frame(
    param = c(
        "(Intercept)", "genotypeG2", "genotypeG3", "genotypeG4",
        "irrigationoverhead", "irrigationdrip"
    ),
    true_val = c(intercept, beta_G2, beta_G3, beta_G4, beta_overhead, beta_drip)
)
coef_tmb <- merge(coef_tmb, true_tmb, by = "param")

coef_tmb$param <- factor(coef_tmb$param,
    levels = rev(c(
        "(Intercept)", "genotypeG2", "genotypeG3", "genotypeG4",
        "irrigationoverhead", "irrigationdrip"
    ))
)

ggplot(coef_tmb, aes(x = estimate, y = param)) +
    geom_errorbar(aes(xmin = lower, xmax = upper, y = param),
        width = 0, linewidth = 0.9
    ) +
    geom_point(size = 3) +
    geom_point(aes(x = true_val),
        position = position_nudge(y = 0.2),
        colour = "#d7191c", shape = 18, size = 4
    ) +
    geom_vline(xintercept = 0, linetype = "dashed", colour = "grey50") +
    labs(x = "Estimate (log-odds scale)", y = NULL) +
    theme_better(base_size = 13)

cs_lme4 <- coef(summary(fit_lme4))
coef_lme4 <- data.frame(
    param = rownames(cs_lme4),
    estimate = cs_lme4[, "Estimate"],
    lower = cs_lme4[, "Estimate"] - 1.96 * cs_lme4[, "Std. Error"],
    upper = cs_lme4[, "Estimate"] + 1.96 * cs_lme4[, "Std. Error"]
)

true_lme4 <- data.frame(
    param = c(
        "(Intercept)", "genotypeG2", "genotypeG3", "genotypeG4",
        "irrigationoverhead", "irrigationdrip"
    ),
    true_val = c(intercept, beta_G2, beta_G3, beta_G4, beta_overhead, beta_drip)
)
coef_lme4 <- merge(coef_lme4, true_lme4, by = "param")

coef_lme4$param <- factor(coef_lme4$param,
    levels = rev(c(
        "(Intercept)", "genotypeG2", "genotypeG3", "genotypeG4",
        "irrigationoverhead", "irrigationdrip"
    ))
)

ggplot(coef_lme4, aes(x = estimate, y = param)) +
    geom_errorbar(aes(xmin = lower, xmax = upper, y = param),
        width = 0, linewidth = 0.9
    ) +
    geom_point(size = 3) +
    geom_point(aes(x = true_val),
        position = position_nudge(y = 0.2),
        colour = "#d7191c", shape = 18, size = 4
    ) +
    geom_vline(xintercept = 0, linetype = "dashed", colour = "grey50") +
    labs(x = "Estimate (log-odds scale)", y = NULL) +
    theme_better(base_size = 13)

# regenerate data at a much larger scale — same true parameters, more farms
set.seed(999)

n_farms_big <- 1000
plants_per_farm_big <- rep(500, n_farms_big)
farm_irrigation_big <- sample(irrigation_types, n_farms_big, replace = TRUE)
farm_re_big <- rnorm(n_farms_big, mean = 0, sd = farm_re_sd)

plant_list_big <- lapply(seq_len(n_farms_big), function(i) {
    n <- plants_per_farm_big[i]
    genotype <- sample(c("G1", "G2", "G3", "G4"), n, replace = TRUE)
    irr_effect <- switch(farm_irrigation_big[i],
        overhead = beta_overhead,
        drip = beta_drip,
        0
    )
    geno_effect <- ifelse(genotype == "G2", beta_G2,
        ifelse(genotype == "G3", beta_G3,
            ifelse(genotype == "G4", beta_G4, 0)
        )
    )
    eta <- intercept + irr_effect + geno_effect + farm_re_big[i]
    p <- plogis(eta)
    survival <- rbinom(n, size = 1, prob = p)
    data.frame(
        farm = paste0("farm_", i),
        irrigation = farm_irrigation_big[i],
        genotype = genotype,
        survival = survival
    )
})

plants_big <- do.call(rbind, plant_list_big)
plants_big$genotype <- factor(plants_big$genotype, levels = c("G1", "G2", "G3", "G4"))
plants_big$irrigation <- factor(plants_big$irrigation, levels = c("rainfed", "overhead", "drip"))
plants_big$farm <- factor(plants_big$farm)

fit_tmb_big <- glmmTMB::glmmTMB(
    survival ~ genotype + irrigation + (1 | farm),
    data = plants_big,
    family = binomial(link = "logit")
)

cs_big <- coef(summary(fit_tmb_big))$cond
coef_big <- data.frame(
    param = rownames(cs_big),
    estimate = cs_big[, "Estimate"],
    lower = cs_big[, "Estimate"] - 1.96 * cs_big[, "Std. Error"],
    upper = cs_big[, "Estimate"] + 1.96 * cs_big[, "Std. Error"]
)
true_vals <- data.frame(
    param = c(
        "(Intercept)", "genotypeG2", "genotypeG3", "genotypeG4",
        "irrigationoverhead", "irrigationdrip"
    ),
    true_val = c(intercept, beta_G2, beta_G3, beta_G4, beta_overhead, beta_drip)
)
coef_big <- merge(coef_big, true_vals, by = "param")
coef_big$param <- factor(coef_big$param,
    levels = rev(c(
        "(Intercept)", "genotypeG2", "genotypeG3", "genotypeG4",
        "irrigationoverhead", "irrigationdrip"
    ))
)

ggplot(coef_big, aes(x = estimate, y = param)) +
    geom_errorbar(aes(xmin = lower, xmax = upper),
        width = 0, linewidth = 0.9
    ) +
    geom_point(size = 3) +
    geom_point(aes(x = true_val),
        position = position_nudge(y = 0.2),
        colour = "#d7191c", shape = 18, size = 4
    ) +
    geom_vline(xintercept = 0, linetype = "dashed", colour = "grey50") +
    labs(x = "Estimate (log-odds scale)", y = NULL) +
    theme_better(base_size = 13)

# the conditional prediction sets the random effect to zero. the marginal one
# averages the inverse link over the random effect distribution, which is a
# different quantity and not just a noisier version of the same one
beta <- glmmTMB::fixef(fit_tmb)$cond
sigma_farm <- sqrt(glmmTMB::VarCorr(fit_tmb)$cond$farm[1])

# a G2 plant on a drip-irrigated farm, our highest-probability cell
eta <- beta["(Intercept)"] + beta["genotypeG2"] + beta["irrigationdrip"]

set.seed(2024)
b_draws <- rnorm(2e5, mean = 0, sd = sigma_farm)

c(
    conditional = unname(plogis(eta)),
    marginal = unname(mean(plogis(eta + b_draws)))
)

conditional    marginal 
  0.8847107   0.8720378 

# i sweep the random effect sd and recompute the marginal probability at each
# value, holding the fixed effects fixed. the conditional prediction does not
# move at all, because it never sees the random effect distribution
sigma_grid <- seq(0, 3, length.out = 80)

set.seed(2024)
marginal_curve <- sapply(sigma_grid, function(s) {
    mean(plogis(eta + rnorm(5e4, mean = 0, sd = s)))
})

sweep_df <- data.frame(sigma = sigma_grid, marginal = marginal_curve)

ggplot(sweep_df, aes(x = sigma, y = marginal)) +
    geom_hline(
        yintercept = plogis(eta),
        linetype = "dashed",
        colour = "grey40"
    ) +
    geom_vline(
        xintercept = sigma_farm,
        linetype = "dotted",
        colour = "#B4656F"
    ) +
    geom_line(linewidth = 1, colour = "#1D9E75") +
    scale_y_continuous(limits = c(0.5, NA)) +
    labs(
        x = expression(paste("Random effect SD, ", sigma[farm])),
        y = "Predicted survival probability"
    ) +
    theme_better(base_size = 13)

# build prediction grids for genotype and irrigation separately
geno_grid <- data.frame(
    genotype = factor(levels(plants$genotype), levels = levels(plants$genotype)),
    irrigation = factor("rainfed", levels = levels(plants$irrigation))
)
irr_grid <- data.frame(
    genotype = factor("G1", levels = levels(plants$genotype)),
    irrigation = factor(levels(plants$irrigation), levels = levels(plants$irrigation))
)

# use predict() with se.fit = TRUE to get link-scale SEs, then back-transform
geno_link <- predict(fit_tmb,
    newdata = geno_grid, type = "link",
    se.fit = TRUE, re.form = NA
)
geno_grid$estimate <- plogis(geno_link$fit)
geno_grid$lower <- plogis(geno_link$fit - 1.96 * geno_link$se.fit)
geno_grid$upper <- plogis(geno_link$fit + 1.96 * geno_link$se.fit)

irr_link <- predict(fit_tmb,
    newdata = irr_grid, type = "link",
    se.fit = TRUE, re.form = NA
)
irr_grid$estimate <- plogis(irr_link$fit)
irr_grid$lower <- plogis(irr_link$fit - 1.96 * irr_link$se.fit)
irr_grid$upper <- plogis(irr_link$fit + 1.96 * irr_link$se.fit)

# compute true marginal probabilities at re = 0 for overlay
true_geno <- data.frame(
    genotype = factor(levels(plants$genotype),
        levels = levels(plants$genotype)
    ),
    true_prob = plogis(c(
        intercept,
        intercept + beta_G2,
        intercept + beta_G3,
        intercept + beta_G4
    ))
)
true_irr <- data.frame(
    irrigation = factor(levels(plants$irrigation),
        levels = levels(plants$irrigation)
    ),
    true_prob = plogis(c(
        intercept,
        intercept + beta_overhead,
        intercept + beta_drip
    ))
)

p1 <- ggplot(geno_grid, aes(x = genotype, y = estimate)) +
    geom_errorbar(aes(ymin = lower, ymax = upper),
        width = 0, linewidth = 1.2, colour = "#2c7bb6"
    ) +
    geom_point(size = 3.5, colour = "#2c7bb6") +
    geom_point(
        data = true_geno, aes(y = true_prob),
        position = position_nudge(x = 0.2),
        shape = 18, size = 4, colour = "#d7191c"
    ) +
    scale_y_continuous(
        limits = c(0, 1),
        labels = function(x) paste0(round(x * 100), "%")
    ) +
    labs(x = "Genotype", y = "Predicted survival probability") +
    theme_better(base_size = 13)

p2 <- ggplot(irr_grid, aes(x = irrigation, y = estimate)) +
    geom_errorbar(aes(ymin = lower, ymax = upper),
        width = 0, linewidth = 1.2, colour = "#d7191c"
    ) +
    geom_point(size = 3.5, colour = "#d7191c") +
    geom_point(
        data = true_irr, aes(y = true_prob),
        position = position_nudge(x = 0.2),
        shape = 18, size = 4, colour = "#2c7bb6"
    ) +
    scale_y_continuous(
        limits = c(0, 1),
        labels = function(x) paste0(round(x * 100), "%")
    ) +
    labs(x = "Irrigation method", y = "Predicted survival probability") +
    theme_better(base_size = 13)

p1 / p2

# lme4's predict.merMod does not support se.fit, so we use the design matrix
# directly and propagate uncertainty via the delta method
predict_marginal_lme4 <- function(grid) {
    X <- model.matrix(~ genotype + irrigation, data = grid)
    beta <- lme4::fixef(fit_lme4)
    V <- as.matrix(vcov(fit_lme4))
    eta <- as.vector(X %*% beta)
    se_eta <- sqrt(diag(X %*% V %*% t(X)))
    grid$estimate <- plogis(eta)
    grid$lower <- plogis(eta - 1.96 * se_eta)
    grid$upper <- plogis(eta + 1.96 * se_eta)
    grid
}

geno_pred <- predict_marginal_lme4(geno_grid)
irr_pred <- predict_marginal_lme4(irr_grid)

p1 <- ggplot(geno_pred, aes(x = genotype, y = estimate)) +
    geom_errorbar(aes(ymin = lower, ymax = upper),
        width = 0, linewidth = 1.2, colour = "#2c7bb6"
    ) +
    geom_point(size = 3.5, colour = "#2c7bb6") +
    geom_point(
        data = true_geno, aes(y = true_prob),
        position = position_nudge(x = 0.2),
        shape = 18, size = 4, colour = "#d7191c"
    ) +
    scale_y_continuous(
        limits = c(0, 1),
        labels = function(x) paste0(round(x * 100), "%")
    ) +
    labs(x = "Genotype", y = "Predicted survival probability") +
    theme_better(base_size = 13)

p2 <- ggplot(irr_pred, aes(x = irrigation, y = estimate)) +
    geom_errorbar(aes(ymin = lower, ymax = upper),
        width = 0, linewidth = 1.2, colour = "#d7191c"
    ) +
    geom_point(size = 3.5, colour = "#d7191c") +
    geom_point(
        data = true_irr, aes(y = true_prob),
        position = position_nudge(x = 0.2),
        shape = 18, size = 4, colour = "#2c7bb6"
    ) +
    scale_y_continuous(
        limits = c(0, 1),
        labels = function(x) paste0(round(x * 100), "%")
    ) +
    labs(x = "Irrigation method", y = "Predicted survival probability") +
    theme_better(base_size = 13)

p1 / p2

# get the unique farm-irrigation mapping from the data
farm_info <- unique(plants[, c("farm", "irrigation")])

# build a grid of all genotype x farm combinations and attach irrigation
cond_grid <- merge(
    expand.grid(
        genotype = levels(plants$genotype),
        farm = levels(plants$farm)
    ),
    farm_info,
    by = "farm"
)

cond_grid$pred <- predict(
    fit_tmb,
    newdata = cond_grid,
    type = "response",
    re.form = NULL
)

ggplot(cond_grid, aes(
    x = genotype, y = pred,
    group = farm, colour = irrigation
)) +
    geom_line(alpha = 0.6, linewidth = 0.7) +
    scale_colour_manual(values = irr_colours, name = "Irrigation") +
    scale_y_continuous(
        limits = c(0, 1),
        labels = function(x) paste0(round(x * 100), "%")
    ) +
    labs(x = "Genotype", y = "Predicted survival probability") +
    theme_better(base_size = 13)

cond_grid$pred <- predict(
    fit_lme4,
    newdata = cond_grid,
    type = "response",
    re.form = NULL
)

ggplot(cond_grid, aes(
    x = genotype, y = pred,
    group = farm, colour = irrigation
)) +
    geom_line(alpha = 0.6, linewidth = 0.7) +
    scale_colour_manual(values = irr_colours, name = "Irrigation") +
    scale_y_continuous(
        limits = c(0, 1),
        labels = function(x) paste0(round(x * 100), "%")
    ) +
    labs(x = "Genotype", y = "Predicted survival probability") +
    theme_better(base_size = 13)

# extract the BLUPs from the glmmTMB fit and merge in irrigation for colour
re_vals <- glmmTMB::ranef(fit_tmb)$cond$farm
re_df <- data.frame(
    farm = rownames(re_vals),
    blup = re_vals[, "(Intercept)"]
)
farm_info <- unique(plants[, c("farm", "irrigation")])
re_df <- merge(re_df, farm_info, by = "farm")
re_df <- re_df[order(re_df$blup), ]
re_df$farm_order <- seq_len(nrow(re_df))

ggplot(re_df, aes(x = farm_order, y = blup, colour = irrigation)) +
    geom_hline(yintercept = 0, linetype = "dashed", colour = "grey50") +
    geom_point(size = 2.5) +
    scale_colour_manual(values = irr_colours, name = "Irrigation") +
    scale_y_continuous(limits = c(-2, 2)) +
    labs(
        x = "Farm (sorted by estimated random intercept)",
        y = "Estimated random intercept (log-odds scale)"
    ) +
    theme_better(base_size = 13) +
    theme(
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank()
    )

The Same Model in Julia

Everything above has been R. That’s a defensible default for mixed models, but it isn’t the only option, and MixedModels.jl is genuinely excellent — fast, actively developed, and written by Douglas Bates, who also wrote lme4. The formula syntax is close enough that most of what you’ve learned transfers directly.

We fit the identical model here: same data, same fixed effects, same random intercept. The point is to show you what the translation looks like, not to argue for one language over the other.

using CSV, DataFrames, CategoricalArrays, MixedModels, StatsModels

# reading the same csv the r code above reads. this is deliberate --- the two
# languages share a file rather than each simulating their own data, because
# matching random number streams across languages is not a thing you can do
plants = CSV.read(
    joinpath(@__DIR__, "assets", "data", "plants.csv"),
    DataFrame
);

first(plants, 5)

5×4 DataFrame
 Row │ farm     irrigation  genotype  survival
     │ String7  String15    String3   Int64
─────┼─────────────────────────────────────────
   1 │ farm_1   drip        G1               1
   2 │ farm_1   drip        G4               1
   3 │ farm_1   drip        G3               1
   4 │ farm_1   drip        G3               1
   5 │ farm_1   drip        G3               0

One thing to get right before fitting. R decides the reference level of a factor by alphabetical order unless you tell it otherwise; Julia uses the first level of a CategoricalArray. Since we set the levels explicitly in the R code, we do the same here so the coefficients are directly comparable.

plants.genotype = categorical(plants.genotype; levels = ["G1", "G2", "G3", "G4"]);
plants.irrigation = categorical(plants.irrigation; levels = ["rainfed", "overhead", "drip"]);
plants.farm = categorical(plants.farm);

Now the model. The formula is nearly identical to the R version, with @formula wrapping it and the response distribution passed as a positional argument rather than through family.

form = @formula(survival ~ 1 + genotype + irrigation + (1 | farm));

gm = fit(MixedModel, form, plants, Bernoulli())

Generalized Linear Mixed Model fit by maximum likelihood (nAGQ = 1)
  survival ~ 1 + genotype + irrigation + (1 | farm)
  Distribution: Bernoulli{Float64}
  Link: LogitLink()

   logLik    deviance     AIC       AICc        BIC    
 -1267.0529  2534.1058  2548.1058  2548.1553  2588.1893

Variance components:
        Column   Variance Std.Dev. 
farm (Intercept)  0.332777 0.576868

 Number of obs: 2267; levels of grouping factors: 25

Fixed-effects parameters:
────────────────────────────────────────────────────────────
                          Coef.  Std. Error      z  Pr(>|z|)
────────────────────────────────────────────────────────────
(Intercept)            0.333338    0.197688   1.69    0.0918
genotype: G2           0.723789    0.148269   4.88    <1e-05
genotype: G3          -0.671094    0.131045  -5.12    <1e-06
genotype: G4           0.314116    0.13809    2.27    0.0229
irrigation: overhead   0.869689    0.335245   2.59    0.0095
irrigation: drip       0.980707    0.288525   3.40    0.0007
────────────────────────────────────────────────────────────

Bernoulli() implies a logit link, so we don’t need to say so, though you can pass LogitLink() explicitly if you prefer being obvious about it. As in lme4 and glmmTMB, estimation is by Laplace approximation. If you want adaptive quadrature instead, nAGQ is available for a single scalar random effect, exactly the restriction lme4 imposes and for the same reason (see the appendix on estimation).

The variance component and the fixed effects come out through functions you’ll recognize:

VarCorr(gm)

Variance components:
        Column   Variance Std.Dev. 
farm (Intercept)  0.332777 0.576868

DataFrame(term = coefnames(gm), estimate = fixef(gm), se = stderror(gm))

6×3 DataFrame
 Row │ term                  estimate   se
     │ String                Float64    Float64
─────┼───────────────────────────────────────────
   1 │ (Intercept)            0.333338  0.197688
   2 │ genotype: G2           0.723789  0.148269
   3 │ genotype: G3          -0.671094  0.131045
   4 │ genotype: G4           0.314116  0.13809
   5 │ irrigation: overhead   0.869689  0.335245
   6 │ irrigation: drip       0.980707  0.288525

Compare these against the glmmTMB output from earlier and they should agree to several decimal places. That agreement is worth more than it might seem: the two packages were written independently, by different people, in different languages, and both are approximating the same intractable integral. When they land in the same place, that is real evidence neither has a bug in it. When they don’t, something is wrong and it’s worth finding out what.

Overdispersion and Zero Inflation

Everything above has worked with a binary response, which is convenient for exposition but hides two of the most common problems in applied GLMM work. Neither can be demonstrated with Bernoulli data, and the reason is structural rather than incidental: for a Bernoulli response the variance is \(p(1-p)\), entirely determined by the mean, so there is no free variance parameter for the data to exceed. Overdispersion at the observation level is not merely absent from our plant data, it is impossible.

Both problems live in count data, so for these two sections we switch to a companion dataset. While the plants were growing, we also counted the insects found on each one. Same plants, same farms, same genotypes and irrigation methods, but now the response is a count rather than a 0 or 1. The simulation code is in the appendix.

head(insects)

# A tibble: 6 × 4
  farm   irrigation genotype count
  <fct>  <fct>      <fct>    <dbl>
1 farm_1 drip       G1           1
2 farm_1 drip       G4           2
3 farm_1 drip       G3           0
4 farm_1 drip       G3           0
5 farm_1 drip       G3           0
6 farm_1 drip       G2           0

table(insects$count == 0)

FALSE  TRUE 
 1542   725 

Roughly a fifth of plants have no insects at all, which is our first hint that something beyond a simple Poisson model may be needed. Let us fit the obvious model and see what happens.

fit_pois <- glmmTMB::glmmTMB(
    count ~ genotype + irrigation + (1 | farm),
    data = insects,
    family = poisson(link = "log")
)
summary(fit_pois)

 Family: poisson  ( log )
Formula:          count ~ genotype + irrigation + (1 | farm)
Data: insects

      AIC       BIC    logLik -2*log(L)  df.resid 
  14347.1   14387.2   -7166.6   14333.1      2260 

Random effects:

Conditional model:
 Groups Name        Variance Std.Dev.
 farm   (Intercept) 0.4105   0.6407  
Number of obs: 2267, groups:  farm, 25

Conditional model:
                   Estimate Std. Error z value Pr(>|z|)    
(Intercept)         1.23957    0.18684   6.634 3.26e-11 ***
genotypeG2          0.33874    0.02907  11.651  < 2e-16 ***
genotypeG3         -0.31461    0.03427  -9.180  < 2e-16 ***
genotypeG4          0.14932    0.03019   4.946 7.57e-07 ***
irrigationoverhead  0.24814    0.34233   0.725  0.46855    
irrigationdrip     -0.87085    0.29477  -2.954  0.00313 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The model runs without complaint. This is worth sitting with for a moment, because it is the central practical lesson of these two sections. A Poisson GLMM will fit overdispersed, zero-inflated data quite happily, converge cleanly, and hand you a coefficient table with reassuringly small standard errors. Nothing in the output tells you it is wrong. You have to go looking.

sim_pois <- DHARMa::simulateResiduals(fit_pois)
DHARMa::testDispersion(sim_pois)

    DHARMa nonparametric dispersion test via sd of residuals fitted vs.
    simulated

data:  simulationOutput
dispersion = 2.8918, p-value < 2.2e-16
alternative hypothesis: two.sided

fit_nb <- glmmTMB::glmmTMB(
    count ~ genotype + irrigation + (1 | farm),
    data = insects,
    family = glmmTMB::nbinom2(link = "log")
)

AIC(fit_pois, fit_nb)

         df      AIC
fit_pois  7 14347.12
fit_nb    8 10301.18

insects$obs <- factor(seq_len(nrow(insects)))
fit_olre <- glmmTMB::glmmTMB(
    count ~ genotype + irrigation + (1 | farm) + (1 | obs),
    data = insects,
    family = poisson(link = "log")
)

sim_nb <- DHARMa::simulateResiduals(fit_nb)
DHARMa::testZeroInflation(sim_nb)

    DHARMa zero-inflation test via comparison to expected zeros with
    simulation under H0 = fitted model

data:  simulationOutput
ratioObsSim = 1.1084, p-value < 2.2e-16
alternative hypothesis: two.sided

fit_zi <- glmmTMB::glmmTMB(
    count ~ genotype + irrigation + (1 | farm),
    ziformula = ~irrigation,
    data = insects,
    family = glmmTMB::nbinom2(link = "log")
)

AIC(fit_nb, fit_zi)

       df      AIC
fit_nb  8 10301.18
fit_zi 11 10182.90

summary(fit_zi)

 Family: nbinom2  ( log )
Formula:          count ~ genotype + irrigation + (1 | farm)
Zero inflation:         ~irrigation
Data: insects

      AIC       BIC    logLik -2*log(L)  df.resid 
  10182.9   10245.9   -5080.5   10160.9      2256 

Random effects:

Conditional model:
 Groups Name        Variance Std.Dev.
 farm   (Intercept) 0.3905   0.6249  
Number of obs: 2267, groups:  farm, 25

Dispersion parameter for nbinom2 family ():  1.8 

Conditional model:
                   Estimate Std. Error z value Pr(>|z|)    
(Intercept)         1.35880    0.18761   7.243 4.40e-13 ***
genotypeG2          0.34014    0.06273   5.422 5.88e-08 ***
genotypeG3         -0.29693    0.06623  -4.483 7.36e-06 ***
genotypeG4          0.13046    0.06292   2.073   0.0381 *  
irrigationoverhead  0.23836    0.33755   0.706   0.4801    
irrigationdrip     -0.56283    0.29261  -1.923   0.0544 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Zero-inflation model:
                   Estimate Std. Error z value Pr(>|z|)    
(Intercept)         -2.0870     0.1645 -12.690  < 2e-16 ***
irrigationoverhead  -0.2695     0.2900  -0.929    0.353    
irrigationdrip       1.4318     0.1929   7.423 1.14e-13 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

fit_zi_intercept <- glmmTMB::glmmTMB(
    count ~ genotype + irrigation + (1 | farm),
    ziformula = ~1,
    data = insects,
    family = glmmTMB::nbinom2(link = "log")
)

AIC(fit_nb, fit_zi_intercept, fit_zi)

                 df      AIC
fit_nb            8 10301.18
fit_zi_intercept  9 10244.43
fit_zi           11 10182.90

# pull the conditional fixed effects out of both fits and stack them so the
# poisson and the corrected model can sit side by side on the same axis
coef_frame <- function(fit, label) {
    tab <- summary(fit)$coefficients$cond
    data.frame(
        term = rownames(tab),
        estimate = tab[, "Estimate"],
        se = tab[, "Std. Error"],
        model = label
    )
}

count_coefs <- rbind(
    coef_frame(fit_pois, "Poisson"),
    coef_frame(fit_zi, "NB + zero-inflation")
)
count_coefs$lower <- count_coefs$estimate - 1.96 * count_coefs$se
count_coefs$upper <- count_coefs$estimate + 1.96 * count_coefs$se

true_vals <- data.frame(
    term = c(
        "(Intercept)", "genotypeG2", "genotypeG3", "genotypeG4",
        "irrigationoverhead", "irrigationdrip"
    ),
    truth = c(
        insect_intercept, insect_beta_G2, insect_beta_G3, insect_beta_G4,
        insect_beta_overhead, insect_beta_drip
    )
)

count_coefs$term <- factor(count_coefs$term, levels = rev(true_vals$term))
true_vals$term <- factor(true_vals$term, levels = rev(true_vals$term))
count_coefs$model <- factor(count_coefs$model,
    levels = c("Poisson", "NB + zero-inflation")
)

ggplot(count_coefs, aes(x = estimate, y = term, colour = model)) +
    geom_vline(xintercept = 0, linetype = "dashed", colour = "grey50") +
    geom_errorbarh(
        aes(xmin = lower, xmax = upper),
        height = 0,
        linewidth = 0.7,
        position = position_dodge(width = 0.5)
    ) +
    geom_point(size = 2.5, position = position_dodge(width = 0.5)) +
    geom_point(
        data = true_vals,
        aes(x = truth, y = term),
        inherit.aes = FALSE,
        shape = 18,
        size = 3.5,
        colour = "red",
        position = position_nudge(y = 0.32)
    ) +
    scale_colour_manual(values = c(
        "Poisson" = "#B4656F",
        "NB + zero-inflation" = "#1D9E75"
    ), name = NULL) +
    labs(x = "Estimate (log scale)", y = NULL) +
    theme_better(base_size = 13) +
    theme(legend.position = "bottom")

Warning: `geom_errorbarh()` was deprecated in ggplot2 4.0.0.
ℹ Please use the `orientation` argument of `geom_errorbar()` instead.

`height` was translated to `width`.

What distribution should I use for my response?

The right choice of response distribution always comes back to the data-generating process: what distribution would your response variable plausibly follow if you had many replicates? The table below covers the most common cases. DHARMa residual diagnostics are a reliable way to check whether a chosen family is a good fit regardless of which you pick.

A note on overdispersion and zero-inflation: these are distinct problems and are worth distinguishing before reaching for a more complex model. Overdispersion is a general variance problem — the observed variance exceeds what the distribution predicts — and is addressed by switching to a negative binomial family. Zero-inflation is specific to count data where some zeros arise from a structural source (a species genuinely absent from a site, a patient who cannot experience the event) rather than from random chance. DHARMa::testDispersion() and DHARMa::testZeroInflation() are useful first diagnostics for each. For a fuller treatment of count model selection, we address both (though not hugely in-depth) in: Section 3.7.1 and Section 3.7.2

Can I add random slopes to my model?

Yes, and it is often a good idea when you have good reason to expect that the effect of a predictor varies across groups. A random slope allows each group its own slope for a given predictor, in addition to its own intercept.

# random intercept and slope for x, by group
fit_rs <- glmmTMB::glmmTMB(
    response ~ x + (1 + x | group),
    data = mydata,
    family = binomial()
)

The (1 + x | group) term estimates three parameters: the variance of the random intercepts, the variance of the random slopes, and the correlation between them. This is a substantially more complex model than a random-intercept-only specification, and it requires considerably more data to estimate reliably — at least 10–15 groups with enough observations per group to inform both the intercept and slope. Random slope models fail to converge more frequently, and singular fit warnings are especially common with small datasets. Before adding a random slope, ask whether the research question genuinely requires group-specific slopes, or whether a fixed-effect interaction would answer the question more parsimoniously. If you want a random slope without the intercept-slope correlation, use (1 | group) + (0 + x | group).

How do I specify crossed vs nested random effects?

This is probably the most common source of confusion in the (1 | ...) notation, and the distinction matters because the two structures answer very different questions about your data.

Nested random effects arise when lower-level units are uniquely embedded within higher-level units — plots within fields, students within classrooms. If plot labels are not unique across fields (plot 1 in field A is a different place from plot 1 in field B), you need to tell the model that explicitly:

# nested: i specify plots within fields using the / operator
fit_nested <- glmmTMB::glmmTMB(
    response ~ predictor + (1 | field/plot),
    data = mydata,
    family = binomial()
)

Crossed random effects arise when every combination of two grouping factors is (at least potentially) observed — every patient might be measured by multiple observers, and every observer measures multiple patients. Here the two random effects are fully crossed and neither is nested in the other:

# crossed: i add two separate random intercept terms
fit_crossed <- glmmTMB::glmmTMB(
    response ~ predictor + (1 | patient) + (1 | observer),
    data = mydata,
    family = binomial()
)

The easiest way to diagnose which structure you need is to look at your grouping variable labels. If plot “1” appears in multiple fields and refers to different things, you need nesting. If your observers have globally unique IDs regardless of which patient they see, use crossed effects.

How many groups do I need for the random effect to be reliable?

We wish we knew! There is no universally agreed minimum, and the question is better framed as “how many groups do I need to estimate the random effect variance with useful precision?” The answer is more than most applied researchers expect 😉.

As a rough guide, five groups is generally considered an absolute floor below which a random effect estimate is essentially uninterpretable. Below about 10–15 groups, estimates of the random effect variance can be highly variable, and singular fits are common. For reliable inference on variance components, 20–30 groups is a more comfortable starting point. The number of observations per group also matters since groups with very few observations contribute little information to the variance estimate and can destabilize estimation.

For our farm dataset, 25 farms is borderline comfortable. With a group-level predictor like irrigation, more farms would give meaningfully tighter estimates — as we saw in the large-sample demonstration above. If you are designing a study from scratch, this is an argument for prioritizing more groups over more observations within groups, particularly when the key inferential targets are group-level predictors.

My model returned a singular fit warning. What does that mean?

A singular fit occurs when the optimizer drives one or more variance components to (or very near) zero, or when a correlation parameter hits its boundary of \(\pm 1\). In lme4 this produces the warning boundary (singular) fit: see help('isSingular'), and in glmmTMB the equivalent warning references a Hessian⁷ that is not positive definite. The most common cause in practice is a random effects structure that is too complex for the data to support — either too few groups, too few observations per group, or a genuine absence of group-level variation.

The chunk below fits our model on a much-reduced subset of the data — six farms with only eight plants each — which gives the optimizer far too little information to estimate the farm variance reliably.

# take a small subset to deliberately provoke a singular fit
set.seed(42)
plants_small <- do.call(rbind, lapply(levels(plants$farm)[1:6], function(f) {
    sub <- plants[plants$farm == f, ]
    sub[sample(nrow(sub), 8), ]
}))
plants_small$farm <- droplevels(plants_small$farm)

# suppress the warning here just so it doesn't clutter the rendered output —
# in your own work, let these warnings surface and take them seriously
fit_singular <- suppressWarnings(
    lme4::glmer(
        survival ~ genotype + irrigation + (1 | farm),
        data = plants_small,
        family = binomial(link = "logit")
    )
)

boundary (singular) fit: see help('isSingular')

lme4::isSingular(fit_singular)

[1] TRUE

# let's do what the nice error message says!
# help('isSingular')
# if you do this on your own computer (assuming you've run `library(lme4)`),
# you'll get a help page pop up that has some very useful info!
lme4::VarCorr(fit_singular)

 Groups Name        Std.Dev.
 farm   (Intercept) 0       

When isSingular() returns TRUE and VarCorr() shows a variance estimate at or very near zero, the farm random effect has effectively been dropped from the model. In this situation the standard errors of the fixed effects will be underestimated, because the model is no longer accounting for the within-farm correlation of observations. The appropriate response is usually to simplify the random effects structure, gather more data, or reconsider whether the grouping factor genuinely introduces the kind of correlation you think it does.

My model failed to converge. What should I do?

Welcome to the club! Convergence failures are common and almost always trace back to one of a few sources. In rough order of how often they occur:

1. Is your random effects structure too complex for your data? Probably. Try simplifying by removing random slopes, reducing the number of grouping factors, or checking whether any grouping factor has too few levels. lme4::isSingular() is your first diagnostic — see the singular fit FAQ above.

2. Are your predictors on very different scales? This can cause numerical problems during optimization. Centering and scaling your continuous predictors before fitting rarely changes the substantive interpretation of coefficients (they can be back-transformed) but frequently stabilizes the optimizer.

3. Is the optimizer stuck in a local minimum? Both glmmTMB and lme4 support alternative optimizers. For lme4, glmerControl(optimizer = "bobyqa") is often the first thing to try when glmer fails. For glmmTMB, the control argument accepts a glmmTMBControl() object where you can switch from the default nlminb to alternatives. You should read up on these different optimizers before placing blind faith in them, as it will (a) make sure you know what you’re asking the model to do, and/or (b) may give you a needed reminder that what’s going on under the hoods of these models is actually super complicated and if they’re tools you’re planning on using, you should invest more time than it will take to read this guide to actually figure out what the hell the computer is doing 😎.

4. Do you have collinearity between predictors? This can inflate standard errors and can destabilize the Hessian matrix. Check variance inflation factors with car::vif() before fitting.

If none of these help, it may be that the model you have specified is genuinely not identified by the data you have. The appropriate response in that case is to reconsider the model, not to try more optimizers.

Are there Bayesian or other newer estimation approaches I should know about?

Yes! And they’re very good! But also, you should know what you’re doing to use them! For certain problems they are meaningfully better than the frequentist approaches covered in this guide.

brms is probably the most accessible entry point into fully Bayesian GLMMs in R. It uses Stan under the hood and provides a formula interface essentially identical to glmmTMB and lme4, making the transition from a frequentist workflow relatively painless. The Bayesian framework handles small group sizes and complex random effects structures more gracefully, because prior distributions regularize the estimates rather than allowing them to degenerate to the boundary. The cost is that prior specification requires thought, and chain diagnostics add another layer of complexity.

# a bayesian glmm with brms — the formula interface is almost identical
# i recommend weakly informative priors rather than flat ones
fit_brms <- brms::brm(
    survival ~ genotype + irrigation + (1 | farm),
    data = plants,
    family = bernoulli(link = "logit"),
    chains = 4,
    iter = 2000
)

INLA (Integrated Nested Laplace Approximations, implemented in the R-INLA package) offers fast approximate Bayesian inference for latent Gaussian models, which includes the GLMM as a special case. See a forthcoming guide on this tool for spatial usage! INLA is particularly popular in spatial statistics and epidemiology, where models with structured spatial random effects would be impractical to fit with MCMC. The learning curve is steeper than brms, but for large datasets with complex dependence structures it can be orders of magnitude faster.

Variational inference and Monte Carlo EM are two further alternatives on the frequentist side, the latter discussed in the appendix on estimation, though this guide’s author is not overly familiar with them. Both are active areas of research, and implementations are available in several R packages for those who want to explore them. The glmmTMB documentation and the lme4 vignettes are good starting points for understanding when these approaches might be preferable to the Laplace approximation.

When would I ever use PQL?

Penalized quasi-likelihood (PQL), introduced by N. E. Breslow and Clayton (1993), sidesteps the intractable likelihood integral entirely by working with a quasi-likelihood instead — a function that behaves like a likelihood in the sense that it can be maximised to obtain estimates, but does not require specifying a full probability model for the data. The key idea is to linearize the link function around current estimates of the random effects using a first-order Taylor expansion, which transforms the GLMM into an approximate linear mixed model that can be fit with standard REML. This approximate linear mixed model is then fit, the random effects are updated, and the process repeats until convergence.

The appeal of PQL is its simplicity and speed. Because it reduces the problem to a sequence of linear mixed model fits, it is computationally cheap and easy to implement. However, this convenience comes at a cost. The linearization is only an approximation, and the approximation is poor when data are sparse (few observations per group) or when the response distribution is far from Gaussian — exactly the conditions that tend to arise in practice with Bernoulli or Poisson data. PQL estimates of both fixed effects and variance components are known to be biased in these settings (Norman E. Breslow and Lin 1995; Bolker 2015), and the quasi-likelihood it maximises is not a true likelihood, meaning that likelihood-ratio tests and AIC are not directly applicable.

In short, with modern computing there is almost never a good reason to reach for PQL. The Laplace approximation is available in both glmmTMB and lme4 at negligible additional cost, and is substantially more accurate for binary data. PQL is available via MASS::glmmPQL if you need it, but its use with Bernoulli or Poisson responses should be approached with real caution.

What is emmeans and when should I use it?

The emmeans package (estimated marginal means, sometimes called least-squares means) is a post-hoc tool for computing model-based predictions and contrasts at specified covariate values, with proper propagation of uncertainty. For GLMMs it is particularly useful in two situations:

First, when you want predicted probabilities (or counts, or means) averaged over the observed distribution of covariates rather than held at a single reference value. predict() with re.form = NA sets the random effect to zero on the link scale, which is not the same as averaging over the random effect distribution — the distinction we worked through above. emmeans can perform the integration over the random effect distribution more rigorously, though the two approaches agree closely when the random effect variance is small. Jiye’s logit/probit guide covers the same territory for models without random effects, where it appears as the gap between odds ratios and average marginal effects, and is worth reading alongside this.

Second, when you want pairwise contrasts between factor levels with multiplicity adjustment. Rather than computing contrasts by hand from the coefficient table, emmeans handles this cleanly and outputs results on the response scale if requested.

library(emmeans)
# compute predicted survival probabilities for each genotype,
# averaged over irrigation and marginalizing over farms
emm <- emmeans::emmeans(fit_tmb, ~ genotype, type = "response")
pairs(emm)

For a GLMM with a non-identity link, always specify type = "response" if you want the output on the probability (or count, etc.) scale rather than the link scale.

An opinionated aside on the “labeling epidemic” in discipline-specific methodology

In many applied disciplines, there is a tendency to label classes of statistical models with names that communicate a cluster of information, such as the underlying probability model, the estimation procedure, and sometimes even the software likely used, all at once. These paraphyletic taxonomies are not so much well-defined classes of models as overlapping sets that are used to communicate within disciplinary boundaries. The result is that two researchers from different fields may use the same term to mean different things, or different terms to mean the same thing.

To be clear, we are not arguing that the term “GLMM” falls into this category. In fact, since the term “GLMM” was developed by mathematicians arising from a coherent theoretical framework, it refers unambiguously to a specific class of probability models characterized by a non-Gaussian response distribution, a link function, and at least one random effect — and writing one down in full generality is tractable, as we have done throughout this guide. The label and the math refer to the same thing, which is exactly the property the discipline-specific naming conventions so often lack.

Ultimately here we are not making the point that discipline-specific names for things are wrong out-and-out, but more that one needs to be careful when discussing these terms. For example, we have used the term “fixed effect” throughout to mean a particular thing. There are a whole class of models called “fixed effects” models that are mean something slightly different.

This is ultimately why we recommend, wherever possible, writing down your model explicitly. A probability model written in full is unambiguous in a way that a name never quite is. If you can write down the likelihood, the link function, and the random effects structure, then you and your collaborators are guaranteed to be talking about the same thing — regardless of what you call it. This habit also has a practical payoff in that it forces you to be precise about what you are actually assuming, which makes it much easier to catch misspecifications early and to explain your choices to reviewers or colleagues from other fields. Throughout this guide we have tried to model this practice, and we encourage you to carry it forward into your own work.

Estimation Procedures in Detail

This section works through the three main estimation strategies for GLMMs, and then takes a closer look at the object all of them are trying to optimize. The material here is not necessary for fitting a model, but it explains where the singular fits, convergence failures, and Hessian warnings that show up in practice actually come from.

Matrix Notation for Linear Models

First of all, a quick review on when the same things look different.

In a standard linear model, we observe \(n\) response values and \(p\) predictors for each observation. Rather than writing out a separate equation for each observation, we can stack everything into a single compact expression. Each individual observation follows:

\[\begin{equation} y_i = \sum_{j=1}^p x_{ij}\beta_j + \epsilon_i \end{equation}\]

but writing this out for every \(i = 1, \ldots, n\) is unwieldy. Instead, we collect all response values into a vector \(\mathbf{y}\), all predictor values into a matrix \(\mathbf{X}\), and all errors into a vector \(\boldsymbol{\epsilon}\), allowing the entire system to be expressed as:

\[\begin{equation} \mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\epsilon} \end{equation}\]

This is purely a notational convenience — it says the same thing as the sum above, just for all observations at once. The errors \(\boldsymbol{\epsilon}\) are assumed to be independent and identically distributed Gaussian random variables with mean 0 and variance \(\sigma^2\).

It’s not always clear why we choose one naming convention over another, so here are some general rules:

• Lowercase scalar (\(y_i\), \(x_{ij}\)) — a single number, either observed data or a specific element of a vector or matrix
• Lowercase bold (\(\mathbf{y}\), \(\mathbf{x}_i\), \(\boldsymbol{\epsilon}\)) — a vector of scalars; typically observed data or unknown errors
• Uppercase bold (\(\mathbf{X}\)) — a matrix; capitalised by linear algebra convention, not because it is a random variable
• Uppercase non-bold (\(Y_i\)) — a random variable; used in the probability model form to signal that we are reasoning about a distribution rather than a fixed observed value
• Uppercase bold (\(\mathbf{Y}\)) — a random vector; the full response vector treated as a collection of random variables, as in the probability model \(\mathbf{Y}|\mathbf{X} \sim \mathcal{N}(\mathbf{X}\boldsymbol{\beta}, \sigma^2\mathbf{I})\)

The key thing to keep in mind is that \(\mathbf{X}\) and \(\mathbf{Y}\) can look superficially similar but are capitalised for entirely different reasons — \(\mathbf{X}\) is capitalised because it is a matrix, \(\mathbf{Y}\) is capitalised because it is a random variable.

Equivalent Commands at a Glance

R and Julia are the two languages worked through in this guide. Python and Stata users will find the mapping below useful for getting oriented, though it is a translation table rather than a promise of equivalence — see the caveats underneath, which are not incidental.

^† These fit by variational Bayes or a Laplace approximation to the posterior rather than by maximum likelihood, and support independent random effects only. They are not drop-in equivalents to glmer and the numbers will not match. Pymer4 provides a genuine equivalent, but it does so by calling lme4 through rpy2, so it requires an R installation.

^‡ Stata’s gsem can be coaxed into some zero-inflated mixed specifications through latent class syntax, but not through a ziformula-style interface.

The pattern in this table is worth naming rather than leaving for the reader to infer. Basic GLMM fitting is well served everywhere. Once you need a negative binomial with random effects, zero inflation, or simulation-based residual diagnostics, glmmTMB and the R ecosystem around it are close to the only maximum-likelihood option available. That is a genuine reason to work in R for this class of model, and it is worth knowing before you commit to a language for a project.

Simulation code

The code below defines the function used to generate plants.csv. Copy it into an R script, then call

plants <- simulate_plants()

to reproduce the dataset. Because the function calls set.seed(123) internally, this will regenerate the file byte-for-byte on any machine running the same R version.

#' Simulate plant survival data with farm-level clustering
#'
#' Generates a synthetic plant-level dataset for demonstrating generalized
#' linear mixed models (GLMMs) with a binary outcome. Plants are nested
#' within farms; each farm is randomly assigned one irrigation method
#' (rainfed, overhead, or drip), and each plant is randomly assigned one of
#' four genotypes (G1-G4). Survival is generated from a logistic model with
#' fixed effects for irrigation and genotype, plus a farm-level random
#' intercept, so that farms sharing an irrigation method still differ from
#' one another in overall survival probability.
#'
#' @param seed Integer. Random seed passed
#'
#' @return A data frame with one row per plant and the following columns:
#' \describe{
#'   \item{farm}{Factor with 25 levels (\code{"farm_1"} through
#'     \code{"farm_25"}), identifying the farm a plant belongs to.}
#'   \item{irrigation}{Factor with levels \code{"rainfed"},
#'     \code{"overhead"}, \code{"drip"} (in that order, with
#'     \code{"rainfed"} as the reference level). Constant within farm.}
#'   \item{genotype}{Factor with levels \code{"G1"}, \code{"G2"},
#'     \code{"G3"}, \code{"G4"} (in that order, with \code{"G1"} as the
#'     reference level). Varies plant-by-plant within a farm.}
#'   \item{survival}{Integer, \code{0} or \code{1}, indicating whether the
#'     plant survived.}
#' }
#'
#' @examples
#' plants <- simulate_plants()
#' head(plants)
#'
#' @export
simulate_plants <- function(seed = 123) {
    set.seed(seed)

    # simulation parameters
    n_farms <- 25
    plants_per_farm <- sample(80:100, n_farms, replace = TRUE)

    irrigation_types <- c("rainfed", "overhead", "drip")
    farm_irrigation <- sample(irrigation_types, n_farms, replace = TRUE)

    intercept <- 0.5
    beta_overhead <- 0.20
    beta_drip <- 0.80
    beta_G2 <- 0.50
    beta_G3 <- -0.75
    beta_G4 <- 0.20

    farm_re_sd <- 0.8
    farm_re <- rnorm(n_farms, mean = 0, sd = farm_re_sd)

    plant_list <- lapply(seq_len(n_farms), function(i) {
        n <- plants_per_farm[i]
        genotype <- sample(c("G1", "G2", "G3", "G4"), n, replace = TRUE)

        irr_effect <- switch(farm_irrigation[i],
            overhead = beta_overhead,
            drip = beta_drip,
            0
        )

        geno_effect <- ifelse(genotype == "G2", beta_G2,
            ifelse(genotype == "G3", beta_G3,
                ifelse(genotype == "G4", beta_G4, 0)
            )
        )

        eta <- intercept + irr_effect + geno_effect + farm_re[i]
        p <- plogis(eta)
        survival <- rbinom(n, size = 1, prob = p)

        data.frame(
            farm = paste0("farm_", i),
            irrigation = farm_irrigation[i],
            genotype = genotype,
            survival = survival
        )
    })

    plants <- do.call(rbind, plant_list)

    plants$genotype <- factor(plants$genotype,
        levels = c("G1", "G2", "G3", "G4")
    )
    plants$irrigation <- factor(plants$irrigation,
        levels = c("rainfed", "overhead", "drip")
    )
    plants$farm <- factor(plants$farm)

    plants
}