Fitting Logistic Mixed Models with glmer

Why fitting a mixed logistic model requires care

When a binary clinical outcome is recorded for patients treated at several clinics, patients from the same clinic may have related outcomes. A logistic mixed model represents this clustering by allowing each clinic to have its own random intercept. The fitting task is not only to call software: the response must be encoded correctly, the formula must match the study design, and the fitted object must be checked before its estimates are used.

This tutorial fits one binary-outcome, random-intercept model with lme4 in R. It assumes the outcome and cluster identifiers have already been audited and the data are in one row per patient. The earlier material on translating a study design into R formula syntax supplies the basic formula conventions used here.

Model: For patient $i$ in clinic $j$, let $Y_{ij}$ equal 1 if the specified event occurs and 0 otherwise. A random-intercept logistic model is $$ Y_{ij}\mid b_j \sim \operatorname{Bernoulli}(p_{ij}), \qquad \operatorname{logit}(p_{ij}) = \beta_0 + \boldsymbol{x}_{ij}^{\mathsf T}\boldsymbol{\beta} + b_j, \qquad b_j \sim N(0,\sigma_b^2). $$ Here $b_j$ is the clinic-specific intercept deviation, and $\sigma_b^2$ is the between-clinic random-intercept variance.

The model assumes that, conditional on the clinic random effect and included predictors, patient outcomes are independent Bernoulli observations with probabilities described by the specified logit mean. Clinic random effects are assumed independent across clinics and normally distributed with mean zero and common variance; the usual random-effects formulation also assumes they are independent of the included predictors. The model must be scientifically reasonable, and the data must contain enough information across clinics and outcome levels to estimate its parameters. These are assumptions to assess, not facts guaranteed by fitting software.

From the model to the glmer call

For a given clinic effect $b_j$, the conditional likelihood contribution from clinic $j$ is the product of that clinic’s Bernoulli probabilities. The likelihood for the observed data must account for the fact that $b_j$ is unobserved. It does so by integrating the conditional likelihood over the normal distribution of $b_j$:

L(\boldsymbol{\beta},\sigma_b^2) = \prod_{j=1}^{J} \int \left\{ \prod_{i=1}^{n_j} p_{ij}(b_j)^{y_{ij}} [1-p_{ij}(b_j)]^{1-y_{ij}} \right\} \phi(b_j;0,\sigma_b^2)\,db_j.

The integration explains why this is not simply an ordinary logistic regression with an extra predictor: the clinic effects are latent, and the likelihood averages over their assumed distribution. glmer fits generalized linear mixed models by approximating this likelihood calculation and optimizing the resulting objective. With its default settings, it uses the Laplace approximation. This is a numerical fitting method; it does not remove the need to check whether the model is supported by the data.

In an R formula, fixed effects appear as ordinary predictors, while (1 | clinic_id) requests a random intercept for each level of the clinic identifier. The binomial family is used for the 0/1 response; in this setting its default link is the logit.

Worked Example: Fit a Clinical Random-Intercept Model

Suppose a hypothetical clinical dataset contains 240 patients from 12 clinics. The outcome event12 records whether a patient experienced a specified clinical event within 12 months, with values 1 for an event and 0 otherwise. Treatment is either Control or Drug, and age is recorded in years. The data include 120 patients in each treatment group and 78 events in total. These totals describe the hypothetical data; they do not determine the fitted coefficients by themselves.

First ensure that the outcome is genuinely coded as 0 and 1. If it is already numeric, do not convert it needlessly. If it is a factor or character vector whose labels are exactly “0” and “1”, convert the labels—not the factor’s internal level codes:

if (is.factor(clinical$event12) || is.character(clinical$event12)) { stopifnot(all(na.omit(unique(as.character(clinical$event12))) %in% c("0", "1"))); clinical$event12 = as.integer(as.character(clinical$event12)) }

Then check the resulting values and prepare the predictors and grouping variable. The age variable age10 is centered at 60 years and measured in decades: age10 equals (age in years minus 60) divided by 10.

