Simulating experience-sampling data and planning a design

The generating model

simulate_ema() generates data from a two-level VAR(1): person i has a mean vector mu_i drawn from a multivariate normal with mean mu and covariance Sigma_mu, and the momentary states follow

\[X_{it} = \mu_i + \Phi (X_{i,t-1} - \mu_i) + \varepsilon_{it}, \qquad \varepsilon_{it} \sim N(0, \Psi).\]

The population parameters default to default_params(): four variables (negative affect, positive affect, stress, fatigue) with autoregressive coefficients between .30 and .45, a few cross-lagged effects, a sparse contemporaneous network in the innovations, and stationary variances scaled to one.

p <- default_params()
p$Phi
#>          NegA  PosA Stress Fatigue
#> NegA     0.40  0.00   0.15     0.0
#> PosA    -0.10  0.35   0.00     0.0
#> Stress   0.15  0.00   0.45     0.1
#> Fatigue  0.00 -0.10   0.00     0.3
round(partial_cors(p$Psi), 2)
#>         NegA PosA Stress Fatigue
#> NegA     1.0 -0.3    0.3     0.0
#> PosA    -0.3  1.0    0.0     0.0
#> Stress   0.3  0.0    1.0     0.2
#> Fatigue  0.0  0.0    0.2     1.0
round(diag(stationary_cov(p$Phi, p$Psi)), 3)
#> [1] 1 1 1 1

Whole prompts are then deleted according to one or more missingness motifs. The response model is a probit,

\[P(R_{it} = 1) = \Phi_N(\alpha_0 + a_i + \text{motif terms}),\]

and the intercept alpha₀ is calibrated numerically so that the realized response rate over the recorded prompts equals compliance. This calibration is what makes cells with different mechanisms comparable at the same response rate.

sim <- simulate_ema(N = 50, n_prompts = 30, motifs = "M0", compliance = 0.75, seed = 1)
sim
#> Simulated EMA data: 50 persons x 30 prompts, 4 variables
#>   motifs: M0   realized response rate: 0.75
head(sim$data)
#>   id time R     NegA     PosA    Stress  Fatigue
#> 1  1    1 1 2.124971 4.754357 2.5806521 4.874073
#> 2  1    2 1 2.616620 4.103303 1.9819488 4.480194
#> 3  1    3 1 1.365294 4.011420 1.3698565 2.887101
#> 4  1    4 1 1.775906 5.339526 2.1860586 2.831509
#> 5  1    5 1 1.791277 4.532241 1.7700848 5.390225
#> 6  1    6 1 1.150537 4.339816 0.8237349 3.062805

The result holds the observed data (data, with NA at skipped prompts), the complete states before deletion (full), the true person means (mu_i), the calibrated intercept (alpha0) and the settings of the response model. Together with the population parameters in default_params(), this makes it easy to check any estimator against the truth:

f <- fit_pairs(sim$data, sim$vars)
round(cbind(estimate = f$Phi[, 1], truth = p$Phi[, 1]), 3)
#>         estimate truth
#> NegA       0.409  0.40
#> PosA      -0.112 -0.10
#> Stress     0.136  0.15
#> Fatigue    0.039  0.00

Every motif and its parameters

The table lists the argument that controls each motif and the term it adds to the probit index (or to the state equation).

Motif	Argument(s)	Term
M1 lagged-state dependence	`gamma`	`gamma' X_{t-1}` (a scalar applies to the first variable)
M2 self-censoring	`delta`	`delta' X_t` (negative: high states are skipped)
M3 burden	`kappa_R`, `burden`	`kappa_R (R_{t-1} - compliance)` plus `burden * trend_t`, the trend running from -0.5 to 0.5 over burn-in and study (negative `burden` = declining compliance)
M4 person propensity	`rho_propensity`, `sd_propensity`	`a_i`, with correlation `rho_propensity` with the person’s mean on the first variable
M5 context	`p_context`, `rho_context`, `gamma_C`, `kappa_C`	binary `C_t` shifts the states by `gamma_C` and the index by `-kappa_C (C_t - p_context)`
M6 reactivity	`rho_react`	the state equation gains `rho_react (R_{t-1} - compliance)`

Both gamma and delta act on the raw (uncentered) state, so a self-censoring person with a high typical level also responds less often overall; this is deliberate, because it is how self-censoring behaves in practice.

