Generalized Linear Mixed Model R

Getting GLMMs Running in R Without Losing Your Mind

Most people pick up generalized linear mixed models because they have clustered data and their data analyst says "just fit a mixed model." The packages available are glmmTMB, lme4, and brms. lme4 is the default answer. brms is for people who want Bayesian posteriors. glmmTMB is for people who actually need zero-inflated distributions or complicated random effect structures. I use glmmTMB most of the time because it handles things lme4 chokes on. A GLMM extends a generalized linear model by adding random effects to the linear predictor. The fixed effects tell you the population-level relationship. The random effects account for grouping structure in your data. If you're modeling count data from multiple hospitals, the hospital-level intercepts are random. The coefficient for your main exposure variable stays fixed across all hospitals. The link function connects the linear predictor to the mean of your response. Binomial data uses a logit link by default. Poisson data uses a log link. Negative binomial data also uses a log link. You specify the family and the link when you call the model function. This is standard GLM stuff plus the random effects layer on top.

The math involves integrating out the random effects. For continuous responses this is often doable with quadrature. For binary and count responses, the integral has no closed form. glmmTMB uses adaptive Gauss-Hermite quadrature with up to 15 points by default. More points are more accurate but slower. Fewer points can give biased estimates on messy data. I usually run with the default and check convergence diagnostics rather than cranking the integration points up unless something looks wrong.

Installation and Setup

You need the glmmTMB package. It pulls in TMB as a dependency, which is written in C++ and handles the automatic differentiation internally. The installation itself is straightforward on most systems but occasionally fails on Windows when the Rtools path is wrong. If you hit a compilation error, check that Rtools 4.x is installed and on your PATH. On macOS, you typically don't need anything extra. On Linux, make sure gcc and g++ are present. Install it with install.packages("glmmTMB"). Load it after. You also need tibble for clean data handling and dplyr for preprocessing. Here is a minimal setup block: library(glmmTMB)
library(tibble)
library(dplyr)

I keep those three loaded for nearly every analysis. The rest depends on what you're doing.

Fitting Your First Model

Let me show you a real example from a dataset I worked with last year. I had binary outcome measurements repeated on patients across clinics. The outcome was whether patients completed a follow-up visit within 90 days. There were about 4,000 patients spread across 60 clinics. The fixed effects included age, treatment group, and baseline severity score. The random effect was a clinic-level intercept. The formula looked like this: fit <- glmmTMB(followup ~ age + treatment + severity + (1 | clinic),
family = binomial,
data = mydata)

That is the basic syntax. The tilde separates the response specification from the predictors. The fixed effects sit on the left side of the formula. The random effects go inside the parentheses on the right. The pipe character before clinic tells glmmTMB that clinic is a grouping factor for the random intercept. You could also specify random slopes like (age | clinic) but that requires more clusters to estimate properly. Sixty clusters can handle a random intercept and slope, but it gets unstable closer to twenty. After fitting, I ran summary(fit). The output gives you the fixed effect coefficients with standard errors and z-values. It also gives you the variance components for the random effects. For the binomial family, the dispersion parameter is fixed at one. That is normal. If you need to relax that assumption, you use a different family or add an observation-level random effect, which I will get to shortly.

A Problem That Nearly Drove Me Crazy

Here is a specific edge case I encountered that is worth remembering. I was fitting a Poisson GLMM for readmission counts with a zero-inflated structure. The model converged but the summary reported a singular fit warning. The random effect variance for one of the grouping factors was effectively zero. The optimizer had pushed it to the boundary. I spent about four hours checking codes, re-specifying the model, and running diagnostics before I realized the grouping factor had too many levels with only a handful of observations per level. Forty-seven groups with a median of three observations each. A random effect with that much sparsity cannot be estimated reliably. The workaround was to aggregate the data to the group level and fit a zero-inflated negative binomial model with glmmTMB using the aggregated counts as the response. I dropped the random effect entirely and added an observation-level random effect to account for overdispersion. The results were stable and the inference was practically identical to what I would have gotten from the original model if it had converged properly. Aggregating before modeling is not always the right move, but in cases like this it saved the analysis. Another issue that trips people up is complete separation in binary GLMMs. If a predictor perfectly predicts the outcome within a cluster, the fixed effect coefficient runs to infinity. lme4 throws a warning and stops. glmmTMB may attempt to fit it and produce huge coefficients with enormous standard errors. I check for separation before fitting by running a table of the outcome crossed with the suspicious predictor inside each cluster. If any cell is empty in a way that matters, I drop the predictor or add a weakly informative prior using the glmmTMB option for that.

Zero-Inflation and Overdispersion

Zero-inflation is common in ecological count data and in clinical data where a large proportion of patients never experience the event of interest. glmmTMB handles zero-inflation natively with the zi ~ formula argument. You specify a separate submodel for the zero-inflation process. The default is a logistic submodel for the probability of being in the structural zero component. A typical zero-inflated Poisson model looks like this: fit_zi <- glmmTMB(counts ~ exposure + (1 | site),
ziformula = ~ 1,
family = poisson,
data = ecodata)

