6  Data analysis

Analysis choices follow from the research question (Section 2.1). Here we start with the workhorse, the multilevel regression model, then we’ll consider temporal dependency, lagged and dynamic models, variability, non-normal outcomes, person-specific models, missing data, power, and causal inference.

6.1 Why multilevel models?

Prompts from the same person are more similar to each other than prompts from different people (the ICC for negative affect in our example dataset is 0.47). Ordinary regression treats all 6,294 responses as independent, which makes standard errors far too small for person-level effects. Averaging each person’s responses into one score avoids that but throws away all within-person information. Multilevel models (also called mixed-effects, or random-effects models) keep every observation while modelling the nesting:

  • they estimate within-person and between-person effects in one model;
  • they allow effects to differ across people (random slopes) and explain those differences with person-level variables;
  • they handle unequal numbers of observations per person without listwise deletion of people;
  • they give valid inference under the missing-at-random assumption for missing outcomes.

6.2 Multilevel regression, step by step

Considering the example dataset, the running question: When people are more stressed than usual, is their negative affect (NA) higher, and is this link stronger in people high in neuroticism?

\[ \begin{aligned} \text{Level 1:}\quad \text{NA}_{ti} &= \beta_{0i} + \beta_{1i}\left(\text{stress}_{ti} - \overline{\text{stress}}_{i}\right) + e_{ti} \\ \text{Level 2:}\quad \beta_{0i} &= \gamma_{00} + \gamma_{01}\,\overline{\text{stress}}_{i} + \gamma_{02}\,\text{neuro}_{i} + u_{0i} \\ \beta_{1i} &= \gamma_{10} + \gamma_{11}\,\text{neuro}_{i} + u_{1i} \end{aligned} \tag{6.1}\]

\(e_{ti} \sim N(0, \sigma^2)\) is moment-level residual variation. \((u_{0i}, u_{1i})\) are person-specific deviations from the average intercept and slope, assumed bivariate normal with variances \(\tau_{00}\), \(\tau_{11}\) and covariance \(\tau_{01}\). \(\gamma_{10}\) is the average within-person effect, \(\gamma_{01}\) the between-person effect, and \(\gamma_{11}\) the cross-level interaction.

Random intercepts and random slopes

The “multilevel” part of the model lives in the random-effects term, written (1 + stress_cw | id) in lme4.

A random intercept gives every person their own baseline level of the outcome. People differ in how negative they feel on average; the random intercept \(u_{0i}\) captures that spread (its variance is \(\tau_{00}\)), so the model does not pretend everyone shares one grand mean. It is also what accounts for the nesting: because repeated prompts from one person share that person’s intercept, the model treats them as non-independent and keeps the standard errors honest. You almost always want at least a random intercept for the grouping variable — it is the minimum that makes a model “multilevel”. In lme4 it is written (1 | id), and it answers the question how much do people differ in their average level?

A random slope goes further and gives every person their own effect of a within-person predictor. If stress reactivity varies — some people’s negative affect leaps when they are stressed, others’ barely moves — then the slope of stress_cw is itself a person-level quantity, \(u_{1i}\), with variance \(\tau_{11}\). In lme4 you add it inside the same bracket: (1 + stress_cw | id) means “let both the intercept and the stress_cw slope vary across people.” It answers how much do people differ in the effect? — often the more interesting question in ESM, and the one a cross-level moderator (here, neuroticism) tries to explain.

Two practical points follow:

  • Why you usually need the random slope. If a within-person effect genuinely varies across people and you omit its random slope, the model underestimates the uncertainty in the corresponding fixed effect and in any cross-level interaction, giving standard errors that are too small and p-values too optimistic. As a default, include a random slope for any within-person predictor your hypothesis concerns.
  • When one or both. Always start with the random intercept. Add random slopes for the within-person predictors you care about, budget permitting: each costs a variance, plus a covariance with the intercept (\(\tau_{01}\)) that is itself informative — for instance, whether people with higher baseline NA also react more strongly. Level-2 predictors that do not vary within a person (neuroticism, the person-mean of stress) cannot have a random slope; they enter as fixed effects only. When the data cannot support every random slope, simplify in a pre-specified order (drop the correlation first, then the least central slope) rather than by trial and error.