rate_after <- function(motifs, ...) {
  s <- simulate_ema(N = 50, n_prompts = 30, motifs = motifs, compliance = 0.7, seed = 2, ...)
  fc <- fatigue_check(s$data)
  c(response_rate = round(mean(s$data$R), 3), persistence = round(fc$difference, 3))
}
rbind(M0 = rate_after("M0"), M1 = rate_after("M1"), M2 = rate_after("M2", delta = -1),
      M3 = rate_after("M3", kappa_R = 1), M4 = rate_after("M4"), M6 = rate_after("M6"))
#>    response_rate persistence
#> M0           0.7       0.044
#> M1           0.7       0.126
#> M2           0.7       0.230
#> M3           0.7       0.376
#> M4           0.7       0.196
#> M6           0.7       0.044

Every mechanism hits the target response rate; they differ in how the skips are arranged in time and in who skips.

Design features

Three optional channels correspond to the calibration designs of calibrate_delta():

sensor_cor adds an always-observed sensor S correlated with the first state variable (within person);
p_probe adds randomized probe prompts Z that force a response;
the context of M5 is always recorded in column C, so an analysis may treat it as observed (fit_pairs(covariates = "C"), fit_ipw()) or ignore it (latent context).

days adds a day column, which the estimators use to form pairs within days only when it is passed to them as day = "day".

sim2 <- simulate_ema(N = 30, n_prompts = 20, motifs = c("M2", "M5"), delta = -1, sensor_cor = 0.6,
                     p_probe = 0.1, rho_context = 0.5, days = 5, seed = 3)
names(sim2$data)
#>  [1] "id"      "time"    "day"     "R"       "NegA"    "PosA"    "Stress" 
#>  [8] "Fatigue" "C"       "S"       "Z"
c(probe_share = mean(sim2$data$Z), answered_at_probes = mean(sim2$data$R[sim2$data$Z == 1]),
  context_rate = mean(sim2$data$C))
#>        probe_share answered_at_probes       context_rate 
#>          0.1266667          1.0000000          0.2916667

Reproducibility

seed sets the random-number seed for the simulation and restores the caller’s random-number state afterwards, so a seeded call inside a larger script does not disturb the stream of that script:

set.seed(100); a <- runif(1)
set.seed(100); invisible(simulate_ema(N = 5, n_prompts = 5, seed = 7)); b <- runif(1)
identical(a, b)
#> [1] TRUE

The same convention holds for simulate_from_fit(), silence_test(se = "dayblock") and calibrate_delta().

Simulating from a fitted model

simulate_from_fit() generates data from a fitted tilt model: the fitted VAR(1), the fitted person means and, with propensity = "person", the fitted response intercepts, resampled jointly so that the association between a person’s level and propensity is preserved, and the probit self-censoring model at the fit’s sensitivity value. It is the engine of the post-skip calibration and a convenient parametric bootstrap.

sim3 <- simulate_ema(N = 60, n_prompts = 40, motifs = "M2", compliance = 0.7, delta = -1, seed = 4)
ft <- fit_tilt(sim3$data, sim3$vars, delta = -1)
d <- simulate_from_fit(ft, N = 60, n_prompts = 40, seed = 5)
c(response_rate = mean(d$R), persons = length(unique(d$id)))
#> response_rate       persons 
#>     0.7029167    60.0000000
# the simulated data carry the self-censoring signature of the fitted model
silence_test(d, sim3$vars)$coef_R[1]
#> [1] -0.488318

The post-skip contrast of a single simulated data set is noisy; calibrate_delta(method = "postskip") averages it over n_sim data sets at every grid value. An optional burden term (kappa_R) and a day structure (days) are available for the burden-aware calibration.

Planning a design: three small studies

The simulator makes design questions concrete. Each study below is deliberately small; a real planning exercise would use more replications and the sample size of the planned study.

How much does self-censoring bias the default analysis?

