Your First Mixed Model with mixeff

The data

lme4::sleepstudy records reaction times (ms) for 18 subjects over ten days of sleep deprivation. Each subject was deprived of sleep starting on Day 0; by Day 9 most show substantial slowing.

sleep <- lme4::sleepstudy
str(sleep)
#> 'data.frame':    180 obs. of  3 variables:
#>  $ Reaction: num  250 259 251 321 357 ...
#>  $ Days    : num  0 1 2 3 4 5 6 7 8 9 ...
#>  $ Subject : Factor w/ 18 levels "308","309","310",..: 1 1 1 1 1 1 1 1 1 1 ...

The key structure: Subject appears 10 times in the data — once per day. Those 10 observations are not independent. They share the subject’s baseline speed, their individual sensitivity to sleep loss, and any other between-person trait that we haven’t measured. If you fit a plain OLS regression, those correlations inflate your Type I error rate and make standard errors too small.

A mixed model handles this by giving each subject its own intercept and its own slope — their personal baseline and their personal rate of slowing — and then estimating a population distribution over those person-level parameters.

library(ggplot2)
ggplot(sleep, aes(Days, Reaction, group = Subject)) +
  geom_line(alpha = 0.4) +
  geom_smooth(aes(group = 1), method = "lm", se = FALSE, colour = "steelblue") +
  labs(x = "Days of sleep deprivation", y = "Reaction time (ms)",
       title = "Individual subject trajectories") +
  theme_minimal()

The grey lines are individuals; the blue line is the population average. Subjects vary in their starting point and in how fast they slow down — which is exactly what a random-intercept-and-slope model captures.

Step 1: Compile the model

Before fitting, compile the model to see what mixeff understands about your formula. This separates formula interpretation from optimization; you can catch mis-specified random effects without paying the cost of a full fit.

spec <- compile_model(Reaction ~ Days + (Days | Subject), sleep)
explain_model(spec)
#> Random effects explanation:
#>   formula: Reaction ~ 1 + Days + (1 + Days | Subject)
#> 
#> Random effects:
#>   r0:
#>     wrote:      (Days | Subject)
#>     canonical:  (1 + Days | Subject)
#>     named form: re(group = Subject, intercept = TRUE, slopes = Days, cov = "full")
#>     scope:      `Subject` units differ in baseline and `Days` slope; the model estimates whether these are associated.
#>     covariance: full; theta parameters: 3
#>     support:    sufficient; group levels: 18; min rows/group: 10; median rows/group: 10
#>     variation:  Days=present; intercept=not_assessed

This pre-fit explanation confirms:

• one fixed effect (Days) plus an intercept,
• one random-effects block with a correlated intercept and slope per subject,
• the formula has been canonicalized to the explicit (1 + Days | Subject) form that the optimizer will receive.

Use audit(spec) when you want the deeper design audit. Its compact print starts with the audit summary and requested model; print(audit(spec), full = TRUE) shows the complete upstream report.

If you had written (1 | Subject) by mistake (random intercepts only), the explanation would show only one random-effects column per subject — useful to verify before a long fit.

Step 2: Fit

fit <- lmm(Reaction ~ Days + (Days | Subject), sleep)

lmm() returns an mm_lmm object. The actual optimization is performed by the bundled Rust engine; lmm() is the R entry point and result container.

Step 3: Read the summary

summary(fit, tests = "coefficients")
#> Linear mixed model fit by REML
#> Formula: Reaction ~ Days + (Days | Subject)
#> Fit status: converged_interior
#> 
#> Variance components:
#>    group        name variance  std_dev correlation
#>  Subject (Intercept) 612.0900 24.74050            
#>  Subject        Days  35.0718  5.92215       +0.07
#> Residual std. dev.: 25.5918
#> 
#> Fixed effects:
#>              Estimate Std. Error       df   t value  Pr(>|t|)        method
#> (Intercept) 251.40510   6.824557 17.00085 36.838304   < 1e-16 satterthwaite
#> Days         10.46729   1.545792 16.99989  6.771473 3.264e-06 satterthwaite
#> 
#> Inference status:
#>         term        method    status reliability
#>  (Intercept) satterthwaite available    moderate
#>         Days satterthwaite available    moderate
#>                             reliability_reason
#>  satterthwaite_finite_difference_approximation
#>  satterthwaite_finite_difference_approximation
#> 
#> Notes:
#>   Satterthwaite denominator df computed from finite-difference vcov_beta Jacobian and deviance Hessian over varpar

The summary has four blocks, in print order.