Writing and running the model in lmer

It helps to read the lmer formula as the equation above, assembled one term at a time. Start from the empty model and build up:

  • na ~ 1 + (1 | id) — the null model: a grand mean plus a random intercept (this is the model used for the ICC in Section 5.8).
  • na ~ stress_cw + (1 | id) — add the within-person predictor as a fixed effect: the average effect of being more stressed than one’s own usual.
  • na ~ stress_cw + (1 + stress_cw | id) — let that effect vary across people (the random slope).
  • na ~ stress_cw + stress_pm_gc + neuro + (1 + stress_cw | id) — add the between-person predictors: the person-mean of stress and the neuroticism trait.
  • na ~ stress_cw + stress_pm_gc + neuro + stress_cw:neuro + (1 + stress_cw | id) — add the cross-level interaction stress_cw:neuro, which asks whether the within-person stress slope depends on neuroticism.

The rule for reading the formula: everything inside the brackets and to the left of the | is a random effect (allowed to vary across the units on the right of the |, here id); everything outside the brackets is a fixed effect. stress_cw:neuro is a product term of a level-1 and a level-2 predictor, hence a cross-level interaction. lmerTest wraps lme4 to add p-values; REML = TRUE gives less biased variance estimates for the final model, and bobyqa is a robust optimiser. Putting the full formula together and fitting it, we also compare it with a three-level version that adds a day-level random intercept:

library(lmerTest)
m1 <- lmer(na ~ stress_cw + stress_pm_gc + neuro + stress_cw:neuro +
             (1 + stress_cw | id),
           data = esm, REML = TRUE,
           control = lmerControl(optimizer = "bobyqa"))
summary(m1)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: na ~ stress_cw + stress_pm_gc + neuro + stress_cw:neuro + (1 +  
    stress_cw | id)
   Data: esm
Control: lmerControl(optimizer = "bobyqa")

REML criterion at convergence: 45885.9

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.4516 -0.6841 -0.0265  0.6589  3.4473 

Random effects:
 Groups   Name        Variance Std.Dev. Corr 
 id       (Intercept) 51.734   7.193         
          stress_cw    1.558   1.248    0.02 
 Residual             79.668   8.926         
Number of obs: 6294, groups:  id, 100

Fixed effects:
                Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)      26.4952     0.7319 97.2413  36.203  < 2e-16 ***
stress_cw         2.7452     0.1571 96.0692  17.470  < 2e-16 ***
stress_pm_gc      4.8827     0.7329 97.5430   6.662 1.61e-09 ***
neuro             3.4292     0.7332 97.2347   4.677 9.40e-06 ***
stress_cw:neuro   1.0733     0.1575 98.5247   6.814 7.63e-10 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) strss_ strs__ neuro
stress_cw   0.019                     
strss_pm_gc 0.010  0.000              
neuro       0.098  0.002  0.097       
strss_cw:nr 0.002  0.104  0.000  0.018
m1b <- lmer(na ~ stress_cw + stress_pm_gc + neuro + stress_cw:neuro +
              (1 + stress_cw | id) + (1 | id:day),
            data = esm, REML = TRUE,
            control = lmerControl(optimizer = "bobyqa"))
anova(m1, m1b, refit = FALSE)
Data: esm
Models:
m1: na ~ stress_cw + stress_pm_gc + neuro + stress_cw:neuro + (1 + stress_cw | id)
m1b: na ~ stress_cw + stress_pm_gc + neuro + stress_cw:neuro + (1 + stress_cw | id) + (1 | id:day)
    npar   AIC   BIC logLik -2*log(L)  Chisq Df Pr(>Chisq)    