stopifnot(is.numeric(clinical$event12), all(na.omit(clinical$event12) %in% c(0, 1)))

clinical$treatment = relevel(factor(clinical$treatment), ref = "Control"); clinical$clinic_id = factor(clinical$clinic_id); clinical$age10 = (clinical$age_years - 60) / 10

Now fit the model. The Control group is the treatment reference category, age10 is a fixed-effect predictor, and each clinic receives a random intercept.

fit = lme4::glmer(event12 ~ treatment + age10 + (1 | clinic_id), data = clinical, family = binomial)

For this hypothetical fit, suppose the relevant output is an intercept estimate of -1.200, a Drug coefficient of -0.550, an age10 coefficient of 0.180, and a clinic random-intercept standard deviation of 0.450. These are illustrative fitted values, not estimates calculated from the three totals above. The corresponding random-intercept variance is $0.450^2 = 0.2025$.

The intercept is the fitted log odds for a Control patient aged 60 in a clinic whose random intercept is zero. For that reference profile, the fitted probability is $1/(1+\exp(1.200)) \approx 0.2315$. For a Drug patient of the same age in that same zero-random-effect clinic, the linear predictor is $-1.200-0.550=-1.750$, giving probability $1/(1+\exp(1.750)) \approx 0.1480$. These calculations illustrate how the fitted model produces probabilities conditional on a specified random effect; they are not a full interpretation of treatment effects.

Basic quantities can be extracted without relying on manually transcribing a printed table: fixef(fit) returns fixed-effect estimates, as.data.frame(VarCorr(fit)) displays the random-effect variance information, and isSingular(fit, tol = 1e-4) checks for a singular fit. Also verify the number of observations used with nobs(fit). The broader interpretation of the summary output is a separate task.

Response coding and the fitted object

The response coding check is essential because R factors store integer level codes internally. For example, if a factor has levels “0” and “1”, directly applying as.integer() can return 1 and 2, not the intended 0 and 1. A binomial fit may then fail or represent the outcome incorrectly. Converting through the labels, as in the example, preserves the intended coding when the labels truly are “0” and “1”. If labels instead say “No” and “Yes”, explicitly map those labels to 0 and 1 after checking which label represents an event.

The fitted object, here called fit, stores the model call, estimates, random-effect information, and fitting diagnostics. Retain it as a named object rather than fitting the model only to print its output. This makes it possible to check the model specification, extract results reproducibly, and compare the fitted object with the analysis dataset.

At a minimum, check that the formula and data in the call are the ones intended, the response has only the valid values, the observed number of rows used is expected, and the fit did not produce convergence warnings. Inspect the random-intercept variance and the singular-fit check as well. A small estimated variance or singular-fit warning is not a reason to silently delete the random effect: it is evidence that the data may not support the requested random-effects structure, and the model specification needs review.

Worked Example: Verify Factor-to-Binary Recoding

Consider a short hypothetical response vector stored as a factor with labels “0” and “1”: factor(c("0", "1", "0", "1"), levels = c("0", "1")). The intended event indicators are 0, 1, 0, 1.

Applying as.integer() directly returns 1, 2, 1, 2, because those are the internal codes for the first and second factor levels. That result does not satisfy the required 0/1 coding. Instead, first recover the labels and then convert them: as.integer(as.character(x)) returns 0, 1, 0, 1. A subsequent check such as all(x %in% c(0, 1)) confirms that every nonmissing value is valid.

The label check matters: this conversion is appropriate only when the labels themselves are “0” and “1”. If the factor labels are “No” and “Yes”, converting them this way will not work. After verifying the labels, use an explicit mapping such as ifelse(x == "Yes", 1, 0); first confirm that the only observed nonmissing labels are the two expected categories and that “Yes” denotes the event. This prevents a coding shortcut from changing the meaning of the response.

Checks, failure modes, and alternatives

A successful function call is not by itself evidence of an adequate fit. Check for warnings about convergence and inspect whether the variance estimate is on the boundary, including a singular-fit result. Also check whether the intended clinics and patients are represented, whether missing data unexpectedly removed observations, and whether the binary outcome has enough variation within and across clinics. Very sparse event patterns or near-separation can make estimates unstable. Diagnostics should be interpreted alongside the design and the outcome pattern, not used as automatic pass/fail labels.

