My Mixed Model Won’t Converge with Crossed Random Effects

Symptoms

• A crossed subject-by-item model reports a gradient or Hessian warning.
• A slope variance is near zero, a correlation approaches ±1, or optimizer changes appear to solve the problem.

What this usually means

A convergence warning concerns numerical checks; a singular fit concerns a covariance matrix on the boundary. Neither tells the whole story. A random slope can also be unsupported by the design without triggering either message. In a crossed experiment, subjects and items are distinct grouping factors, not automatically nested.

Common causes

• Predictors have very different scales or fixed-effect columns are redundant.
• A predictor does not vary within the grouping factor assigned its random slope.
• There are too few groups or too little replication for the requested covariance structure.

Run these checks

1. Inspect the exact warning, fixed-effect rank, predictor scales, group counts, and replication.
2. Tabulate predictor variation within each grouping factor.
3. Inspect variance components and use isSingular(); absence of a flag is not proof of identifiability.
4. Refit with a different optimizer or stricter tolerances and compare coefficients and likelihoods. The lme4 documentation describes allFit() for broader comparison.
5. Document any scientifically justified simplification and its effect on conclusions.

What not to do

Do not increase iteration limits until a warning disappears and call that validation. Do not delete item effects just because subjects are the familiar unit, or choose the structure giving the preferred p-value.

Treatment options

Center or rescale predictors when warranted, correct coding errors, and compare optimizers. Retain slopes supported by within-group variation and the design. If a covariance structure is too rich, consider an explicitly justified simpler structure or regularization; explain the resulting assumptions.

Worked example

The synthetic data cross 30 subjects with 20 items, giving 600 observations. Treatment is between subjects: each subject has exactly one predictor value. The attempted subject-specific treatment slope cannot be separated from the subject intercept without additional constraints. For predictor values ±0.5, the data identify at most two group-specific intercept variances, while an intercept–slope covariance matrix asks for three parameters.

Here, isSingular() returned FALSE even for the unsupported structure. That is the teaching point: software diagnostics cannot replace a design check. Removing the unsupported subject slope while retaining both crossed intercepts gives a treatment coefficient of 0.7950. The default and bobyqa optimizers agree to the displayed precision, and the supported model reports no convergence messages.

This example diagnoses structural support rather than manufacturing a convergence failure. Your warning may instead be numerical or boundary-related; follow the checks before applying the same remedy.

What to tell the reviewers

We checked predictor variation within subjects and items. Because treatment was constant within each subject, we removed the unsupported subject-specific treatment slope while retaining crossed subject and item intercepts. Alternative optimizers agreed on the fitted treatment coefficient. We report the revised structure and numerical diagnostics.

See it in R and Python

Python implements maximum likelihood specifically for this complete balanced Gaussian crossed-intercept design using NumPy and SciPy. It checks the unsupported subject slope and compares two optimizers for the supported model. It is not a general mixed-model package and does not reproduce lme4’s isSingular() diagnostic.

Python dependencies: NumPy and SciPy. Install with python -m pip install numpy scipy.

Download the shared synthetic CSV beside the Python script before running it. Identical seeds in R and Python do not generally generate identical samples.

# Requires lme4; no automatic installation. Synthetic crossed design.
if (!requireNamespace('lme4', quietly=TRUE)) stop('Install lme4 before running this example.')
library(lme4)
set.seed(6006)
d <- expand.grid(subject=factor(1:30),item=factor(1:20))
# Between-subject predictor: no within-subject variation. A subject-specific
# slope for it is confounded with that subject's intercept.
d$x <- rep(rep(c(-.5,.5),each=15), times=20)
d$y <- .7*d$x + rnorm(30)[as.integer(d$subject)] + rnorm(20,sd=.5)[as.integer(d$item)] + rnorm(nrow(d))
within <- tapply(d$x,d$subject,function(z) length(unique(z)))
cat(sprintf('Rows=%d; subjects=%d; items=%d; max within-subject x values=%d\n',nrow(d),nlevels(d$subject),nlevels(d$item),max(within)))
bad <- lmer(y ~ x + (1+x|subject) + (1|item),d,REML=FALSE)
fit <- lmer(y ~ x + (1|subject) + (1|item),d,REML=FALSE)
other <- update(fit,control=lmerControl(optimizer='bobyqa',optCtrl=list(maxfun=100000)))
cat(sprintf('Unsupported subject slope singular=%s; supported model singular=%s\n',isSingular(bad),isSingular(fit)))
cat(sprintf('Supported x=%.4f; alternative optimizer x=%.4f; absolute difference=%.8f\n',fixef(fit)['x'],fixef(other)['x'],abs(fixef(fit)['x']-fixef(other)['x'])))
cat('Supported model convergence messages:\n'); print(fit@optinfo$conv$lme4$messages)
cat('lme4 version:',as.character(packageVersion('lme4')),'\n')
stopifnot(all(within==1),abs(fixef(fit)['x']-fixef(other)['x'])<1e-4)

