Mixed-Effects Models Exercises in R: 18 Practice Problems

Explanation: The notation (1 | Subject) adds a separate intercept deviation for every level of Subject, while the Days coefficient stays a single population-level slope. The Subject standard deviation of 37.1 ms is far larger than zero, signalling that ignoring subject-to-subject baseline differences (i.e., a plain lm()) would underestimate uncertainty for the Days effect. REML is the default estimator because it gives less biased variance components than ML when group counts are modest.

Explanation: VarCorr() returns a list with one entry per grouping factor plus a residual scale stored as the "sc" attribute. Squaring the standard deviations gives the variances on the original response scale. An ICC near 0.59 means almost 60% of the unexplained variance after accounting for Days is between-subject, so any inference that ignores grouping (pseudo-replication) will badly understate standard errors. The performance::icc() helper computes the same quantity in one call.

Explanation: ranef() returns shrinkage estimates: subject deviations are pulled toward zero in proportion to how little data each subject contributes. These are not the same as fitting 18 separate intercepts with lm(): BLUPs borrow strength across subjects and produce smaller mean squared error for new predictions. Coercing to a tibble keeps the subject IDs as a column instead of leaving them buried in rownames, which is the form most downstream plotting and joining steps want.

Explanation: (Days | Subject) is shorthand for (1 + Days | Subject): an intercept and a slope per subject, with a covariance between them estimated from the data. The residual SD drops from 31 to 26 ms because per-subject variation in the rate of decline is no longer absorbed into the noise term. The fixed-effect estimates (251.4, 10.47) are unchanged because the design is balanced; in unbalanced data they would shift.

Explanation: The random-effects correlation matrix is attached as an attribute on each grouping factor's variance-covariance matrix; the off-diagonal element [1, 2] is the corr between intercept and Days. A value near 0.07 is effectively zero given the small subject count, meaning baseline reaction time tells you little about how steeply someone deteriorates. If this correlation were estimated as exactly ±1 it would signal a singular fit; consider (Days || Subject) (the double-bar form, which forces zero correlation) or refitting with a simpler random structure.

Explanation: A likelihood ratio test for nested models requires ML, not REML, because REML log-likelihoods compare differently parameterised mean structures incorrectly. anova.merMod() quietly refits both models with REML = FALSE before computing the chi-square. Note that the LRT is known to be conservative on the boundary of the parameter space (testing whether a variance is zero), so the true p-value is likely even smaller; for borderline cases use parametric bootstrap with pbkrtest::PBmodcomp().

Explanation: Each Chick only ever receives one Diet, so Chick is implicitly nested inside Diet; you do not need Diet/Chick in the random part because the diet effect is captured as a fixed factor. The near-perfect negative correlation (-0.99) between random intercept and slope is a classic warning sign: the model is on the boundary, often because the time variable is not centred. Re-running with Time centred at its mean usually returns a sensible correlation and a non-degenerate covariance matrix.

Explanation: The shorthand school_id/class_id expands to school_id + school_id:class_id, which creates the interaction grouping factor needed when class labels are reused across schools. Writing (1 | school_id) + (1 | class_id) instead would tell lmer() that class "1" is the same unit in every school, which is wrong. The cleaner long-term fix is to mint globally unique class IDs (paste(school_id, class_id, sep = "_")) and use (1 | school_id) + (1 | unique_class) so the structure is explicit in the data, not implicit in the formula.

Explanation: Crossed effects appear when every level of one grouping factor (subjects) is observed with every level of another (items). Treating items as a fixed factor would burn 19 degrees of freedom and prevent generalization to new items; treating them as random instead lets you generalize the condition effect to the broader population of stimuli. Barr et al.'s "keep it maximal" advice would also add random slopes for condition by both subject and item, but that often fails to converge with small designs.

Explanation: REML estimates variance components after profiling out the fixed effects, so the REML log-likelihood depends on the fixed-effect design matrix and is not comparable between models with different fixed parts. Switching to ML with REML = FALSE fixes this. For random-effect tests (testing whether a variance is zero) the opposite advice applies: use REML, because ML biases variance estimates downward, especially with few groups.

Explanation: lmerTest masks lme4::lmer() and attaches Satterthwaite-approximate degrees of freedom and t-tests to the summary. Df of 17 reflects the 18-subject design (n_groups - 1). Use Kenward-Roger (lmerTest::lmer(..., ddf = "Kenward-Roger")) for smaller designs where Satterthwaite under-covers; it is slower but typically more accurate. Loading lmerTest after lme4 is the canonical order, since you want the p-value-aware lmer() to win the name conflict.

Explanation: Parametric bootstrap simulates new responses from the fitted model and refits each replicate; the empirical quantiles of the replicate coefficient form the CI. This is more honest than Wald intervals when the sampling distribution is asymmetric, which often happens for variance components and GLMM coefficients. The seed argument makes the run reproducible; without it, two readers running the same code see slightly different CIs.

Explanation: A patient who is "compliant" tends to stay compliant across all 8 visits; a random intercept absorbs that within-person correlation so the arm coefficient gets honest standard errors. glmer defaults to the Laplace approximation, fast and adequate for binary outcomes with reasonable group sizes. For more accuracy use nAGQ = 10 or higher (adaptive Gauss-Hermite quadrature), but it only works for single-grouping-factor models.

Explanation: Including log(read_total) as an offset (coefficient fixed at 1) reparameterises counts as rates: the linear predictor returns log(hits / read_total). Without the offset, deeply sequenced samples would have systematically higher counts and the tissue coefficients would absorb that noise. The random intercept on sample handles residual library-prep variation beyond what depth captures; with only one observation per sample it is identified by the Poisson mean-variance link.

Explanation: A dispersion of about 1.0 says the Poisson mean-variance assumption fits the data, partly because the (1 | sample) random intercept acts as an observation-level random effect that already soaks up extra variance. If this number came back at 3 you would refit with glmer.nb() (negative binomial) or add an explicit observation-level random effect like (1 | obs_id) where obs_id is unique per row. performance::check_overdispersion() automates the whole check.

Explanation: Mixed-effects fits return residuals conditional on the random effects, so a flat scatter centred on zero is the diagnostic of choice. A funnel shape (variance growing with the mean) would suggest a log-transform on the response; a curved trend would suggest a non-linear time effect, e.g., a spline on Days. DHARMa::simulateResiduals() produces simulation-based residuals that are uniform on [0, 1] under a correct model and is more reliable for GLMM diagnostics.

Explanation: Setting re.form = NA (or equivalently ~0) zeros out random effects, giving the marginal expectation: what the average subject would do. Using re.form = NULL (the default) would still need a Subject value in newdata to add the subject-specific deviations. For uncertainty around these population predictions, use merTools::predictInterval() or a parametric bootstrap, since predict.merMod() does not return standard errors.

Explanation: Nakagawa-Schielzeth marginal R² (R2m) is the share of variance the fixed predictors alone explain; conditional R² (R2c) adds the variance soaked up by random effects. The gap (about 52 percentage points here) is the share captured by subject-to-subject heterogeneity, not by Days. For GLMMs you get two flavours of each (theoretical and delta-method), so the function returns a 2x2 matrix; the theoretical row is usually the one to report.