Fit status. Check this line first. converged_interior is the good outcome: the optimizer found a clean solution and every variance component is comfortably positive, so you can read the rest of the summary at face value. The status only changes when something needs your attention — for example, if a variance component collapses to zero (the data show no detectable variation for that term), the status reports a boundary fit, the affected row is flagged in the variance-components table, and the p-values below switch to a more conservative, clearly labeled method. In short: a healthy fit says converged_interior, and an unhealthy one tells you what went wrong instead of leaving you to notice.

Variance components. Subject (Intercept) is the between-subject spread in baseline reaction time (SD ≈ 25 ms); Subject Days is the between-subject spread in sensitivity to sleep loss (SD ≈ 6 ms per day); correlation is their association (+0.07 — essentially none). The residual standard deviation is the within-subject noise left over.

Fixed effects. Estimate is the population-level coefficient. Days = 10.47 means that, on average, reaction time increases by about 10.5 ms per day of sleep deprivation. The method column names how each p-value was computed — here satterthwaite, a finite-sample t test whose df column (≈ 17) is on the scale of the 18 subjects, not the 180 raw rows. mixeff never reports a number without naming the method behind it. On fits where Satterthwaite degrees of freedom cannot be computed (for example a variance component at the boundary), the summary shows clearly labeled asymptotic Wald z rows instead; the method column always tells you which one you got.

Inference status. One row per coefficient stating how much to trust the test: status says whether it was computed, reliability grades it, and reliability_reason names the engine’s warrant for that grade — here satterthwaite_finite_difference_approximation, meaning the degrees of freedom come from a finite-difference approximation, the standard route for this method (graded moderate). Any further engine notes print underneath. Grades and warrants are authored by the Rust engine, not by R-side heuristics.

Step 4: Extract components

The lme4-compatible extractors work on mm_lmm objects:

fixef(fit)
#> (Intercept)        Days 
#>   251.40510    10.46729

VarCorr(fit)
#> Variance components:
#>    group        name variance  std_dev correlation
#>  Subject (Intercept) 612.0900 24.74050            
#>  Subject        Days  35.0718  5.92215       +0.07
#> Residual std. dev.: 25.5918

confint(fit, method = "asymptotic")
#> Confidence intervals:
#>                 2.5 %    97.5 %
#> (Intercept) 238.02922 264.78099
#> Days          7.43759  13.49698
#> method: wald_asymptotic_from_stored_standard_errors
#> status: Wald (asymptotic) intervals from stored standard errors (engine-certified profile intervals: method = "profile")

confint() with method = "asymptotic" gives Wald intervals (fast; use "bootstrap" for small samples). The (Intercept) interval is the population mean reaction time on Day 0; the Days interval is the range of plausible slopes.

Step 5: Test a specific claim

contrast() evaluates a linear combination of fixed-effect coefficients — the mixeff equivalent of lme4::fixef() one-row hypothesis tests, but with method labelling.

ct <- contrast(fit, c("(Intercept)" = 0, Days = 1))
ct$table[, c("estimate", "std_error", "p_value", "method", "status", "reliability")]
#>   estimate std_error      p_value        method    status reliability
#> 1 10.46729  1.545792 3.263971e-06 satterthwaite available    moderate

The estimate is the Days slope; status = "available" and reliability = "certified" confirm that the covariance payload was present and the inference method was verified as appropriate for this fit. If a method cannot be certified (e.g. you request Kenward-Roger on a singular fit), reliability becomes "indicative" and reason names the problem.

Step 6: Compare specific conditions

mm_comparisons() makes it easy to answer “how much slower is a subject on Day 9 compared to Day 0?” without computing the arithmetic yourself.

cmp <- mm_comparisons(fit, specs = "Days", at = list(Days = c(0, 9)))
cmp$table[, c("label", "estimate", "conf_low", "conf_high", "p_value")]
#>             label estimate conf_low conf_high      p_value
#> 1 Days=9 - Days=0 94.20557 64.85354  123.5576 3.263971e-06

The at argument pins Days to exactly those two values. The difference Days=9 - Days=0 is 9 × the slope: about 94 ms, with a confidence interval.

Step 7: Save and reload

Mixed-model fits can take minutes to hours on large data sets. mixeff stores everything needed to revive the result later.

tmp <- tempfile(fileext = ".rds")
saveRDS(fit, tmp)

fit2 <- revive(readRDS(tmp))
stopifnot(isTRUE(all.equal(fixef(fit), fixef(fit2))))

revive() reconnects the R object to the Rust handle so the full inference surface is available. Plain readRDS() without revive() is sufficient when you only need the JSON-stored extractor values (coefficients, variance components, fitted values); call revive() when you want to run new contrast(), mm_means(), or other live-inference calls on the reloaded fit. See vignette("saving-and-reviving") for the full discussion.

Where to go next