Variance Inference Guide¶

After fitting a mixed model the natural questions are: are my variance components reliably estimated? and how uncertain are they? This page covers three tools for answering those questions:

is_singular / isSingular() — detect boundary fits before reporting
result.confint() — profile likelihood confidence intervals for variance parameters
VarCorr() — R-style variance-covariance summary

Boundary and singular fits¶

A singular fit occurs when one or more variance components are estimated at exactly zero — the boundary of the parameter space. This is not a crash or an error; it means the data contain little or no evidence for that grouping factor’s contribution. But boundary fits have practical consequences:

Standard errors for fixed effects may be inflated or poorly calibrated
Profile CIs for the collapsed component will hit zero on the lower bound
The model may be over-parameterised for the available data

Detecting a singular fit¶

import interlace
from interlace import isSingular

result = interlace.fit("score ~ hours_studied", data=df,
                       groups=["student_id", "school_id"])

# Function form (mirrors lme4)
if isSingular(result):
    print("Singular fit — check variance_components")

# Property shorthand
print(result.is_singular)          # True or False

# Per-factor flags
print(result.boundary_flags)
# {'student_id': False, 'school_id': True}

fit() issues a ConvergenceWarning automatically when a boundary fit is detected, so you will usually see the warning before you check manually.

What to do¶

Situation	Action
Near-zero variance for a grouping factor	Consider dropping it; compare AIC with `anova()`
Pre-registered random effect	Report the boundary result and flag the limitation
Convergence warning but non-zero variance	Try BOBYQA: `optimizer="bobyqa"`
Random slope variance collapsed	Try the independent (`

The tolerance for declaring a component “effectively zero” is tol=1e-4 (the lme4 default), matching the diagonal of the relative covariance factor Lambda_theta. Pass a custom tolerance to isSingular() if needed:

isSingular(result, tol=1e-5)   # stricter

Profile likelihood CIs for variance parameters¶

Wald confidence intervals (estimate ± 1.96 SE) assume asymptotic normality, which is a poor approximation for variance components — especially near the boundary or with few groups. Profile likelihood CIs are more accurate because they follow the actual curvature of the log-likelihood surface.

result.confint() computes these by fixing each variance parameter in turn, profiling out the others, and finding the points where 2 × (L_max − L(θ)) = χ²(level, df=1).

Basic usage¶

ci = result.confint()           # 95 % by default
print(ci)
#                 estimate   2.5 %  97.5 %
# school_id          0.412   0.201   0.731
# student_id         0.638   0.389   1.042
# residual           1.000   1.000   1.000

ci_90 = result.confint(level=0.90)

Interpreting the theta scale¶

CIs are reported on the theta scale — the diagonal of the relative Cholesky factor Lambda_theta. For an intercept-only random effect:

sigma_b ≈ theta × sqrt(sigma2_hat)

where sigma2_hat is result.scale. To convert:

sigma2 = result.scale
ci_sd = ci[["2.5 %", "97.5 %"]] * sigma2 ** 0.5
print(ci_sd)   # CIs on the standard-deviation scale

Boundary lower bounds¶

If the profile drops to its target before theta reaches zero, the lower bound is reported as 0.0. This is expected for singular or near-singular fits and means the data are compatible with the variance component being zero.

When to use profile vs Wald CIs¶

	Profile	Wald
Accuracy for variance params	Good	Poor near boundary
Speed	Slower (1D optimisation per param)	Instant
Appropriate for	Final reporting, small n_groups	Quick checks, large samples

`VarCorr()` — R-style variance-covariance summary¶

VarCorr() returns the variance components in the same format as R’s as.data.frame(VarCorr(fit)) — useful for direct comparison with lme4 output or for embedding in a tidy report.

from interlace import VarCorr

vc = VarCorr(result)
df_vc = vc.as_dataframe()
print(df_vc)
#         grp           var1   var2    vcov    sdcor
# 0  school_id  (Intercept)   None   0.412    0.642
# 1   Residual          None  None   1.284    1.133

Columns:

grp — grouping factor name (or "Residual")
var1, var2 — term names; var2 is non-null only for covariance entries in random-slope models
vcov — variance (diagonal) or covariance (off-diagonal)
sdcor — standard deviation (diagonal) or correlation (off-diagonal)

Random slopes¶

For a model with correlated random slopes, as_dataframe() includes both the variance and covariance rows:

result_slopes = interlace.fit(
    "rt ~ condition",
    data=df,
    random=["(1 + condition | subject)"],
)
print(VarCorr(result_slopes).as_dataframe())
#        grp           var1       var2    vcov   sdcor
# 0  subject  (Intercept)       None  45.20   6.723
# 1  subject  (Intercept)  condition  -8.40  -0.307
# 2  subject    condition       None   3.10   1.761
# 3  Residual        None       None  12.80   3.578

`VarCorr` vs `variance_components`¶

	`VarCorr().as_dataframe()`	`result.variance_components`
Format	Long-form DataFrame (R-compatible)	Dict keyed by group name
Covariances	Included	Not directly
SDs / correlations	Included as `sdcor`	Not directly
Best for	Reporting, R parity	Programmatic access

Putting it together¶

A complete post-fit variance diagnostic workflow:

import interlace
from interlace import isSingular, VarCorr

result = interlace.fit(
    "score ~ hours_studied + prior_gpa",
    data=df,
    groups=["student_id", "school_id"],
)

# 1. Check for boundary
if isSingular(result):
    print("Singular:", result.boundary_flags)

# 2. Inspect variance components in R-compatible format
print(VarCorr(result).as_dataframe())

# 3. Profile CIs (skip if singular — lower bound will be 0)
if not result.is_singular:
    ci = result.confint()
    print(ci)