m1     9 45904 45965 -22943     45886                         
m1b   10 45754 45822 -22867     45734 151.88  1  < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Reading the output

  • Within-person effect (stress_cw = 2.75): at moments when a person with average neuroticism is 10 points more stressed than their own average, their NA is 2.75 points higher.
  • Between-person effect (stress_pm_gc = 4.88): people whose average stress is 10 points higher have average NA that is 4.88 points higher. The contextual effect, 4.88 − 2.75 = 2.14, is the extra between-person association beyond what within-person processes alone would imply.
  • Cross-level interaction (1.07): per SD of neuroticism, the within-person stress slope is about 1.07 points steeper.
  • Random slope SD (1.25): people differ in stress reactivity. For people with average neuroticism, roughly 95% of individual slopes lie within 2.75 ± 1.96 × 1.25, i.e. about 0.3 to 5.2. This spread is often more informative than the average.
  • Degrees of freedom (Satterthwaite, via lmerTest) are close to the number of participants, not the number of observations: for effects that vary across people, the person is the effective unit of replication.
  • Day-level variance: adding a random intercept for days (m1b) significantly improves fit, so negative affect also varies systematically between days within a person. A three-level model is appropriate here.

Model-building advice

  • Convergence warnings and singular fits usually mean the random-effects structure is too complex for the data, or variables are on very different scales. Rescale predictors, try another optimiser (lmerControl(optimizer = "bobyqa")), and simplify the random-effects structure in a pre-specified order (e.g., drop random-effect correlations before dropping slopes). Bayesian estimation with weakly informative priors (e.g., brms) is another route.
  • REML vs ML. Use REML for final estimates of variance components. Compare models with different fixed effects using ML (anova() refits automatically) and models with different random effects on REML fits.
  • Time. Consider study day, time of day, and weekday/weekend as covariates. A trend shared by two variables can produce a spurious within-person association.

Why report unstandardised effects

Report coefficients in their original units (in the example, points of negative affect per 10 points of stress) rather than standardised (\(\beta\)) coefficients. There are three reasons.

First, in a multilevel model it is unclear which standard deviation to divide by. A predictor like stress_cw has a within-person SD, a between-person SD, and a total SD, and they differ; standardising by the total SD blends the very levels you worked to separate (Section 5.7), and there is no single agreed standardised effect size for within-person effects.

Second, standardised coefficients are hostage to the sample’s variances. Two studies with different compliance, sampling ranges, or populations will have different SDs, so the same underlying process yields different \(\beta\)s — which makes standardised effects hard to compare across ESM studies, the opposite of what standardisation is meant to achieve.

Third, unstandardised coefficients are directly interpretable and reusable: a reader can reason about “3 points of NA per 10 points of stress” on the actual scale, and a future power analysis can take it as an input. To convey magnitude, pair the unstandardised effects with variance explained at each level (Rights & Sterba, 2019, in the r2mlm package) and with the random-slope SD, which shows how much the effect varies across people.

Fitting the same model with brms

This guide uses lme4 throughout because it is fast, familiar, and enough for most ESM questions. It is worth knowing the Bayesian alternative, brms, which fits the same models with a nearly identical formula but estimates them with Stan. It is slower, but it comes into its own for the harder cases in this chapter — non-normal outcomes (Section 6.8), location-scale (variability) models (Section 6.7), and small samples or complex random-effects structures where maximum likelihood struggles to converge — and it returns full posterior distributions and credible intervals for every quantity.

library(brms)
m1_brms <- brm(
  na ~ stress_cw + stress_pm_gc + neuro + stress_cw:neuro +
    (1 + stress_cw | id),
  data = esm, chains = 4, cores = 4, seed = 2026,
  prior = prior(normal(0, 5), class = "b"))   # weakly informative priors on fixed effects
summary(m1_brms)

The formula is exactly the model fitted above; prior() adds weakly informative priors on the fixed effects, and summary() reports posterior means with 95% credible intervals in place of p-values. For latent (rather than observed) person-mean centering, and for the dynamic models later in this chapter, Bayesian estimation via brms or DSEM (Section 6.5) is often the more natural choice.

6.3 Autocorrelated residuals

Momentary states carry over, so residuals of close prompts are correlated. A standard multilevel model assumes residuals are independent within people after accounting for random effects. Ignoring autocorrelation mostly affects standard errors of level-1 effects, and it matters more with closely spaced prompts. The nlme package can model a first-order autoregressive residual structure:

library(nlme)
m2 <- lme(na ~ stress_cw + stress_pm_gc + neuro + stress_cw:neuro,
          random = ~ 1 + stress_cw | id,
          correlation = corAR1(form = ~ beep | id/day),
          data = esm, na.action = na.omit)
