Answer

What does a singular fit mean in a mixed model?

Inspired by a question on Cross Validated · November 8, 2026

multilevel modelsmodel selection

What the warning means

A mixed model estimates a covariance matrix for its random effects. Variances cannot be negative and correlations must lie between -1 and 1, so the estimates live inside a bounded region. A fit is singular when the estimate lands on the edge of that region. In lme4 you see the message boundary (singular) fit and isSingular() returns TRUE.

• One random effect: its variance is estimated at exactly zero.
• Two correlated random effects: a variance is zero or the correlation is -1 or +1.
• Three or more: the matrix can be singular with no zero variance and no perfect correlation, because some combination of the effects has no variance. rePCA() in lme4 shows this.

This is different from a convergence failure. The optimizer has usually done its job, and the best value allowed really is on the boundary.

Why the estimate lands on zero

Take balanced data with k groups of n observations. The between-group mean square has an expected value of the residual variance plus n times the group variance. So the natural estimate of the group variance is (MS between - MS within) / n, and by chance MS between can come out smaller than MS within. That would make the estimate negative, which is impossible, so REML reports zero.

That is common whenever there are few groups or the group variance is small compared with the noise. With 8 groups of 5, a true group SD of 0.4 and a residual SD of 1, the exact chance of a zero estimate is 21.4%, even though the groups really do differ. Models that ask for several random slopes and correlations from a modest number of groups run into the same limit, only in more dimensions.

What it does to your results

• Zero is not proof of no variance. It means the data are compatible with little or none, but a confidence interval for the variance would usually reach well above zero.
• The term stops doing anything. When a single variance is estimated at zero, the model gives the same answer as the model without that term. In the example below the mixed model and plain regression agree to four decimals.
• Standard errors can be too small. If the true variance is above zero, the fit ignores the clustering. In the example the reported standard error is 0.1777, while the true variances imply 0.2121.
• Boundary tests need care. A chi-square test of whether a variance is zero is conservative, because zero is at the edge of the allowed range.

Can you keep only the random slope?

Sometimes a model with only random intercepts is singular, a slope-only structure such as (0 + x | id) fits better, and dropping the intercept looks like the obvious fix. Be careful, for three reasons.

The two models do not answer that question. They are not nested, and if they have the same number of parameters a likelihood ratio test has zero degrees of freedom and means nothing, regardless of any boundary. AIC can rank them, but it only says which of those two works better. A zero intercept variance in a model with no slope variance tells you little about the intercept variance in a model that has both, since unmodeled slope differences can distort it. Fit (1 + x | id) first.

A slope-only model depends on where x = 0 is. With correlated effects, the variance between people at a given x is the intercept variance + 2x times the covariance + x² times the slope variance. Removing the intercept variance and the covariance forces people's lines to meet exactly at x = 0. If that point is arbitrary, the assumption is too. The uncorrelated form (1 + x || id) also changes with centering, for a similar reason, though it is less restrictive.

It can be right when the design says so. If x is time since a baseline at which everyone starts from the same value, for example because scores were set to zero there, lines that fan out from x = 0 are a real claim you can defend. This decision differs from removing the intercept from a regression, which forces the fixed line through the origin. A slope-only random structure keeps the fixed intercept and removes only person-to-person variation around it.

What to do instead

1. Check the data first. Look at the number of groups, the observations per group, and a plot of each group's fitted line. Coding errors and a predictor that barely varies within groups also cause this.
2. Rule out the optimizer. allFit() refits with other optimizers. If they agree, the boundary is real; it will not remove a genuine singularity.
3. Keep the term if the design requires it. A random effect that reflects how the data were collected can stay in a confirmatory model (Barr et al., 2013). With one variance at zero, the fit is unchanged either way.
4. Or simplify for a stated reason. Try uncorrelated effects, then remove the terms the data cannot support (Bates et al., 2015; Matuschek et al., 2017). Say that you did this and why.
5. Or use a weak prior. blme or brms keep the variance away from exactly zero (Chung et al., 2013), which gives more honest uncertainty when there are few groups.

