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

## ----sim----------------------------------------------------------------------
sim <- simulate_ema(N = 80, n_prompts = 40, motifs = "M2", compliance = 0.7, delta = -1,
                    sensor_cor = 0.6, p_probe = 0.10, days = 5, seed = 21)
sim
truth <- default_params()$Phi

## ----tilt---------------------------------------------------------------------
f0 <- fit_tilt(sim$data, sim$vars, delta = 0, day = "day", probe = "Z")
f1 <- fit_tilt(sim$data, sim$vars, delta = -1, day = "day", probe = "Z")
rbind(untilted = f0$Phi[1, ], tilted = f1$Phi[1, ], truth = truth[1, ])

## ----means--------------------------------------------------------------------
rbind(untilted = f0$mu, tilted = f1$mu, truth = default_params()$mu)

## ----ess----------------------------------------------------------------------
c(pairs = f1$n_pairs, ess = round(f1$ess), max_weight = round(f1$max_weight, 2),
  iterations = f1$iterations, converged = f1$converged)

## ----profile------------------------------------------------------------------
prof <- tilt_profile(sim$data, sim$vars, delta_grid = seq(-2, 0.5, by = 0.5), day = "day", probe = "Z")
prof

## ----profile-parts------------------------------------------------------------
subset(prof$Phi, coef == "Stress<-NegA")[, c("delta", "estimate", "lower", "upper")]
subset(prof$pcor, edge == "NegA--Stress")
round(prof$mu, 2)

## ----plot-profile-------------------------------------------------------------
plot(prof, plausible = c(-1.5, -0.5), calibrated = -1,
     main = "Effects of NegA at t - 1 along the sensitivity grid")

## ----breakeven----------------------------------------------------------------
be <- break_even(prof, delta_max = 1.5)
be[be$coef %in% c("NegA<-NegA", "Stress<-NegA", "PosA<-NegA"), ]

## ----bounds-------------------------------------------------------------------
range(prof$mu[abs(prof$delta_grid + 1) <= 0.5, "NegA"])          # identified set of the mean over [-1.5, -0.5]
b <- bounds_support(sim$data, sim$vars, L = 1, U = 7)
summary(b$upper[b$variable == "NegA"] - b$lower[b$variable == "NegA"])   # per-person support-bound widths

## ----cal-sensor---------------------------------------------------------------
cs <- calibrate_delta(prof, sim$data, method = "sensor", day = "day")
cs
plot(cs, main = "Sensor calibration")

## ----cal-probe----------------------------------------------------------------
cp <- calibrate_delta(prof, sim$data, method = "probe", day = "day")
cp

## ----cal-postskip-------------------------------------------------------------
ck <- calibrate_delta(prof, sim$data, method = "postskip", day = "day", n_sim = 5, seed = 2)
ck

## ----cal-compare--------------------------------------------------------------
data.frame(method = c("sensor", "probe", "postskip"),
           delta_hat = c(cs$delta_hat, cp$delta_hat, ck$delta_hat),
           lower = c(cs$interval[1], cp$interval[1], ck$interval[1]),
           upper = c(cs$interval[2], cp$interval[2], ck$interval[2]))

## ----burden-------------------------------------------------------------------
sim23 <- simulate_ema(N = 60, n_prompts = 40, motifs = c("M2", "M3"), compliance = 0.7, delta = -1,
                      kappa_R = 1, days = 5, seed = 12)
prof23 <- tilt_profile(sim23$data, sim23$vars, delta_grid = c(-2, -1.5, -1, -0.5, 0), day = "day")
c_none <- calibrate_delta(prof23, sim23$data, method = "postskip", day = "day",
                          n_sim = 4, seed = 2, burden = "none")
c_fit <- calibrate_delta(prof23, sim23$data, method = "postskip", day = "day",
                         n_sim = 4, seed = 2, burden = "fit", kappa_grid = c(0, 0.5, 1, 1.5, 2))
c(M2_only = c_none$delta_hat, burden_fitted = c_fit$delta_hat)
c_fit$curve

## ----boot, eval = FALSE-------------------------------------------------------
#  calibrate_delta(prof, sim$data, method = "sensor", day = "day", boot = 200, seed = 3)

## ----declare------------------------------------------------------------------
g <- dm_graph("M2", sensor = TRUE, probe = TRUE)
st <- silence_test(sim$data, sim$vars, day = "day")
sg <- sensor_gap_test(sim$data, sim$vars, sensor = "S", day = "day")
decl <- missingness_declaration(g, silence = st, sensor_gap = sg, profile = prof,
                                calibration = list(cs, cp, ck), plausible = c(-1.5, -0.5))
cat(decl[grep("^## 5", decl):length(decl)], sep = "\n")