## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 4.2)
library(silentema)

## ----params-------------------------------------------------------------------
p <- default_params()
p$Phi
round(partial_cors(p$Psi), 2)
round(diag(stationary_cov(p$Phi, p$Psi)), 3)

## ----basic--------------------------------------------------------------------
sim <- simulate_ema(N = 50, n_prompts = 30, motifs = "M0", compliance = 0.75, seed = 1)
sim
head(sim$data)

## ----full---------------------------------------------------------------------
f <- fit_pairs(sim$data, sim$vars)
round(cbind(estimate = f$Phi[, 1], truth = p$Phi[, 1]), 3)

## ----motifs-------------------------------------------------------------------
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"))

## ----channels-----------------------------------------------------------------
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)
c(probe_share = mean(sim2$data$Z), answered_at_probes = mean(sim2$data$R[sim2$data$Z == 1]),
  context_rate = mean(sim2$data$C))

## ----seed---------------------------------------------------------------------
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)

## ----fromfit------------------------------------------------------------------
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)))
# the simulated data carry the self-censoring signature of the fitted model
silence_test(d, sim3$vars)$coef_R[1]

## ----bias---------------------------------------------------------------------
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)

## ----sensor-design------------------------------------------------------------
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))

## ----probe-design-------------------------------------------------------------
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))

## ----estimators---------------------------------------------------------------
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)