Not sure a mixed model is the right analysis? Try the test chooser, and the DASS blog has more on regression and model building. For the estimation side, see REML or ML for likelihood ratio tests.

See it in R and Python

Part 1 simulates the 8-groups-of-5 design. Part 2 simulates 30 people measured at x = 0 to 4 who differ only in slope, so their lines fan out from x = 0, and compares slope-only and full models using raw and centered x. R uses lme4. The Python tab writes out the REML and ML likelihoods with NumPy and SciPy, because statsmodels is not used here.

library(lme4)

# Part 1: 8 groups of 5, true group SD 0.4, residual SD 1
set.seed(8)
d <- data.frame(g = factor(rep(1:8, each = 5)))
d$y <- 10 + rep(rnorm(8, 0, 0.4), each = 5) + rnorm(40)

fit <- lmer(y ~ 1 + (1 | g), data = d)        # REML
print(VarCorr(fit))
print(isSingular(fit))

# The reason: between-group mean square is below the within-group one
print(round(anova(lm(y ~ g, data = d))[["Mean Sq"]], 3))

# With the variance at zero the fit is the ordinary regression
print(round(c(lmer_est = fixef(fit), lmer_se = sqrt(vcov(fit)[1, 1])), 4))
print(round(coef(summary(lm(y ~ 1, data = d)))[1:2], 4))

# SE of the mean if the true variances were known
print(round(sqrt(0.4^2 / 8 + 1 / 40), 4))

# Exact chance of a zero estimate in this design
print(round(pf(1 / (1 + 5 * 0.4^2), 7, 32), 3))

# Part 2: 30 people, 5 occasions; people differ only in slope, fanning out from x = 0
set.seed(2)
long <- expand.grid(x = 0:4, id = factor(1:30))
b <- rnorm(30, 0, 0.5)
long$y <- 5 + long$x + b[long$id] * long$x + rnorm(150)
long$xc <- long$x - 2                         # same data, x centred

ll <- function(f) round(as.numeric(logLik(lmer(f, data = long, REML = FALSE))), 2)
print(c(slope_only_raw     = ll(y ~ x  + (0 + x  | id)),
        slope_only_centred = ll(y ~ xc + (0 + xc | id)),
        full_raw           = ll(y ~ x  + (1 + x  | id)),
        full_centred       = ll(y ~ xc + (1 + xc | id))))

import numpy as np
from scipy import stats
from scipy.optimize import minimize


def deviance(theta, y, X, Z, g, reml):
    # V = I + Z(LL')Z' within groups; profiled over the residual variance
    L = np.array([[theta[0]]]) if Z.shape[1] == 1 else np.array([[theta[0], 0], [theta[1], theta[2]]])
    W = Z @ L
    n, p = X.shape
    V = np.eye(n) + (W @ W.T) * (g[:, None] == g[None, :])
    Vi = np.linalg.inv(V)
    XVX = X.T @ Vi @ X
    r = y - X @ np.linalg.solve(XVX, X.T @ Vi @ y)
    df = n - p if reml else n
    dev = df * (1 + np.log(2 * np.pi * (r @ Vi @ r) / df)) + np.linalg.slogdet(V)[1]
    return dev + (np.linalg.slogdet(XVX)[1] if reml else 0)


def fit(y, X, Z, g, reml=True):
    q = Z.shape[1]
    start, bnds = ([0.5], [(0, None)]) if q == 1 else ([0.5, 0, 0.5], [(0, None), (None, None), (0, None)])
    return minimize(deviance, start, args=(y, X, Z, g, reml), bounds=bnds, method="L-BFGS-B")


# Part 1: 8 groups of 5, true group SD 0.4, residual SD 1
rng = np.random.default_rng(5)
g = np.repeat(np.arange(8), 5)
y = 10 + rng.normal(0, 0.4, 8)[g] + rng.normal(size=40)
one = np.ones((40, 1))
print("group SD / sigma ratio:", round(float(fit(y, one, one, g).x[0]), 4))  # 0 means singular

