Volume 5 · twelve days

Your rows are not independent.

Diary entries within people, trials within participants, pupils within classes, patients within therapists: behavioural data is nested almost by default, and a linear model that ignores it reports standard errors that are too small and p-values that are too impressive. Mixed-effects models handle it natively — unbalanced designs, missing occasions and all — and this volume takes you from the ICC to growth curves and a full ESM analysis.

12
days, ICC to growth curves
2
sources of variance in every model
0
participants dropped for one missing occasion
Day 41 45 minutes · the error that inflates every p

Dependence, and the size of the damage

Ordinary regression assumes each row contributes one independent piece of information. When forty participants each give twenty diary entries, you have 800 rows but nowhere near 800 independent observations — people are more like themselves than like each other. Ignore that and your standard errors shrink by a factor that depends on cluster size and the intraclass correlation. The false-positive rate can rise from 5% to over 40%.

# Simulate data with NO true effect, but strong clustering. set.seed(7) sim_one <- function() { id <- rep(1:40, each = 20) u <- rnorm(40, sd = 1.2)[id] # person intercepts: the clustering x <- rnorm(800) # predictor unrelated to y y <- 0 * x + u + rnorm(800, sd = 1) # true slope is exactly zero coef(summary(lm(y ~ x)))["x", "Pr(>|t|)"] } ps <- replicate(500, sim_one()) mean(ps < .05) # should be .05 if inference were valid#> [1] 0.052 # ...now make the predictor a person-level variable, which is the common case: # x <- rnorm(40)[id] #> [1] 0.416 ← 42% false positives at a nominal 5%
Pseudo-replication, named
Treating 800 correlated rows as 800 independent observations is called pseudo-replication, and it is the single most common fatal flaw in student ESM and trial-level analyses. It does not always inflate error — for purely within-cluster predictors the damage can be small — but for between-cluster predictors (condition assigned per person, therapist characteristics, classroom-level interventions) it is severe, because the effective n is the number of clusters, not the number of rows.
repeated measures
Same person measured several times: trials, sessions, diary beeps, days.
hierarchical sampling
Pupils in classes in schools; patients within therapists within clinics.
crossed designs
Every participant sees every stimulus: two crossing sources of dependence.
longitudinal growth
Repeated measures where time itself is the predictor of interest.
Old solutionWhy mixed models replaced it
Aggregate to one row per personThrows away within-person variance — usually the theoretically interesting part — and weights a person with 3 observations like one with 30.
Repeated-measures ANOVARequires complete cases and sphericity; cannot take continuous time or person-level moderators comfortably.
Cluster-robust standard errorsA defensible fix for inference, but gives you no variance components, no random slopes and no cluster-level predictions.
Do this now · 15 minutes
Write down your data's nesting structure explicitly: what is a row, what is a cluster, how many clusters, how many rows per cluster, and which predictors vary within versus between clusters. That sentence determines every model in this volume.
Day 42 55 minutes · one new piece of syntax

lmer(): the random intercept, and how to read the output

The whole syntactic addition is a term in parentheses. (1 | id) means "let every participant have their own intercept, drawn from a normal distribution whose variance I want estimated". Everything else — the formula, the predictors, the interpretation of fixed effects — is what you already learned in volume 4.

library(lme4); library(lmerTest) # lmerTest adds p-values to the summary m0 <- lmer(mood ~ 1 + (1 | id), data = esm) # intercept-only: the ICC model m1 <- lmer(mood ~ stress + (1 | id), data = esm) # one fixed predictor summary(m1)#> Linear mixed model fit by REML. t-tests use Satterthwaite's method #> Formula: mood ~ stress + (1 | id) #> #> Random effects: #> Groups Name Variance Std.Dev. #> id (Intercept) 0.6412 0.8008 #> Residual 0.9033 0.9504 #> Number of obs: 2814, groups: id, 94 #> #> Fixed effects: #> Estimate Std. Error df t value Pr(>|t|) #> (Intercept) 4.81204 0.09112 92.41000 52.809 < 2e-16 *** #> stress -0.31182 0.01884 2716.30 -16.55 < 2e-16 ***
LineWhat it tells you
id (Intercept) 0.641How much people differ from each other in average mood. This is variance you previously threw into the error term.
Residual 0.903Moment-to-moment variance within a person, after the fixed effects.
groups: id, 94Your effective sample size for person-level predictors is 94, not 2,814. Check this number every time.
stress -0.312One point more stress goes with 0.31 less mood — a within-person association here, because stress varies within people. Day 44 makes that explicit.
df 2716.3Fractional because Satterthwaite approximates the denominator df. Mixed models have no single exact df.
Fixed or random?
A factor is fixed if its specific levels are what you want to draw conclusions about (treatment versus control; three explicit time points) and random if the levels are an interchangeable sample from a population you want to generalise over (these 94 participants, these 32 stimuli, these 18 therapists). The practical test: would replacing the levels with a different sample of the same kind leave your hypothesis intact? If yes, random.
Do this now · 25 minutes
Refit the repeated-measures analysis from day 32 as lmer(y ~ time + (1 | id)). Compare the number of participants used against what aov_ez() dropped. That difference is the practical argument for this volume.
Day 43 55 minutes · partitioning variance