If the optimizer reports a convergence warning, first verify the data, factor levels, model formula, and predictor scales. Rescaling a predictor can sometimes improve numerical behavior. Changing an optimizer or its settings may be considered after these checks, but it is not a cure for a misspecified model or insufficient information; record the change and confirm that the resulting estimates are stable. Do not suppress a warning and report the fit as problem-free.

A standard logistic regression fit with glm is simpler, but it does not model the clinic-level dependence represented by the random intercept. Treating clustered patient records as independent may lead to misleading uncertainty estimates. At the other extreme, representing every clinic with fixed indicator terms is a different modeling choice: it estimates a separate clinic effect rather than a distribution of clinic deviations and may be awkward when there are many clinics. Generalized estimating equations offer a population-averaged approach with a different target and assumptions. None of these alternatives is interchangeable with the random-intercept model solely because the outcome is binary.

The random-intercept model also has limits. It assumes a particular distribution for clinic effects and a common random-intercept variance, and it assumes the included predictors describe the mean structure adequately. With few clinics, the variance estimate can be imprecise. A model that converges may still be poorly chosen, and changing the random-effects structure changes what is being modeled. These issues should be addressed by study design, scientific reasoning, and diagnostic evidence.

Worked Example: Obtain a Clinic-Conditional Prediction

Suppose the fitted model from the clinical example estimates a random intercept of $0.300$ for one clinic. For a Drug patient aged 60, age10 is 0, so the linear predictor for this clinic is the intercept plus the treatment coefficient plus the estimated clinic effect:

\widehat{\eta} = -1.200 - 0.550 + 0.300 = -1.450.

The corresponding fitted probability is

\widehat{p} = \frac{\exp(-1.450)}{1+\exp(-1.450)} \approx 0.1900.

This is a point prediction conditional on the estimated random effect for that clinic and the stated patient covariates. It is not an average probability over clinics. In R, predictions that exclude random effects (for example, using predict(fit, re.form = NA, type = "response")) evaluate the fixed-effects part with random effects set to zero; that also should not be described as integrating over the clinic distribution to obtain a population-average probability.

A fitting workflow to retain

Check the analysis columns.
Confirm one row per intended observation, valid 0/1 response values, meaningful treatment reference levels, and a cluster identifier with the expected levels.

Write the model in design terms.
State the binary outcome, fixed predictors, and clustering structure before translating them into the formula.

Fit and retain the object.
Use glmer with the binomial family, the intended dataset, and a named fitted object.

Verify before interpreting.
Check the observations used, warnings, fixed-effect estimates, random-intercept variance, and singular-fit status; then assess diagnostics and scientific plausibility.

Common pitfall: A factor labelled “0” and “1” is not safely converted to a binary response by applying as.integer() directly. The factor’s internal codes are 1 and 2. Convert verified labels with as.integer(as.character(x)), or explicitly map verified categories, and check that the resulting nonmissing values are exactly 0 and 1.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions below.

In the formula event12 ~ treatment + (1 | clinic_id), which term requests a clinic-specific random intercept?
A factor response labelled “0” and “1” is converted with as.integer(). What values can result, and how should the conversion be checked?
What does a singular-fit check assess, and why should a singular result not be ignored automatically?
A prediction sets the random effect to zero. Is it generally the same as a population-average probability across clinics?

Solutions. (1) (1 | clinic_id) requests a random intercept for each clinic. (2) Direct conversion can produce 1 and 2, the internal factor codes; convert verified labels through character values, then check that all nonmissing values are in {0, 1}. (3) It flags a random-effects fit whose estimated structure is on or near a boundary, such as a variance effectively zero; this may indicate limited information or an unsupported structure and calls for review. (4) No. Setting the random effect to zero gives a prediction at that value of the random effect, whereas a population-average probability requires averaging over the random-effect distribution.