# Maintainer-only export; run from repository root. Synthetic data only.
if (Sys.getenv("CLINIC_EXPORT_DATA") == "1") {
write.csv(d,'analysis-clinic/examples/case-006-shared-data.csv',row.names=FALSE,quote=FALSE)
}

# DASS Analysis Clinic Case 006; Python dependencies: numpy, scipy.
from pathlib import Path
import numpy as np
from scipy import stats, optimize
HERE = Path(__file__).resolve().parent

def ols(X, y):
    b = np.linalg.lstsq(X, y, rcond=None)[0]
    residual = y - X @ b
    df = len(y) - X.shape[1]
    cov = (residual @ residual / df) * np.linalg.inv(X.T @ X)
    se = np.sqrt(np.diag(cov))
    ci = np.column_stack((b-stats.t.ppf(.975, df)*se,b+stats.t.ppf(.975,df)*se))
    return b, cov, ci

def power(n, d):
    # Matches R power.t.test(strict=FALSE): rejection tail in effect direction.
    return stats.nct.sf(stats.t.ppf(.975,2*n-2),2*n-2,d*np.sqrt(n/2))

# Specialized ML for this complete balanced Gaussian crossed-intercept design.
# Not a general mixed-model package, and not lme4's convergence/singularity test.
d=np.genfromtxt(HERE/'case-006-shared-data.csv',delimiter=',',names=True)
y=d['y']; x=d['x']; subject=d['subject']; item=d['item']
subjects=np.unique(subject); items=np.unique(item); ns=len(subjects); ni=len(items)
assert len(y)==ns*ni and len(set(zip(subject,item)))==ns*ni
assert all(len(np.unique(x[subject==s]))==1 for s in subjects)
X=np.column_stack((np.ones(len(y)),x)); beta=np.linalg.lstsq(X,y,rcond=None)[0]
r=y-X@beta
grand=np.full(len(y),r.mean())
s=np.array([r[subject==g].mean() for g in subject])-grand
i=np.array([r[item==g].mean() for g in item])-grand
e=r-grand-s-i
ss=np.array([e@e,s@s,i@i,grand@grand]); dims=np.array([(ns-1)*(ni-1),ns-1,ni-1,1])
def objective(logvar):
 vs,vi,ve=np.exp(logvar)
 eigen=np.array([ve,ve+ni*vs,ve+ns*vi,ve+ni*vs+ns*vi])
 return .5*np.sum(dims*np.log(eigen)+ss/eigen)
a=optimize.minimize(objective,np.zeros(3),method='L-BFGS-B',bounds=[(-16,8)]*3)
b=optimize.minimize(objective,np.zeros(3),method='Powell',bounds=[(-16,8)]*3)
print('Within-subject x values: 1; unsupported subject slope is confounded with intercept.')
print('Supported treatment coefficient:',beta[1]); print('Variance components subject/item/residual:',np.exp(a.x))
print('Optimizer success:',a.success,b.success,'objective difference:',abs(a.fun-b.fun))
assert a.success and b.success and abs(a.fun-b.fun)<1e-4 and np.isclose(beta[1],.7950,atol=.00005)

Download Python script · Download shared synthetic CSV

Sources

• lme4 authors. Assessing convergence for fitted models.
• lme4 authors. Test fitted model for (near) singularity.
• lme4 authors. Refit a fitted model with available optimizers.