The ICC, and letting effects differ by person

The intraclass correlation is the share of total variance that lies between clusters — how much of the outcome is explained by who you are rather than when you were measured. It tells you how badly clustering would have bitten, and it is a reportable descriptive in its own right. The next step is allowing not just the intercept but the slope to vary: some people react more strongly to stress than others, and that variation is often the hypothesis.

library(performance) icc(m0) # from the intercept-only model # By hand, so the formula is not magic: # ICC = between-person variance / (between + within) 0.6412 / (0.6412 + 0.9033) VarCorr(m1) # variance components with correlations#> # Intraclass Correlation Coefficient #> Adjusted ICC: 0.415 #> Unadjusted ICC: 0.415 #> [1] 0.4151 #> #> ← 42% of mood variance is between people. Substantial clustering; lm was never an option. # Random slopes: (1 + stress | id) = intercept AND stress effect vary by person m2 <- lmer(mood ~ stress + (1 + stress | id), data = esm) summary(m2) anova(m1, m2, refit = FALSE) # is the slope variance worth estimating? # Per-person deviations, for plotting library(tidyverse) ranef(m2)$id |> as_tibble(rownames = "id") |> head()#> Random effects: #> Groups Name Variance Std.Dev. Corr #> id (Intercept) 0.71204 0.8438 #> stress 0.02183 0.1478 -0.31 #> Residual 0.84121 0.9172 #> #> Models: #> m1: mood ~ stress + (1 | id) #> m2: mood ~ stress + (1 + stress | id) #> npar AIC BIC logLik deviance Chisq Df Pr(>Chisq) #> m1 4 8214.3 8238.1 -4103.2 8206.3 #> m2 6 8163.7 8199.4 -4075.8 8151.7 54.712 2 1.316e-12 ***
Reading a slope SD, and the correlation term
An SD of 0.148 on the stress slope means that around the average effect of −0.31, individual slopes span roughly −0.61 to −0.02 (±2 SD): almost everyone responds negatively to stress, but the strength varies threefold. The Corr = -0.31 says people with higher average mood have more negative stress slopes — a substantive finding sitting in the random-effects table that most papers never mention. Report it when it is interpretable.
# The figure that makes random slopes obvious library(ggplot2) esm |> ggplot(aes(stress, mood, group = id)) + geom_line(stat = "smooth", method = "lm", se = FALSE, alpha = 0.18, linewidth = 0.4) + geom_smooth(aes(group = 1), method = "lm", colour = "#e26d5c", linewidth = 1.3) + labs(x = "Momentary stress", y = "Momentary mood", caption = "Thin lines: 94 individual slopes. Thick line: fixed effect.")
Do this now · 25 minutes
Compute the ICC for your own outcome and report it as a sentence. Then fit the random-slope version of your key predictor, compare with anova(..., refit = FALSE), and produce the spaghetti figure.
Day 44 55 minutes · the confound hiding in plain sight

Within-person and between-person effects are different effects

"Stress predicts mood" is two claims. On days you are more stressed than usual, is your mood worse (a within-person process)? And are generally more stressed people generally in worse moods (a between-person difference)? A raw predictor in a mixed model returns an uninterpretable blend of the two. Splitting it is one mutate() and it changes what your paper is about.

