Cookbook: common tasks

Short recipes for common tasks. Every recipe uses the simulated data below (continuous treatment X, mediator M, outcome Y, covariate C). Chunks are not evaluated when the site is built; the code was run against version 0.4.1.

library(probmed)
set.seed(123)
n <- 300
d <- data.frame(X = rnorm(n), C = rnorm(n))
d$M <- 0.5 * d$X + 0.3 * d$C + rnorm(n)
d$Y <- 0.4 * d$M + 0.2 * d$X + 0.2 * d$C + rnorm(n)

Point estimate only

fit <- pmed(Y ~ X + M + C, formula_m = M ~ X + C, data = d,
            treatment = "X", mediator = "M", method = "plugin", seed = 1)
fit@estimate

"plugin" returns the point estimate with no interval. Use it for quick checks.

Interval by method

All methods share one call; only method changes.

args <- list(Y ~ X + M + C, formula_m = M ~ X + C, data = d,
             treatment = "X", mediator = "M", seed = 1)
pb   <- do.call(pmed, c(args, method = "parametric_bootstrap", n_boot = 1000))
npb  <- do.call(pmed, c(args, method = "nonparametric_bootstrap", n_boot = 1000))
mbco <- do.call(pmed, c(args, method = "mbco"))

See Choosing a confidence interval for when to use which.

Read the result

mbco@estimate; mbco@ci_lower; mbco@ci_upper   # P_med and its interval
mbco@ie_estimate                               # indirect effect a*b
mbco@diagnostics$pmed                          # interval status, near_double_null flag

For method = "mbco", check that the printed interval status reads ok. A missing bound means a verified likelihood-ratio endpoint was not found; the method reports that instead of a guess.

Binary mediator or outcome with covariates

Use the model-based bootstrap, which refits both equations on simulated data.

pmed(Y_bin ~ X + M + C, formula_m = M ~ X + C, data = d_bin,
     treatment = "X", mediator = "M",
     method = "model_bootstrap", n_boot = 500, seed = 1)

A rare outcome (about 5%) at small N often returns bootstrap_incomplete instead of an interval. That is the method declining to report, not an error. Poisson, ordinal and other families, interactions and transformed treatment or mediator terms are rejected.

From a fitted lavaan or mediation object

extract_mediation() returns a MediationData object, and pmed() has a method for it. See lavaan integration and mediation integration.

Two mediators

For parallel mediators use the joint P_med; see Parallel mediation.

Tipping point for unmeasured M-Y confounding (mean-based shares)

pmed_sensitivity(indirect = 0.30, total = 0.50, threshold = 0, evalue = TRUE)

More in Sensitivity analysis.

Weak identification check

ward_residual(data.frame(A = as.integer(d$X > 0), M = d$M, Y = d$Y, C = d$C))

ward_residual() expects a binary 0/1 treatment A and covariates named in covars (default "C"). It reports gates for weak identification; see Gauge residual.

Missing values

pmed() stops if any variable named in either formula has a missing value, before fitting anything. Drop incomplete rows first so both models use the same sample.

d_complete <- d[stats::complete.cases(d[, c("X", "M", "Y", "C")]), ]

Columns that neither formula uses may contain missing values. Since 0.4.1 this check runs on the data you pass, so it also holds with medfit 0.10 and later.

Reproducibility

Pass seed to every call. Plug-in calls honor it, MBCO is deterministic, and bootstrap methods draw from seed.