summary(m2)$tTable
                    Value Std.Error   DF   t-value       p-value
(Intercept)     26.490898 0.7341754 6192 36.082521 6.166463e-259
stress_cw        2.718165 0.1532864 6192 17.732592  1.145251e-68
stress_pm_gc     4.856963 0.7344632   97  6.612942  2.058318e-09
neuro            3.411521 0.7355716   97  4.637918  1.099766e-05
stress_cw:neuro  1.050667 0.1537188 6192  6.834996  8.985605e-12
m2$modelStruct$corStruct
Correlation structure of class corARMA representing
    Phi1 
0.258382

corAR1 assumes equally spaced integer occasions; because the beep number records each prompt’s scheduled position, a missed prompt correctly counts as a gap of two lags. For unequal spacing, corCAR1(form = ~ hours | id/day) models correlation as a function of continuous time. Note that nlme reports level-1 degrees of freedom based on observations, which is liberal for cross-level terms.

6.4 Lagged (temporal) effects

Lagged models predict a variable at prompt t from variables at prompt t − 1. Including the lagged outcome as a predictor turns the model into a question about change: does stress now predict NA beyond what NA at the previous prompt already predicts? The coefficient of the lagged outcome is the autoregressive effect, often interpreted as inertia.

m3 <- lmer(na ~ na_lag_cw + stress_cw + stress_pm_gc +
             (1 + na_lag_cw + stress_cw | id),
           data = esm, REML = TRUE,
           control = lmerControl(optimizer = "bobyqa"))
summary(m3)$coefficients
               Estimate Std. Error        df   t value     Pr(>|t|)
(Intercept)  26.2864819 0.81269083  98.66756 32.344996 2.530129e-54
na_lag_cw     0.2177502 0.02019549 103.86254 10.782119 1.224700e-18
stress_cw     2.5750677 0.19815234  91.21513 12.995394 1.809546e-22
stress_pm_gc  4.3026390 0.79082950 100.41162  5.440666 3.756581e-07

Things to keep in mind:

  • Cross-lagged effects. To test whether stress predicts later NA, use stress_lag (person-mean centred) as a predictor while controlling for na_lag. Also test the reverse direction.
  • Nickell’s bias. With observed person-mean centering, autoregressive estimates are biased toward zero, and the bias grows as the number of observations per person shrinks. Latent centering in DSEM avoids this.
  • Measurement error attenuates autoregressive effects and can create spurious cross-lagged effects (Schuurman & Hamaker, 2019). Multi-item measures or models with a measurement-error term help.
  • The lag is a time interval. A lag-1 effect with prompts every two hours is a different quantity from a lag-1 effect with daily diaries. Report the average interval, and see Section 6.12 on continuous-time models.

6.5 Dynamic structural equation modelling

Dynamic SEM (DSEM, Asparouhov et al., 2018) combines multilevel modelling, time-series analysis and SEM, using Bayesian estimation. For ESM it has several advantages over lagged lmer models:

  • latent person-mean centering, which removes Nickell’s bias and the bias in between-person effects;
  • random autoregressive and cross-lagged effects and random residual variances, so individual differences in inertia, reactivity and variability can be estimated and predicted in one model;
  • missing data and time handled within the model, including placing observations on a time grid to approximate unequal intervals;
  • measurement models: latent variables at each level, and measurement error in single indicators.

\[ \begin{aligned} \text{Within:}\quad \text{NA}^{w}_{ti} &= \phi_{i}\,\text{NA}^{w}_{t-1,i} + \beta_{i}\,\text{stress}^{w}_{ti} + \zeta_{ti}, \qquad \zeta_{ti} \sim N\!\left(0, e^{\omega_i}\right) \\ \text{Between:}\quad (\mu_{i}, \phi_{i}, \beta_{i}, \omega_{i}) &\sim \text{MVN}, \text{ each optionally regressed on } \text{neuro}_{i} \end{aligned} \tag{6.2}\]