library(tidyverse) esm <- esm |> group_by(id) |> mutate(stress_pm = mean(stress, na.rm = TRUE), # person mean: between stress_cw = stress - stress_pm) |> # deviation: within ungroup() |> mutate(stress_pm_c = stress_pm - mean(stress_pm, na.rm = TRUE)) m3 <- lmer(mood ~ stress_cw + stress_pm_c + (1 + stress_cw | id), data = esm) summary(m3)#> Fixed effects: #> Estimate Std. Error df t value Pr(>|t|) #> (Intercept) 4.79812 0.08914 91.204 53.827 < 2e-16 *** #> stress_cw -0.28714 0.02012 88.771 -14.272 < 2e-16 *** #> stress_pm_c -0.51028 0.11841 92.013 -4.310 4.16e-05 *** #> #> ← the two effects differ in size by nearly a factor of two, and answer #> two different research questions
TermQuestion it answers
stress_cwWhen this person is more stressed than their own average, is their mood worse? The process claim. Cannot be confounded by stable person characteristics.
stress_pm_cDo more-stressed people have worse mood on average? A between-person, cross-sectional claim with all the usual confounding.
Raw stress, uncentredA weighted mixture of the two, weights depending on the ICC and cluster sizes. Interpretable as neither.
Why this is the single highest-value hour in the volume
Reviewers of intensive longitudinal work check for this first. The within-person coefficient is the one your theory almost always predicts, and it is the one immune to trait-level confounding — which is the main methodological selling point of ESM designs. Reporting a raw uncentred coefficient throws away that advantage and invites the comment "the authors conflate within- and between-person effects". Grand-mean-centre the person means too, so the intercept is the average person's average mood.
# Interactions across levels: does a trait moderate a within-person process? m4 <- lmer(mood ~ stress_cw * neuroticism_c + stress_pm_c + (1 + stress_cw | id), data = esm) library(ggeffects) ggpredict(m4, terms = c("stress_cw", "neuroticism_c [-1, 0, 1]")) |> plot()#> Fixed effects: #> Estimate Std. Error df t value Pr(>|t|) #> stress_cw -0.28812 0.01984 87.412 -14.521 < 2e-16 *** #> neuroticism_c -0.31204 0.09812 91.883 -3.181 0.00202 ** #> stress_cw:neuroticism_c -0.08114 0.02013 88.102 -4.031 0.00012 ***
Do this now · 25 minutes
Split your key time-varying predictor into person-mean and within-person deviation, refit, and write both coefficients as separate sentences. Then decide which one your research question was actually about — and say so in the paper.
Day 45 55 minutes · REML, ML, and warnings

Model comparison, convergence, and singular fits

Two practical matters that consume more PhD hours than the statistics: which estimator to compare models under, and what to do when lmer() prints a warning. REML gives better variance estimates and is the default; comparing models that differ in fixed effects requires ML, so set REML = FALSE or let anova() refit for you.

# Comparing RANDOM structures: keep REML, do not refit anova(m1, m2, refit = FALSE) # Comparing FIXED structures: needs ML m_a <- lmer(mood ~ stress_cw + (1 | id), data = esm, REML = FALSE) m_b <- lmer(mood ~ stress_cw + sleep_cw + (1 | id), data = esm, REML = FALSE) anova(m_a, m_b) # Marginal and conditional R-squared, for the reader who asks performance::r2(m3)#> Data: esm #> Models: #> m_a: mood ~ stress_cw + (1 | id) #> m_b: mood ~ stress_cw + sleep_cw + (1 | id) #> npar AIC BIC logLik deviance Chisq Df Pr(>Chisq) #> m_a 4 8214.3 8238.1 -4103.2 8206.3 #> m_b 5 8189.7 8219.4 -4089.9 8179.7 26.712 1 2.36e-07 *** #> #> # R2 for Mixed Models #> Conditional R2: 0.512 ← fixed + random #> Marginal R2: 0.183 ← fixed effects only
WarningWhat it means and what to do
boundary (singular) fitA variance component is estimated at zero (or a correlation at ±1): the data cannot support that random term. Simplify — drop the correlation with ||, or the slope. It is not an error, but reporting a model with a singular fit unremarked is.
Model failed to convergeThe optimiser stopped unhappy. Rescale predictors to similar ranges first (this fixes most cases), then try allFit() to see whether the estimates are stable across optimisers.
nearly unidentifiableUsually collinear predictors or too few clusters for the random structure. Check cluster count per random term.
No warning at allStill check: 5 clusters cannot support a random slope, warning or not.
# When it will not converge: try every optimiser before simplifying library(lme4) af <- allFit(m2) summary(af)$fixef # are the fixed effects the same everywhere? then relax # A common, effective fix m2b <- lmer(mood ~ stress_cw + (1 + stress_cw || id), data = esm) # no correlation control <- lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e5))#> bobyqa : [OK] #> Nelder_Mead : [OK] #> nlminbwrap : [OK] #> nmkbw : [OK] #> optimx.L-BFGS-B : [OK] #> (Intercept) stress_cw #> bobyqa 4.7981 -0.2871 #> Nelder_Mead 4.7981 -0.2871 #> ← identical to four decimals: the warning was cosmetic
Two R-squareds, and why you report both
Marginal R² is the variance explained by fixed effects alone — what a new participant's data would be predicted from. Conditional R² adds the random effects, so it includes "knowing who this person is". The gap between them is a direct measure of how much of your outcome is stable individual difference. Report both, and never report only the conditional one as if it were your model's explanatory power.
Do this now · 25 minutes
Take your most complex model and deliberately over-specify it until it warns. Then work back through the checklist — rescale, drop the correlation, drop the slope — recording what each step changed. You now have the debugging routine you will use for years.
Day 46 50 minutes · why they are hard here

Degrees of freedom, and four defensible options

lme4 deliberately does not print p-values, because in an unbalanced mixed model there is no exact denominator degrees of freedom. That is honest rather than obstructive. Four approaches are accepted, in roughly ascending order of rigour and computation time.

ApproachWhen to use it
Satterthwaite (lmerTest)The practical default. Load lmerTest and summary() gains p-values. Well-calibrated for typical designs.
Kenward-Roger (afex, pbkrtest)Better with few clusters (under ~30) or complex random structures. Slower.
Likelihood-ratio testFor comparing nested fixed structures; needs ML fits. Anti-conservative with few clusters.
Parametric bootstrap (pbkrtest, confint(method="boot"))Gold standard, expensive. Use for the one or two effects your paper rests on.
|t| > 2 as a rule of thumbAcceptable only as a screening heuristic with many clusters; never as reported inference.
library(lmerTest) # Satterthwaite in summary() summary(m3)$coefficients library(afex) mixed(mood ~ stress_cw + stress_pm_c + (1 + stress_cw | id), data = esm, method = "KR") # Kenward-Roger # Bootstrapped CIs for the effects you care most about confint(m3, method = "boot", nsim = 1000, parm = "beta_")#> Mixed Model Anova Table (Type 3 tests, KR-method) #> Effect df F p.value #> 1 stress_cw 1, 88.77 203.68 *** <.001 #> 2 stress_pm_c 1, 92.01 18.57 *** <.001 #> #> Computing bootstrap confidence intervals ... #> 2.5 % 97.5 % #> (Intercept) 4.6212041 4.9740182 #> stress_cw -0.3271408 -0.2478213 #> stress_pm_c -0.7412083 -0.2801442
What to report
Give the estimate, its 95% CI, and the test with its degrees of freedom and method named: "b = −0.29, 95% CI [−0.33, −0.25], t(88.8) = −14.27, p < .001 (Satterthwaite)". Also report the number of clusters and observations, and the random-effects structure in full — a mixed-model result is not reproducible without it. Volume 9 day 87 builds the table.
Do this now · 20 minutes
Get p-values for your main effect via Satterthwaite and Kenward-Roger, and a bootstrapped CI for the same coefficient. If they disagree materially, you have too few clusters for the model — note that as a limitation and prefer the bootstrap.
Day 47 55 minutes · participants AND stimuli

Two crossing populations: the stimulus-as-fixed-effect fallacy

In a lexical decision task every participant sees every word. Participants are a sample from a population of people; words are a sample from a population of words. Treating items as fixed means your conclusions apply only to those 32 words — and inflates false positives sharply. Crossed random effects fix it with one more term in parentheses.

# Both sources of sampling, crossed (not nested) m_cross <- lmer(rt ~ condition + (1 + condition | subject) + (1 | item), data = trials) summary(m_cross)#> Random effects: #> Groups Name Variance Std.Dev. Corr #> subject (Intercept) 2104.31 45.873 #> conditionrelated 381.22 19.525 -0.18 #> item (Intercept) 892.04 29.867 #> Residual 7218.66 84.963 #> Number of obs: 6144, groups: subject, 48; item, 32 #> #> Fixed effects: #> Estimate Std. Error df t value Pr(>|t|) #> conditionrelated -31.4082 4.9121 47.2100 -6.394 6.71e-08 ***
Maximal, then simplify — carefully
The influential advice is to keep the random structure "maximal": random slopes for every effect that varies within each grouping factor. In practice maximal models often fail to converge with realistic trial counts. The defensible compromise: specify the maximal structure your design justifies, and if it will not fit, simplify in a preregistered order — drop random correlations first (||), then the slopes with the smallest estimated variance — and report exactly what you dropped and why. Never simplify in the direction that makes your effect significant.
# Nested vs crossed: get this wrong and the model is silently mis-specified # Pupils in classes in schools — NESTED (class 1 in school A ≠ class 1 in school B) lmer(score ~ treat + (1 | school/class), data = d) # = (1|school) + (1|school:class) # Participants and stimuli — CROSSED (item 1 is the same item for everyone) lmer(rt ~ cond + (1 | subject) + (1 | item), data = d) # Are my grouping factors crossed or nested? Ask the data. with(trials, table(subject, item)) |> (\(t) mean(t > 0))() # 1 = fully crossed#> [1] 1 ← every subject saw every item: fully crossed
Do this now · 25 minutes
If you have trial-level data, fit it with and without (1 | item) and compare the standard error of your main effect. The difference is the size of the fallacy you would have committed. If you have hierarchical sampling instead, verify with a cross-tabulation that your nesting is real.
Day 48 55 minutes · binary and count outcomes, clustered

glmer: accuracy, relapse and counts within people

Trial-level accuracy, day-level binge occurrence, weekly session counts — binary and count outcomes that are also clustered. glmer() combines volume 4 day 40 with everything in this volume. Two warnings: it is much slower and more fragile than lmer(), and coefficients are on the link scale, so report predicted probabilities as well.

library(lme4) m_bin <- glmer(correct ~ condition + (1 + condition | subject) + (1 | item), data = trials, family = binomial, control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e5))) summary(m_bin) exp(fixef(m_bin)) # odds ratios library(marginaleffects) avg_comparisons(m_bin, variables = "condition") # in probability, for readers#> Random effects: #> Groups Name Variance Std.Dev. Corr #> subject (Intercept) 0.41204 0.6419 #> conditionrelated 0.08812 0.2969 -0.24 #> item (Intercept) 0.22104 0.4701 #> Number of obs: 6144, groups: subject, 48; item, 32 #> #> Fixed effects: #> Estimate Std. Error z value Pr(>|z|) #> conditionrelated 0.51204 0.11842 4.323 1.54e-05 *** #> #> Term Contrast Estimate Std. Error z 2.5 % 97.5 % #> condition related - unrelated 0.062 0.014 4.43 0.035 0.089 #> ← 6.2 percentage points more accurate # Counts, and the check you must not skip m_cnt <- glmer(episodes ~ treatment + (1 | id), data = diary, family = poisson) performance::check_overdispersion(m_cnt) # Overdispersed (usual for behavioural counts): negative binomial m_nb <- glmmTMB::glmmTMB(episodes ~ treatment + (1 | id), data = diary, family = nbinom2) # Many structural zeros (days with no episodes at all)? m_zi <- glmmTMB::glmmTMB(episodes ~ treatment + (1 | id), ziformula = ~ 1, data = diary, family = nbinom2) AIC(m_cnt, m_nb, m_zi) # Exposure differs between people (different numbers of observed days)? # offset(log(days_observed)) — models a RATE, not a count#> # Overdispersion test #> dispersion ratio = 3.812 #> p-value = < 0.001 ← Poisson is wrong; SEs far too small #> #> df AIC #> m_cnt 3 4218.412 #> m_nb 4 3611.208 #> m_zi 5 3588.114 ← zero-inflated NB wins clearly
Three GLMM realities
No residual variance term — the binomial and Poisson distributions fix it, which is why overdispersion has to be checked separately. Convergence is harder: rescale predictors, use bobyqa, and expect to simplify the random structure more than with lmer. Effect sizes on the link scale are not comparable across models with different random structures; report marginal predictions in probability or rate units instead.
Do this now · 25 minutes
Fit a binary or count GLMM on your own data, check overdispersion if it is a count, and report both the odds/rate ratio and the marginal effect in natural units. Note the fitting time — that number decides how you plan day 80's power simulation.
Day 49 60 minutes · the design of the decade