The ziformula argument controls the zero-inflation part. ziformula = ~ 1 means the zero-inflation probability is constant across all observations. You can add predictors there too, like ziformula = ~ habitat_type, if you suspect certain sites produce more zeros for reasons unrelated to the main count process. Overdispersion in count data is another frequent problem. A Poisson model assumes the mean equals the variance. Real data rarely satisfies that. The negative binomial family in glmmTMB adds a dispersion parameter. Use family = nbinom1 if the variance scales linearly with the mean. Use family = nbinom2 if the variance scales quadratically. nbinom2 is more common in practice. I check for overdispersion by comparing the residual deviance to the residual degrees of freedom. If the ratio is above two, I switch to negative binomial.

Model Diagnostics and Validation

glmmTMB does not have a full suite of diagnostic plots like lme4. I use DHARMa for residual diagnostics. It simulates residuals from the fitted model and checks for uniformity, overdispersion, and zero-inflation. The workflow is: library(DHARMa)
sim_residuals <- simulateResiduals(fit, nsim = 1000)
plot(sim_residuals) If the simulated residuals show a pattern, the model is misspecified. A U-shaped quantile residual plot indicates underdispersion or a missing term. A reverse U-shape suggests overdispersion. Heavy tails mean the distribution family is wrong. I check all of these before reporting results.

For convergence checks, look at the gradient norm in the summary output. Values below 0.01 are generally fine. Values above 0.1 suggest the optimizer did not find a proper optimum. If the gradient is large, try changing the optimizer. glmmTMB defaults to NLOPT_LN_BOBYQA. I sometimes switch to NLOPT_LN_PRAXIS by adding control = glmmTMBcontrol(optctrl = list(method = "NLOPT_LN_PRAXIS")) if BOBYQA struggles. It is slower but more robust on tricky likelihood surfaces.

Comparing Models and Selecting Structures

People ask about AIC for model selection. You can use AIC(fit) to get the value. Lower is better. Likelihood ratio tests work for nested models when the null hypothesis is not on a boundary. They do not work reliably for comparing zero-inflated versus non-zero-inflated models because the null places the zero-inflation parameter at a boundary. Use AIC or cross-validation instead for those comparisons. I run a stepwise procedure by hand rather than with an automated function. I start with the maximal random effects structure justified by the design. Then I remove random slopes one at a time if the data do not support them. I keep fixed effects that are theoretically meaningful regardless of p-values. Random effect removal is data-driven. Fixed effect removal should be substantive unless the model is clearly overparameterized.

Common Pitfalls and Limitations

glmmTMB is fast but not magic. It struggles with very sparse binary data. If most clusters have zero events or all events, the model cannot estimate the fixed effects properly. Bayesian approaches with brms handle this better because the prior regularizes the estimates. If you have fewer than ten clusters with binary outcomes, consider brms with weakly informative priors instead of fighting with glmmTMB. Another limitation is that glmmTMB does not support complex covariance structures for multivariate random effects as flexibly as nlme. If you need a full unstructured covariance matrix for multiple random slopes across many groups, glmmTMB will attempt it but convergence becomes unreliable past three or four random effects per group. Stick to diagonal or simple structured covariance matrices in those cases. Parallel processing is supported for refitting with bootMer and for some optimization routines, but it is not automatic. You set it up manually with the future package and pass it through the control arguments. This can cut computation time roughly in half for large datasets with many clusters, depending on your machine.

Exporting Results for Reporting

There are several packages for formatting glmmTMB output into tables suitable for manuscripts. emmeans gives you marginal means and pairwise comparisons with correct standard errors that account for the random effects. You call emmeans(fit, ~ treatment) after fitting. It back-transforms from the link scale automatically. For table formatting, the parameters package or the broom.mixed extension works. broom.mixed::tidy(fit) pulls the fixed effects into a clean data frame. broom.mixed::glance(fit) gives you model-level statistics like AIC, BIC, and log-likelihood. I usually combine broom.mixed for the coefficient table and emmeans for the estimated marginal means. That covers almost everything a results section needs. I avoid sjPlot because it introduces dependencies that are harder to manage in reproducible workflows and its theming options are limited.

When to Reach for Something Else

If your data have hierarchical structure at three or more levels with many observations per group and you want a full Bayesian posterior, brms is the better tool. It is slower but handles priors, missing data, and complex hierarchical structures more gracefully. If you have thousands of random effect levels and need speed, glmmTMB is still faster than brms for most standard models. If you need penalized quasi-likelihood for overdispersed binary data and cannot fit a GLMM, consider the glmmPQL function from the MASS package as a fallback, though it is approximate and less reliable than full likelihood methods. For most practical GLMM work in R, glmmTMB covers the range from simple random intercepts to zero-inflated negative binomial models with crossed random effects. The learning curve is moderate. The documentation is adequate but sparse on the details of what happens during optimization. Reading the TMB manual helps if you want to understand what the package is actually doing under the hood. That level of understanding matters when the model refuses to converge and you need to decide whether to change the data, change the model, or change the optimizer. The code examples above assume you have a tidy data frame with your response variable and all predictors ready. Factor variables need to be coded as factors before fitting, or glmmTMB will treat them as numeric and the random effects grouping will be wrong. I spend more time on data preparation than on model fitting. That is normal. The model fitting itself usually takes seconds to minutes. The data wrangling takes hours.