\(\text{NA}^{w}\) and \(\text{stress}^{w}\) are deviations from each person’s latent mean \(\mu_i\). \(\phi_i\) is person-specific inertia, \(\beta_i\) person-specific stress reactivity, and \(\omega_i\) the log innovation variance (moment-to-moment unpredictability).

Mplus
VARIABLE:  NAMES = id day beep na stress neuro;
           USEVARIABLES = na stress neuro;
           CLUSTER = id;
           BETWEEN = neuro;
           LAGGED = na(1);
           TINTERVAL = hour(1);    ! places observations on a time grid
ANALYSIS:  TYPE = TWOLEVEL RANDOM;
           ESTIMATOR = BAYES;
           BITERATIONS = (5000);
           PROCESSORS = 2;
MODEL:     %WITHIN%
           phi | na ON na&1;       ! random inertia
           s   | na ON stress;     ! random stress reactivity
           logv | na;              ! random residual variance
           %BETWEEN%
           na phi s logv ON neuro;
           na phi s logv WITH na phi s logv;

McNeish & Hamaker (2020) give an accessible primer with annotated Mplus code. In R, mlts (multilevel latent time series models via Stan), brms (flexible Bayesian multilevel models, with observed centering), dynr, and ctsem (continuous-time models, Driver et al., 2017) cover much of the same ground. Blimp also supports latent centering and lagged effects with missing data.

6.6 Multilevel VAR and network models

Vector autoregressive (VAR) models extend lagged models to several variables at once: each variable at t is predicted by all variables at t − 1. The multilevel VAR in the mlVAR package (Epskamp et al., 2018) estimates three networks: a temporal network (lagged effects), a contemporaneous network (partial correlations among residuals at the same prompt), and a between-person network (partial correlations among person means).

library(mlVAR)
fit_var <- mlVAR(esm,
                 vars     = c("na", "stress"),
                 idvar    = "id",
                 dayvar   = "day",     # no lags across days
                 beepvar  = "beep",    # missing beeps break lag pairs
                 lags     = 1,
                 temporal = "correlated")
plot(fit_var, "temporal")
plot(fit_var, "contemporaneous")
WarningInterpreting networks

Network models inherit every assumption of lagged regression (stationarity, equal intervals, no unmeasured confounders, reliable measures) and add a large number of parameters. Contemporaneous edges have no direction. Centrality indices are often unstable and their meaning for intervention targets is debated. Treat networks as exploratory descriptions unless the design supports stronger claims.

6.7 Modelling variability

Sometimes the question is not about the level of a variable but about how much it fluctuates: affective instability in borderline personality disorder, for instance. Instead of computing iSD or MSSD per person and entering them into a second-stage regression (which ignores their estimation error), model variability directly. The mixed-effects location-scale model (Hedeker et al., 2008) gives each person their own residual variance, which can be predicted by person-level variables.

library(brms)
fit_ls <- brm(
  bf(na ~ stress_cw + stress_pm_gc + (1 + stress_cw |p| id),
     sigma ~ neuro + (1 |p| id)),   # person-specific residual SD, predicted by neuroticism
  data = esm, chains = 4, cores = 4, seed = 2026)
# |p| lets the location and scale random effects correlate

6.8 Non-normal outcomes

Outcome Model R
Binary (alone yes/no, smoked since last prompt) Multilevel logistic regression lme4::glmer(family = binomial)
Ordinal (1–7 Likert with few used categories) Cumulative link mixed model ordinal::clmm, brms (cumulative)
Counts (cigarettes, drinks) Poisson or negative binomial mixed model glmmTMB, brms
Many zeros (symptoms usually absent) Zero-inflated or hurdle mixed model glmmTMB (ziformula), brms (hurdle_*)
Bounded slider with piles at 0 or 100 Beta, ordered beta, or censored models glmmTMB, brms

In non-linear mixed models, fixed effects are conditional on random effects (effects for a “typical” person), not population-averaged. Plot predicted values on the response scale to communicate results.

6.9 Person-specific (idiographic) models