ESM: structure, lags, and within-person dynamics

Experience sampling is where multilevel modelling earns its keep, and where the data-handling is hardest: beeps nested in days nested in people, missing responses that are not missing at random, and questions about temporal order that require lagged predictors built within person and within day. Get the lag construction wrong and you have silently modelled yesterday's last beep predicting today's first.

library(tidyverse) esm <- esm |> arrange(id, day, beep) |> group_by(id, day) |> # NOT just group_by(id) mutate(stress_lag = lag(stress), # previous beep, same day mood_lag = lag(mood), gap_mins = as.numeric(difftime(time, lag(time), units = "mins"))) |> ungroup() |> filter(is.na(gap_mins) | gap_mins < 240) # drop implausible intervals # Compliance and structure — report these before any model esm |> group_by(id) |> summarise(n_beeps = sum(!is.na(mood))) |> summarise(median = median(n_beeps), min = min(n_beeps), max = max(n_beeps))#> # A tibble: 1 × 3 #> median min max #> <dbl> <int> <int> #> 31 8 42 #> ← report the range: a participant with 8 of 42 beeps contributes little #> and may differ systematically # The canonical ESM model: does stress predict mood at the NEXT beep, # controlling for current mood (autoregression)? m_esm <- lmer(mood ~ stress_lag_cw + mood_lag_cw + stress_pm_c + (1 + stress_lag_cw | id), data = esm) summary(m_esm)#> Fixed effects: #> Estimate Std. Error df t value Pr(>|t|) #> (Intercept) 4.81042 0.08814 90.412 54.578 < 2e-16 *** #> stress_lag_cw -0.11204 0.02184 87.204 -5.130 1.91e-06 *** #> mood_lag_cw 0.28114 0.02012 2402.11 13.974 < 2e-16 *** #> stress_pm_c -0.48812 0.11204 91.883 -4.357 3.61e-05 *** #> #> ← controlling for mood inertia, a within-person stress spike still predicts #> worse mood ~2 hours later
ESM decisionThe defensible choice
Two or three levels?Beeps in days in people is three levels, but (1 | id) + (1 | id:day) is often enough. Fit both; report the one you preregistered.
Unequal intervalsState the average interval and whether beeps were fixed or random. Consider continuous-time models (ctsem) if intervals vary wildly.
Missing beepsReport compliance per person; test whether missingness relates to previous mood. Mixed models tolerate missing occasions but not informative missingness.
AutoregressionInclude mood_lag unless you have a reason not to — without it, a lagged effect may just be mood's own persistence.
Time of day / day of studyControl for both. Diurnal cycles and reactivity effects are real and will otherwise sit in your effect.
Do this now · 30 minutes
Build lagged versions of your key variables grouped by person and day, inspect the gap distribution, and fit the lagged model with an autoregressive term. Report compliance before results — reviewers of ESM papers look for it immediately.
Day 50 60 minutes · change as the outcome