bias <- function(delta, reps = 20) {
  est <- vapply(seq_len(reps), function(r) {
    s <- simulate_ema(N = 60, n_prompts = 40, motifs = "M2", compliance = 0.7, delta = delta, seed = 100 + r)
    f <- fit_pairs(s$data, s$vars)
    c(f$Phi[1, 1], f$mu[1])
  }, numeric(2))
  c(bias_Phi11 = mean(est[1, ]) - p$Phi[1, 1], bias_mu1 = unname(mean(est[2, ]) - p$mu[1]))
}
round(rbind("delta = 0" = bias(0), "delta = -0.5" = bias(-0.5), "delta = -1" = bias(-1)), 3)
#>              bias_Phi11 bias_mu1
#> delta = 0         0.002    0.024
#> delta = -0.5     -0.018   -0.277
#> delta = -1       -0.057   -0.444

The autoregression is attenuated and the person mean is pulled down (people skip their high moments), both increasingly with the strength of self-censoring, with the mean the more damaged of the two. With 20 replications the Monte Carlo standard error of the bias in the autoregression is about .01, so small entries should be read as noise.

Which sensor is good enough to calibrate?

The precision of the sensor calibration depends on how strongly the sensor loads on the self-censored state. The loop compares the width of the calibration interval for two sensor correlations.

width <- function(sensor_cor, reps = 3) {
  w <- vapply(seq_len(reps), function(r) {
    s <- simulate_ema(N = 60, n_prompts = 40, motifs = "M2", compliance = 0.7, delta = -1,
                      sensor_cor = sensor_cor, seed = 200 + r)
    pr <- tilt_profile(s$data, s$vars, delta_grid = c(-2, -1.5, -1, -0.5, 0))
    cs <- calibrate_delta(pr, s$data, method = "sensor")
    diff(cs$interval)
  }, numeric(1))
  c(median_width = median(w, na.rm = TRUE), finite = mean(is.finite(w)))
}
rbind("r = 0.3" = width(0.3), "r = 0.6" = width(0.6))
#>         median_width    finite
#> r = 0.3    0.9136317 0.6666667
#> r = 0.6    0.4451380 1.0000000

A sensor with a modest loading gives wide or open intervals; a loading around .6 gives a usable calibration at this sample size.

How many probes?

probe_se <- function(p_probe, reps = 3) {
  se <- vapply(seq_len(reps), function(r) {
    s <- simulate_ema(N = 60, n_prompts = 40, motifs = "M2", compliance = 0.7, delta = -1,
                      p_probe = p_probe, seed = 300 + r)
    pr <- tilt_profile(s$data, s$vars, delta_grid = c(-2, -1.5, -1, -0.5, 0), probe = "Z")
    calibrate_delta(pr, s$data, method = "probe")$se
  }, numeric(1))
  c(median_se_of_contrast = median(se))
}
rbind("5% probes" = probe_se(0.05), "15% probes" = probe_se(0.15))
#>            median_se_of_contrast
#> 5% probes             0.10274056
#> 15% probes            0.06390391

The standard error of the probe contrast, which the calibration inverts, shrinks roughly with the square root of the number of probe prompts. Probes also recover the transition kernel directly (recoverability(dm_graph("M2", probe = TRUE))), but the probe-only sample is small, so the calibration route is usually the more informative use of them.

Comparing estimators across mechanisms

The last study reproduces, in miniature, the logic of the damage assessment in Yu (2026): the answered-pairs estimator and the MAR likelihood (fit_fiml()) are both fine under a recoverable mechanism and both biased under self-censoring.

compare <- function(motifs, ...) {
  s <- simulate_ema(N = 80, n_prompts = 40, motifs = motifs, compliance = 0.7, seed = 6, ...)
  c(pairs = fit_pairs(s$data, s$vars)$Phi[1, 1], fiml = fit_fiml(s$data, s$vars)$Phi[1, 1])
}
round(rbind(M1 = compare("M1"), M4 = compare("M4"), M2 = compare("M2", delta = -1),
            truth = c(p$Phi[1, 1], p$Phi[1, 1])), 3)
#>       pairs  fiml
#> M1    0.405 0.383
#> M4    0.391 0.397
#> M2    0.306 0.337
#> truth 0.400 0.400

One replication per cell cannot separate bias from sampling error; the pattern over many replications is in the accompanying article and in its archived simulation results.

References

Yu, H.-T. (2026). What skipped prompts hide: Detecting, diagnosing, and correcting informative nonresponse in ecological momentary assessment. Manuscript under review. Materials: https://osf.io/x6d2t/