means = np.array([y[g == j].mean() for j in range(8)])
msb = 5 * np.sum((means - y.mean()) ** 2) / 7
msw = sum(((y[g == j] - means[j]) ** 2).sum() for j in range(8)) / 32
print("MS between, within:", round(msb, 3), round(msw, 3))

# With the variance at zero the fit is the ordinary regression
print("mean, OLS SE:", round(y.mean(), 4), round(y.std(ddof=1) / np.sqrt(40), 4))
print("SE with known variances:", round(np.sqrt(0.4**2 / 8 + 1 / 40), 4))
print("P(zero estimate):", round(stats.f.cdf(1 / (1 + 5 * 0.4**2), 7, 32), 3))

# Part 2: 30 people, 5 occasions; people differ only in slope, fanning out from x = 0
rng = np.random.default_rng(2)
ids = np.repeat(np.arange(30), 5)
x = np.tile(np.arange(5.0), 30)
y = 5 + x + rng.normal(0, 0.5, 30)[ids] * x + rng.normal(size=150)

for label, c in (("raw", 0), ("centred", 2)):
    xs = x - c
    X = np.column_stack([np.ones(150), xs])
    slope_only = -fit(y, X, xs[:, None], ids, reml=False).fun / 2
    full = -fit(y, X, X, ids, reml=False).fun / 2
    print(label, "slope-only:", round(slope_only, 2), "full:", round(full, 2))

R output from one seeded run (same numbers on a second run):
boundary (singular) fit: see help('isSingular')
 Groups   Name        Std.Dev.
 g        (Intercept) 0.0000  
 Residual             1.1238  
[1] TRUE
[1] 0.989 1.323
lmer_est.(Intercept)              lmer_se 
              9.9118               0.1777 
[1] 9.9118 0.1777
[1] 0.2121
[1] 0.214
    slope_only_raw slope_only_centred           full_raw       full_centred 
           -256.17            -290.85            -253.62            -253.62

All figures come from one seeded run of the R code, and the seeds were picked so that part 1 lands on the boundary, which this design does about one time in five. Python's random numbers differ from R's, so its figures differ slightly, but the pattern is the same: a zero variance, MS between below MS within, and the same result for the full model under both codings.

• Part 1. The true group SD is 0.4, but the estimate is 0.0000, because MS between (0.989) is below MS within (1.323). The mean and its standard error match the plain regression (9.9118 and 0.1777).
• Part 2. The full model has a log-likelihood of -253.62 with either coding of x. The slope-only model gets -256.17 with raw x, where the lines really do meet, but -290.85 with centered x. Same data, same number of parameters: only the origin moved.

How to report a singular fit in APA style (7th edition)

In the method section, give the random-effects structure, the software, the estimation method, the warning, and what you did about it. Using part 1 of the example, with placeholders for the interval:

"A linear mixed model with a random intercept for group was fitted by REML with the lme4 package in R. The fit was singular: the between-group variance was estimated at zero (SD = 0.00; residual SD = 1.12). The estimated mean was 9.91 (SE = 0.18, 95% CI [LL, UL]). With only eight groups, we treat the zero estimate as a limit of the data and not as evidence that the groups do not differ."

Put the variance components in the table of model results next to the fixed effects, so that readers can see the zero. If you simplified the random effects, report the structure you started from and the one you ended with.

Related tools and guides

• Multilevel power and MDES calculator
• Fixed vs random effects in mixed models
• When to remove the intercept in regression
• REML or ML for likelihood ratio tests
• lme4 package (CRAN)

More answered questions

• REML or ML: which should a likelihood ratio test of fixed effects use?
• Which interaction terms belong in a regression or mixed model?
• When can you drop the intercept from a linear regression?

Written by AskStats with AI assistance. This is general information, not advice for your specific data or study. When the results matter, check your approach with a qualified statistician.