Growth models: time coding is a substantive choice

When the question is how people change, time becomes a predictor with a random slope, and the coding of time decides what your intercept means. Centre time at baseline and the intercept is the starting level; centre at the midpoint and it is the average level. Neither is wrong, but only one matches your hypothesis, and the intercept-slope correlation — do people who start higher change faster? — is often the interesting result.

# Time coded so 0 = baseline: intercept = starting level long <- long |> mutate(t = wave - 1, # 0, 1, 2, 3 t_sq = t^2) g1 <- lmer(symptoms ~ t + (1 + t | id), data = long) # linear growth g2 <- lmer(symptoms ~ t + t_sq + (1 + t | id), data = long) # curvature anova(g1, g2, refit = FALSE) summary(g1)#> Random effects: #> Groups Name Variance Std.Dev. Corr #> id (Intercept) 18.4120 4.2909 #> t 1.8812 1.3716 -0.42 #> Residual 6.2104 2.4921 #> #> Fixed effects: #> Estimate Std. Error df t value Pr(>|t|) #> (Intercept) 24.81204 0.48121 118.4100 51.562 < 2e-16 *** #> t -2.41082 0.16204 117.2100 -14.877 < 2e-16 *** #> #> ← average 2.41-point decline per wave; SD of slopes 1.37, so some people #> improve three times faster than others, and higher starters improve faster #> (r = -0.42) # Does treatment change the RATE of change? A condition × time interaction. g3 <- lmer(symptoms ~ t * condition + (1 + t | id), data = long) library(emmeans) emtrends(g3, ~ condition, var = "t") # slope per condition, with CIs pairs(emtrends(g3, ~ condition, var = "t")) # difference in slopes = the effect library(ggeffects) ggpredict(g3, terms = c("t", "condition")) |> plot(add.data = TRUE)#> condition t.trend SE df lower.CL upper.CL #> control -1.204 0.2114 117 -1.623 -0.785 #> CBT -3.612 0.2087 117 -4.026 -3.198 #> #> contrast estimate SE df t.ratio p.value #> control - CBT 2.408 0.2971 117 8.105 <.0001
Four decisions to make before fitting
Where is time zero? It defines the intercept. Is time in waves or real elapsed days? With variable assessment dates, real time is more honest and mixed models handle it natively. Linear or curved? Compare, but do not fit a quadratic to three waves — it is saturated. Random slope? Almost always yes: assuming everyone changes at the same rate is exactly the assumption growth modelling exists to relax.
Do this now · 30 minutes
Fit a growth model with time centred at baseline, report the average slope with its CI and the SD of individual slopes, then add your condition-by-time interaction and report the difference in slopes with emtrends(). That contrast is the treatment effect.
Day 51 50 minutes · the table a reviewer accepts