A multilevel model describes an average process and the spread around it. Whether the average describes any individual well depends on ergodicity: group-level and individual-level structures coincide only if every person follows the same stationary process, which is rarely true in psychology (Molenaar, 2004). Options for person-specific analysis:

  • N = 1 time-series models (e.g., a VAR per person with the vars package). They need many observations per person, roughly equal spacing, and stationarity or explicit detrending.
  • Random effects as person-specific estimates. Best linear unbiased predictions (BLUPs; ranef() in lme4) borrow strength from the group and are more stable than separate per-person regressions.
  • GIMME (group iterative multiple model estimation, R package gimme) (Gates & Molenaar, 2012) estimates person-specific models that share group-level paths where the data support them, and can find subgroups with similar dynamics.
  • Clustering dynamics: mixture models or clustering on person-specific parameters to find subtypes.

6.10 Missing data

  • Missing outcomes: likelihood-based multilevel models use all observed responses and are valid if data are missing at random (MAR) given the variables in the model.
  • Missing predictors, including lagged predictors, cause rows to be dropped. A missed prompt removes up to two lag pairs (Section 5.6). Multilevel multiple imputation (mice with 2l.* methods, jomo with mitml, or Blimp) or Bayesian models that treat missing predictors as parameters (DSEM) retain them.
  • Make MAR more plausible by including variables that predict missingness (time of day, study day, weekend, previous-prompt states) as covariates or auxiliary variables.
  • Not-at-random missingness (people skip prompts because they feel bad at that moment) cannot be ruled out from the data. Sensitivity analyses, such as comparing results between high- and low-compliance participants or pattern-mixture approaches, show how much conclusions depend on the assumption.

6.11 Sample size and power

Sample size is a design decision, set before data collection; Section 2.5 covers which sample size (people versus prompts per person) matters for which kind of question. This section is the technical counterpart: how to actually estimate power for the multilevel models above. Beyond the two sample sizes, power here depends on the ICC, the size of random-slope variances, autocorrelation, and expected compliance — quantities you can only pin down once the model is specified.

Closed-form formulas rarely fit ESM designs, so simulation is the standard approach. The PowerAnalysisIL Shiny app (Lafit et al., 2021) performs simulation-based power analysis for common multilevel models with autocorrelation. For custom models, simulate data under plausible parameters (pilot data or published estimates), fit the planned model, and repeat:

# simulate_esm() is in R/simulate_esm.R; prepare_esm() stands for the
# Chapter 5 processing steps wrapped in a function of your own
power_for <- function(n_id, n_day, n_beep, n_sims = 500) {
  hits <- replicate(n_sims, {
    d   <- simulate_esm(n_id, n_day, n_beep, seed = NULL) |>   # new data each time
           prepare_esm()
    fit <- lmerTest::lmer(na ~ stress_cw * neuro + stress_pm_gc +
                            (1 + stress_cw | id), data = d)
    coef(summary(fit))["stress_cw:neuro", "Pr(>|t|)"] < .05
  })
  mean(hits)
}
designs <- expand.grid(n_id = c(60, 100, 150), n_day = c(7, 14), n_beep = 6)
designs$power <- mapply(power_for, designs$n_id, designs$n_day, designs$n_beep)

6.12 From temporal precedence to causal claims

Lagged associations are tempting to read causally. Threats to that reading:

  • Time-varying confounders. Person-mean centering removes stable differences between people, but not a third variable that changes within a person (a bad night’s sleep increasing both stress and NA the next day).
  • The interval problem. Lagged effects depend on the time between measurements. An effect that peaks after 30 minutes may be invisible with a two-hour interval, and effects at different intervals can even differ in sign. Continuous-time models (ctsem) estimate the underlying process and translate it to any interval.
  • Between-person confounding in panel models. Traditional cross-lagged panel models mix stable between-person differences with within-person dynamics; the random-intercept cross-lagged panel model separates them (Hamaker et al., 2015). Multilevel ILD models achieve the same separation through centering.
  • Measurement error and reverse effects, as discussed in Section 6.4.
  • Granger causality is predictive, not causal. “X improves prediction of later Y” is compatible with many causal structures.

The strongest designs manipulate the putative cause in daily life: micro-randomised trials, or randomised momentary interventions (Section 2.6). Where manipulation is impossible, state assumptions explicitly and consider causal inference frameworks for time-varying exposures.