Reporting: fixed effects, variance components, and a figure

A mixed-model result is not reproducible from fixed effects alone. The reporting standard is: the full model formula, sample size at every level, fixed effects with CIs and the df method named, the variance components (and correlations), the estimator, and marginal plus conditional R². Then one figure of model-implied values.

library(sjPlot) tab_model(m3, m_esm, show.ci = 0.95, show.re.var = TRUE, show.icc = TRUE, show.obs = TRUE, show.ngroups = TRUE, dv.labels = c("Momentary mood", "Mood, lagged model"), file = "output/tables/table3_mixed.html") # Alternative, more controllable, straight to Word or LaTeX library(modelsummary) modelsummary(list("Random intercept" = m1, "Random slope" = m2, "Centred" = m3), statistic = "conf.int", gof_map = c("nobs", "icc", "r2.marginal"), output = "output/tables/table3.docx") library(ggeffects); library(ggplot2) p <- ggpredict(m3, terms = "stress_cw [all]") |> plot(add.data = TRUE, dot.alpha = 0.06) + labs(x = "Within-person stress (deviation from own mean)", y = "Predicted mood", title = NULL, caption = "Model-implied values with 95% CI. 2,814 observations from 94 participants.") + theme_paper() ggsave("output/figures/fig3_mixed.png", p, width = 85, height = 70, units = "mm", dpi = 600, bg = "white")
The reporting checklist
1. Full formula, in the text or a table note. 2. N at each level ("2,814 observations from 94 participants, median 31 each, range 8–42"). 3. Fixed effects: estimate, SE, 95% CI, t/z, df and method. 4. Variance components and any random correlations. 5. Estimator (REML) and software versions — renv and sessionInfo(), volume 9 day 92. 6. Marginal and conditional R². 7. Any simplification of the random structure, and why. 8. How missing data were handled. Miss item 7 and a good reviewer will assume the worst.
Do this now · 20 minutes
Generate your mixed-model table with tab_model() and the marginal-effects figure, then check every one of the eight checklist items against your own draft. Write the missing ones.
Day 52 60 minutes · the volume's deliverable

One pipeline: raw diary export to results section

This day is not new content; it is the integration that proves the volume. Take a raw platform export and produce a results section with two tables and two figures, from a script that reruns end to end. If you have no ESM data of your own, simulate it — the code below generates a realistic dataset with known parameters, which also lets you check that your analysis recovers what you put in.

# ---- 00_simulate.R : an ESM dataset with known ground truth ---------- set.seed(2026) n_id <- 90; n_beep <- 35 truth <- list(b_within = -0.30, b_between = -0.50, sd_int = 0.85, sd_slope = 0.15) esm <- expand_grid(id = 1:n_id, beep = 1:n_beep) |> mutate(day = ceiling(beep / 5)) |> group_by(id) |> mutate(trait_stress = rnorm(1, 0, 1), u_int = rnorm(1, 0, truth$sd_int), u_slope = rnorm(1, 0, truth$sd_slope), stress = 3 + trait_stress + rnorm(n(), 0, 1)) |> ungroup() |> group_by(id) |> mutate(stress_pm = mean(stress), stress_cw = stress - stress_pm) |> ungroup() |> mutate(mood = 5 + u_int + (truth$b_within + u_slope) * stress_cw + truth$b_between * (stress_pm - mean(stress_pm)) + rnorm(n(), 0, 0.9), mood = if_else(runif(n()) < 0.18, NA_real_, mood)) # realistic missingness # ---- 01_clean.R ---- esm <- esm |> arrange(id, day, beep) |> group_by(id, day) |> mutate(stress_lag = lag(stress_cw)) |> ungroup() |> filter(!is.na(mood)) # ---- 02_describe.R ---- compliance <- esm |> count(id) |> summarise(median = median(n), min = min(n), max = max(n)) icc_null <- performance::icc(lmer(mood ~ 1 + (1 | id), data = esm)) # ---- 03_model.R ---- m_main <- lmer(mood ~ stress_cw + stress_pm_c + (1 + stress_cw | id), data = esm) m_lag <- lmer(mood ~ stress_lag + stress_pm_c + (1 + stress_lag | id), data = esm) # Did the analysis recover the truth we built in? fixef(m_main)[c("stress_cw", "stress_pm_c")] unlist(truth[c("b_within", "b_between")])#> stress_cw stress_pm_c #> -0.3041 -0.5218 #> b_within b_between #> -0.30 -0.50 ← recovered within Monte Carlo error. The analysis works.
Why simulate first, always
Running your planned analysis on simulated data with known parameters — before the real data arrives — does four things at once: it proves the code runs, it proves the model is identified with your planned n, it produces the power curve for your preregistration (volume 8 day 80), and it means the day your data lands you press one button. This is the single habit that most distinguishes researchers who finish on time.
DeliverableFile
Numbered, rerunnable scriptsR/00_simulate.R … R/03_model.R
Descriptives and compliance tableoutput/tables/table1.docx
Mixed-model tableoutput/tables/table3_mixed.html
Within-person slopes figureoutput/figures/fig3_mixed.png, 85 mm, 600 dpi
Results sectionmanuscript.qmd — volume 9 makes the numbers come from the model
Do this now · 30 minutes
Run the simulation, the cleaning and the two models as four separate numbered scripts, then write the results paragraph from their output. Check the recovery line: if your estimates do not return the parameters you built in, the bug is in your code, not your theory.
Reference

Every random term you will need

The parenthesis is the whole language. Read it as "(what varies | what it varies across)".

TermMeaning
(1 | id)Each person has their own intercept. The baseline mixed model.
(1 + x | id)Intercept and the effect of x vary by person, and their correlation is estimated.
(1 + x || id)Same, but the correlation is fixed at zero. The first simplification when a model will not converge.
(0 + x | id)Slope varies, intercept does not. Rarely what you mean.
(1 | school/class)Nested: expands to (1|school) + (1|school:class).
(1 | subject) + (1 | item)Crossed: two independent populations sampled.
(1 | id) + (1 | id:day)Three levels: beeps within days within people.
(1 + cond | subject) + (1 | item)The standard psycholinguistics design: condition varies within subject, items sampled.
Reference

The order to do things in

1 · structure
Write down what a row is, what the clusters are, and how many of each. Everything follows from this.
2 · empty model
Fit y ~ 1 + (1 | id), report the ICC. Now you know how much clustering there is.
3 · centre
Split every time-varying predictor into person-mean and within-person deviation before any inference.
4 · fixed effects
Add the predictors your hypotheses name. Resist adding others.
5 · random slopes
Add slopes for effects that vary within cluster; compare with refit = FALSE.
6 · diagnose
check_model(), convergence warnings, singular fits, cluster counts per term.
7 · report
Table with variance components, marginal-effects figure, all eight checklist items.
Reference

Install these once

install.packages(c( "lme4", # lmer and glmer: the engine "lmerTest", # Satterthwaite p-values in summary() "afex", # mixed() with Kenward-Roger, and ANOVA-style tables "performance", # icc(), r2(), check_model(), check_overdispersion() "emmeans", # emtrends() for slopes per group "ggeffects", # ggpredict() for model-implied figures "sjPlot", # tab_model() publication tables "modelsummary",# flexible model tables to Word/LaTeX "glmmTMB", # negative binomial, zero-inflated, beta GLMMs "pbkrtest", # parametric bootstrap and Kenward-Roger "broom.mixed" # tidy() for mixed models ))
Checkpoint

Seven questions before volume 6

{{ quizCounter }}
{{ quizScore }}
{{ quizQ }}
{{ quizFb }}

Next: volume 6

Every model so far assumed your variables measure what you claim. Volume 6 tests that: item analysis, reliability beyond alpha, factor structure, confirmatory models, measurement invariance and full SEM — fourteen days on the question a reviewer asks first.

Start day 53 →