| Title: | Probabilistic Effect Sizes for Causal Mediation Analysis |
|---|---|
| Description: | Compute P_med, a scale-free probabilistic effect size for causal mediation analysis. |
| Authors: | Davood Tofighi [aut, cre, cph] (ORCID: <https://orcid.org/0000-0001-8523-7776>) |
| Maintainer: | Davood Tofighi <[email protected]> |
| License: | MIT + file LICENSE |
| Version: | 0.4.1 |
| Built: | 2026-10-31 07:47:57 UTC |
| Source: | https://github.com/Data-Wise/probmed |
Generic function to extract mediation structure from fitted models. This function is provided by the medfit package and re-exported here for convenience.
object |
Fitted model object |
... |
Additional arguments passed to methods |
A medfit::MediationData object
extract_mediation for full documentation
Estimates g_l = NIE_l / TE_l on coarse-grained data by plug-in, with a
cluster-robust standard error from the efficient influence function.
The plug-in is already efficient in the linear-Gaussian case, so no
cross-fitting is performed.
g_hat_eif(Z, cluster, m = 1, ci_level = 0.95)g_hat_eif(Z, cluster, m = 1, ci_level = 0.95)
Z |
Numeric matrix with columns |
cluster |
Cluster labels for the rows of |
m |
Numeric: the aggregation scale |
ci_level |
Numeric: confidence level for the Wald interval. |
An RgFlowResult.
set.seed(1) d <- data.frame( A = rnorm(400), M = rnorm(400), Y = rnorm(400), sch = rep(1:40, each = 10) ) cg <- rg_coarse_grain(d, m = 2, cluster = "sch") fit <- g_hat_eif(cg$Z, cg$cluster, m = 2) fit@gset.seed(1) d <- data.frame( A = rnorm(400), M = rnorm(400), Y = rnorm(400), sch = rep(1:40, each = 10) ) cg <- rg_coarse_grain(d, m = 2, cluster = "sch") fit <- g_hat_eif(cg$Z, cg$cluster, m = 2) fit@g
S7 class for the gauge-calibrated proportion mediated. Reports the
interventional proportion mediated P_med = IIE / OE alongside the
gauge residual W = R / OE, where R = OE - IDE - IIE is the
non-decomposability (treatment-by-mediator interaction) term. A W
significantly different from zero signals that the additive split of the
overall effect fails and the single-number P_med is unreliable.
GaugePmedResult( p_med = integer(0), p_med_ci = integer(0), p_med_fieller = numeric(0), fieller_type = NA_character_, W = integer(0), W_ci = integer(0), W_se = integer(0), W_p = integer(0), W_ci_wald = numeric(0), weak_id = NA, weak_id_ratio = NA_real_, oe_snr = NA_real_, oe_regular = NA, weak_id_ratio_threshold = 3, oe_snr_threshold = 2, OE = integer(0), IDE = integer(0), IIE = integer(0), R = integer(0), theta = integer(0), method = character(0), n = integer(0), ci_level = integer(0), se_method = "analytic", reps = 1L, call = NULL )GaugePmedResult( p_med = integer(0), p_med_ci = integer(0), p_med_fieller = numeric(0), fieller_type = NA_character_, W = integer(0), W_ci = integer(0), W_se = integer(0), W_p = integer(0), W_ci_wald = numeric(0), weak_id = NA, weak_id_ratio = NA_real_, oe_snr = NA_real_, oe_regular = NA, weak_id_ratio_threshold = 3, oe_snr_threshold = 2, OE = integer(0), IDE = integer(0), IIE = integer(0), R = integer(0), theta = integer(0), method = character(0), n = integer(0), ci_level = integer(0), se_method = "analytic", reps = 1L, call = NULL )
p_med |
Numeric: interventional proportion mediated, |
p_med_ci |
Numeric length-2: confidence interval for |
p_med_fieller |
Numeric: Fieller confidence set for the ratio |
fieller_type |
Character: type of Fieller set returned (e.g. bounded,
unbounded, or empty); |
W |
Numeric: gauge residual |
W_ci |
Numeric length-2: confidence interval for |
W_se |
Numeric: standard error of |
W_p |
Numeric: two-sided p-value for |
W_ci_wald |
Numeric length-2: symmetric Wald interval for |
weak_id |
Logical: What it measures (2026-08-22, corrected corner weights;
History. The flag shipped with its 3x threshold "calibrated above the
~2x width gap of the anti-conservative Wald interval" and was validated on a
48,000-rep grid as a "narrowness proxy" whose coverage separation was
explained by Wald width ( |
weak_id_ratio |
Numeric: ratio of the percentile |
oe_snr |
Numeric: signal-to-noise of the denominator, |
oe_regular |
Logical: What |
weak_id_ratio_threshold, oe_snr_threshold
|
Numeric: the threshold values
actually used for the |
OE, IDE, IIE, R
|
Numeric: overall, interventional direct, interventional indirect effects and the remainder. |
theta |
Numeric length-4: corner means |
method |
Character: estimation method. |
n |
Integer: sample size. |
ci_level |
Numeric: confidence level. |
se_method |
Character: |
reps |
Integer: number of repeated cross-fitting fold draws averaged for the point estimate. |
call |
Call: original call. |
Estimate the incremental mediated elasticity P_med^delta(delta) – the
derivative-scale proportion mediated as a function of the treatment-tilt factor
delta. The estimator is cross-fitted and built on the same triply-robust
corner pseudo-outcomes as ward_residual(); standard errors use the
ratio-identity influence function for med / tot at each delta.
incr_pmed( object, deltas = c(1/3, 1/2, 1, 2, 3), covars = "C", K = 5L, ci_level = 0.95, seed = 1L, ... )incr_pmed( object, deltas = c(1/3, 1/2, 1, 2, 3), covars = "C", K = 5L, ci_level = 0.95, seed = 1L, ... )
object |
A data frame with columns |
deltas |
Numeric: tilt factors at which to evaluate the curve. |
covars |
Character: covariate column names (default |
K |
Integer: number of cross-fitting folds (default |
ci_level |
Numeric: confidence level (default |
seed |
Integer: random seed for fold assignment (default |
... |
Unused. |
At each delta the tilted assignment probability is
q = delta * g / (delta * g + 1 - g) with derivative
q' = g (1 - g) / (delta * g + 1 - g)^2, where g = P(A = 1 | C). The
direct and mediated elasticities are weighted contrasts of the four corner
pseudo-outcomes by q and q'; by Theorem 1 (multivariate chain rule) they
sum to the total elasticity exactly, so P_med^delta = med / (dir + med).
Inference. Point estimates use the tilt-derivative weight
q' = g(1 - g) / (delta g + 1 - g)^2 (Kennedy 2019, Corollary 2, Term I).
Standard errors use the full efficient influence function for the ratio
P_med = med / tot, which adds the g-score correction
(dq/dg) * (A - g(C)) * (dir * E[a_med | C] - med * E[a_dir | C]) / tot^2,
where a_med = q' * gamma_med and a_dir = q' * gamma_dir are the same
tilt-weighted corner contrasts used in the point estimate (the g-score
term must carry the q' weight to match the units of dir/med/tot)
and E[. | C] is a linear projection on the covariates, with
dq/dg = delta / (delta g + 1 - g)^2 (Term II). Term II is mean-zero (by
E[A - g(C) | C] = 0), so it does not shift the point estimate but
restores Neyman orthogonality w.r.t. the propensity score, ensuring the CI
is consistent under nonparametric estimation of g. The projection
matters: the per-row a_med/a_dir carry the corner EIF's
inverse-probability noise, and multiplying that by A - g(C) inflated the
se by 1.40x at delta = 2 (fixed 2026-08-22; se/empSD now 0.94-0.97 across
delta).
Scope (exposure type). The incremental intervention is defined only for a
binary treatment A: the odds-ratio tilt q = delta g / (delta g + 1 - g)
requires g = P(A = 1 | C) to be a Bernoulli probability. Continuous or
multi-valued exposures are a distinct object requiring a modified treatment
policy (shift) intervention with its own estimand and efficient influence
function; they are deliberately out of scope here (deferred to future work –
a separate MTP-mediation paper).
An IncrPmedResult object whose @curve slot holds one row per
delta with the direct, mediated, and total elasticities, the proportion
mediated Pmed, and its standard error and Wald interval.
set.seed(1) n <- 1500 C <- rnorm(n) A <- rbinom(n, 1, plogis(-0.2 + 0.8 * C)) M <- 0.6 * A + 0.4 * C + rnorm(n) Y <- 0.5 * A + 0.7 * M + 0.3 * C + rnorm(n) d <- data.frame(A, M, Y, C) incr_pmed(d, deltas = c(0.5, 1, 2))set.seed(1) n <- 1500 C <- rnorm(n) A <- rbinom(n, 1, plogis(-0.2 + 0.8 * C)) M <- 0.6 * A + 0.4 * C + rnorm(n) Y <- 0.5 * A + 0.7 * M + 0.3 * C + rnorm(n) d <- data.frame(A, M, Y, C) incr_pmed(d, deltas = c(0.5, 1, 2))
Per-delta M–Y unmeasured-confounding sensitivity for an IncrPmedResult.
For each tilt factor delta, reports the additive numerator bias that zeroes
the incremental mediated share (tipping) and, optionally, the bias that
drives P_med^delta to a user threshold. Reuses pmed_sensitivity();
the additive-offset model is exact for the Paper 2 (mean-based) share.
incr_sensitivity(object, threshold = 0, ...)incr_sensitivity(object, threshold = 0, ...)
object |
An IncrPmedResult from |
threshold |
Numeric or |
... |
Unused. |
A data frame: one row per delta with columns delta, med,
tot, Pmed, tipping, tipping_threshold.
set.seed(1) n <- 800; C <- rnorm(n) A <- rbinom(n, 1, plogis(-0.2 + 0.8 * C)) M <- 0.6 * A + 0.4 * C + rnorm(n) Y <- 0.5 * A + 0.7 * M + 0.3 * C + rnorm(n) fit <- incr_pmed(data.frame(A, M, Y, C), deltas = c(0.5, 1, 2)) incr_sensitivity(fit, threshold = 0)set.seed(1) n <- 800; C <- rnorm(n) A <- rbinom(n, 1, plogis(-0.2 + 0.8 * C)) M <- 0.6 * A + 0.4 * C + rnorm(n) Y <- 0.5 * A + 0.7 * M + 0.3 * C + rnorm(n) fit <- incr_pmed(data.frame(A, M, Y, C), deltas = c(0.5, 1, 2)) incr_sensitivity(fit, threshold = 0)
S7 class for the incremental mediated elasticity P_med^delta(delta)
(Paper 2 of the P_med-for-modern-estimands program). Rather than a single
number, the estimand is a curve: the share of the marginal (derivative-scale)
total effect that flows through the mediator when the treatment-assignment odds
are tilted by a factor delta.
Writing Theta(delta_1, delta_2) for the mean outcome under a delta_1-tilt of
the direct pathway and a delta_2-tilt of the mediated pathway, the elasticities
are theta'_dir = d/d delta_1 Theta, theta'_med = d/d delta_2 Theta, and
theta'_tot = d/d delta Theta(delta, delta). The multivariate chain rule gives
the exact decomposition theta'_tot = theta'_dir + theta'_med (remainder
R(delta) = 0 for every delta), and P_med^delta(delta) = theta'_med / theta'_tot.
Under a linear-Gaussian model with no A * M interaction the curve is flat and
equals the classical proportion mediated; it bends with delta only when an
interaction is present.
IncrPmedResult( curve = data.frame(), deltas = integer(0), method = character(0), n = integer(0), ci_level = integer(0), call = NULL )IncrPmedResult( curve = data.frame(), deltas = integer(0), method = character(0), n = integer(0), ci_level = integer(0), call = NULL )
curve |
Data frame: one row per |
deltas |
Numeric: the tilt factors at which the curve was evaluated. |
method |
Character: estimation method. |
n |
Integer: sample size. |
ci_level |
Numeric: confidence level. |
call |
Call: original call. |
Publicly documented design metadata for a small set of clustered (multilevel) studies. Mediation analyses in clustered settings are planned against a number of clusters, and the cluster counts that are actually attainable in a real study are often much smaller than a study's nominal sample. This reference table records both, so that planning can start from documented figures rather than round numbers.
multilevel_designsmultilevel_designs
A data frame with 3 rows and 13 columns:
Character. Study or data-source name.
Character. Cohort label, or NA when the source is not
organized by cohort.
Character. "probability sample" for a designed
sample with a documented frame, "teaching data" for a convenience
data set shipped for instruction.
Integer. Number of clusters (schools) in the
original drawn sample, before any substitution, or NA when the
source does not document one.
Integer. Number of sampled clusters found
eligible for data collection, or NA when the source does not
document one.
Integer. Total number of clusters approached,
including any substitutes added after the original sample, or NA
when the source does not document one.
Integer. Number of clusters that actually
participated. This count may include substitute clusters, so it does
not share a base with clusters_sampled.
Integer. Total number of level-1 units (children or students).
Numeric. units_total / clusters_participating. Derived here, not taken from the source; note
that the denominator is the participating cluster count, not the
sampled count.
Numeric. Cluster-level response rate as
reported by the source, or NA when none is given. The ECLS-K:1998-99
value is the weighted school response rate.
Logical. TRUE when a restricted-use data licence
is required to obtain cluster (school) identifiers, which a multilevel
analysis needs.
Character. Source-specific qualifications: which
counts are approximate, which alternative counts the source reports,
and why a field is NA.
Character. Citation, with a URL where the source is available online.
NCES names the second cohort's study "ECLS-K:2011"; its cohort value here
is "2010-11", the kindergarten class year. The two labels refer to the
same cohort.
Three points about how the numbers should be read.
Recruitment chain. clusters_sampled, clusters_eligible,
clusters_recruited, and clusters_participating are successive stages
with different bases. For ECLS-K:2011, 1,352 schools were sampled, 1,264
were eligible, 93 substitutes brought recruitment attempts to 1,446, and
roughly 970 participated. Because the participating count includes
substitutes, clusters_participating / clusters_sampled is not a
response rate; cluster_response_rate is NA on that row for this
reason, and only the rate a source reports directly is carried.
Participating versus sampled. For both ECLS-K cohorts the participating cluster count is well below the sampled count, and it is the participating count that describes the data a user would actually analyze. The NCES Handbook of Survey Methods summary page for ECLS-K:2011 implies roughly 1,310 schools; the base-year user's manual (NCES 2015-074) gives roughly 970 participating schools. The two figures come from different documents and are not interchangeable. This table carries the user's-manual figure.
Approximate counts. The ECLS-K:2011 participating-school and child counts are reported as approximate in the user's manual and are stored here as plain numbers; treat them as approximate.
ECLS-K:1998-99: NCES Handbook of Survey Methods, https://nces.ed.gov/statprog/handbook/ecls_k_dataquality.asp
ECLS-K:2011: User's Manual for the ECLS-K:2011 Kindergarten Data File (NCES 2015-074), https://nces.ed.gov/pubs2015/2015074.pdf
mediation::student: the mediation R package.
multilevel_designs[, c("study", "cohort", "clusters_participating")] # Sources obtainable without a restricted-use licence subset(multilevel_designs, !restricted_use)multilevel_designs[, c("study", "cohort", "clusters_participating")] # Sources obtainable without a restricted-use licence subset(multilevel_designs, !restricted_use)
Compute , a scale-free probabilistic effect size for mediation
analysis, along with the traditional Indirect Effect ().
Provides point estimates and bootstrap confidence intervals.
represents
:
the probability that the outcome under treatment with the mediator at its
treated level exceeds the outcome under treatment with the mediator at its
control level (a one-half correction is added for ties). Both potential
outcomes hold treatment at , so the direct effect cancels and
reflects mediation only; values indicate positive
mediation.
pmed(object, ...)pmed(object, ...)
object |
Either:
|
... |
Additional arguments passed to methods (see method documentation) |
PmedResult object containing:
estimate: P_med point estimate
ci_lower, ci_upper: Confidence interval bounds
ie_estimate: Indirect Effect point estimate
ie_ci_lower, ie_ci_upper: IE confidence interval
boot_estimates: Bootstrap distribution (if applicable)
method: Inference method used
Available inference methods via method argument:
"parametric_bootstrap" (default): Fast, assumes normality
"nonparametric_bootstrap": Robust to assumptions, slower
"plugin": Point estimate only, no CI
extract_mediation for extracting from fitted models
# Basic example with formula interface set.seed(123) n <- 200 data <- data.frame(X = rnorm(n), C = rnorm(n)) data$M <- 0.5 * data$X + 0.3 * data$C + rnorm(n) data$Y <- 0.4 * data$M + 0.2 * data$X + 0.2 * data$C + rnorm(n) # Compute with parametric bootstrap result <- pmed( Y ~ X + M + C, formula_m = M ~ X + C, data = data, treatment = "X", mediator = "M", method = "parametric_bootstrap", n_boot = 500 ) print(result) # View bootstrap distribution summary(result) plot(result)# Basic example with formula interface set.seed(123) n <- 200 data <- data.frame(X = rnorm(n), C = rnorm(n)) data$M <- 0.5 * data$X + 0.3 * data$C + rnorm(n) data$Y <- 0.4 * data$M + 0.2 * data$X + 0.2 * data$C + rnorm(n) # Compute with parametric bootstrap result <- pmed( Y ~ X + M + C, formula_m = M ~ X + C, data = data, treatment = "X", mediator = "M", method = "parametric_bootstrap", n_boot = 500 ) print(result) # View bootstrap distribution summary(result) plot(result)
A shared, minimal sensitivity helper for the P_med program (Papers 1–3).
All three estimands have the ratio form P_med = indirect / total. The
single shared point of failure is assumption A4 (no exposure-induced
mediator–outcome confounding): it biases the cross-world corner mean that
enters the indirect (numerator) component, leaving the total effect
(denominator), identified from observable corner means, robust. This function
propagates a stated numerator bias b to P_med, reporting
P_med(b) = (indirect + b) / total over a grid, the tipping bias that
zeroes the indirect effect, and (optionally) the bias that pushes P_med to
a user threshold. Gives each paper one honest sensitivity row.
pmed_sensitivity( object = NULL, indirect = NULL, total = NULL, bias_grid = NULL, threshold = NA_real_, evalue = FALSE, se = NA_real_, scale = NA_character_, ... )pmed_sensitivity( object = NULL, indirect = NULL, total = NULL, bias_grid = NULL, threshold = NA_real_, evalue = FALSE, se = NA_real_, scale = NA_character_, ... )
object |
A |
indirect, total
|
Numeric: supply the components directly instead of (or
to override) extraction from |
bias_grid |
Numeric: numerator biases |
threshold |
Numeric or |
evalue |
Logical: also compute an EValue E-value summary if the
package is installed (default |
se |
Numeric or |
scale |
Character: label for the |
... |
Unused. |
The bias model is an additive offset on the mean-based indirect component.
This is exact for the difference-scale (Paper 1) and incremental
(Paper 2) shares, whose numerators are differences of (cross-world) means.
For variance-share estimands (Paper 3, Sobol P_med^{sigma^2}) the numerator
is a variance component, not a mean; the additive-offset model is not
the right object there. Use this helper for mean-based shares; for Paper 3,
supply a pre-computed numerator perturbation if you have one, and prefer the
fuller Cinelli–Hazlett-style contour (planned) for a principled
variance-component bias map.
The optional E-value summary (evalue = TRUE) is a coarse RR-oriented
robustness number; it is reported only as a familiar single figure and is
guarded behind the EValue package (Suggests).
A PmedSensitivity object.
## From components (any paper): pmed_sensitivity(indirect = 0.30, total = 0.50, threshold = 0) ## From a Paper 1 fit: set.seed(1) n <- 400; C <- rnorm(n) A <- rbinom(n, 1, plogis(0.5 * C)) M <- 0.5 * A + 0.3 * C + rnorm(n) Y <- 0.4 * A + 0.6 * M + 0.2 * C + rnorm(n) fit <- ward_residual(data.frame(A, M, Y, C)) pmed_sensitivity(fit, threshold = 0)## From components (any paper): pmed_sensitivity(indirect = 0.30, total = 0.50, threshold = 0) ## From a Paper 1 fit: set.seed(1) n <- 400; C <- rnorm(n) A <- rbinom(n, 1, plogis(0.5 * C)) M <- 0.5 * A + 0.3 * C + rnorm(n) Y <- 0.4 * A + 0.6 * M + 0.2 * C + rnorm(n) fit <- ward_residual(data.frame(A, M, Y, C)) pmed_sensitivity(fit, threshold = 0)
Compute P_med from Formula
object |
Formula for outcome model (Y ~ X + M + C) |
formula_m |
Formula for mediator model (M ~ X + C) |
data |
Data frame |
treatment |
Character: treatment variable name |
mediator |
Character: mediator variable name |
family_y |
Family for outcome model (default: gaussian()) |
family_m |
Family for mediator model (default: gaussian()) |
x_ref |
Reference treatment value (default: 0) |
x_value |
Treatment value (default: 1) |
method |
Inference method: "parametric_bootstrap", "nonparametric_bootstrap", "plugin", "mbco", or the opt-in "model_bootstrap" (see Details). |
n_boot |
Number of bootstrap samples (default: 1000) |
ci_level |
Confidence level (default: 0.95) |
seed |
Random seed for reproducibility. For single-mediator plug-in and bootstrap methods, an explicit seed resets the caller's RNG and the computation then advances it. MBCO does not use the seed. |
n_sim |
Number of cross-world pairs drawn per evaluation of the P_med
core (default 10000). Used by the plug-in and bootstrap methods per
evaluation. |
... |
Additional arguments |
The Gaussian MBCO exact null is the union of a = 0 and b = 0;
both reduced-model likelihoods are fitted and their maximum is used.
Single-mediator MBCO uses ML residual variances (SSE/n) for both its
point estimate and likelihood. This can change the point estimate relative
to the earlier residual-df scale convention.
Numerical statuses and scan resolution are in result@diagnostics.
Verified roots use tolerance at most 1e-9. The bounded-domain scan has
spacing at most half the outward-search step and can miss narrower
components; it does not prove global connectedness. Detected disconnected
sets or failed profiles return missing interval bounds. Domain-truncated
bounds remain diagnostic bounds and are not verified LR endpoints.
method = "model_bootstrap" (opt-in; not the default) simulates complete
synthetic datasets from the fitted sequential models with the observed
treatment and covariate rows held fixed: the mediator from its model, then
the outcome using the simulated mediator. Both equations are refitted on the
same rows, so coefficients and residual scales are re-estimated jointly. The
estimate is the fitted-model plug-in; the percentile interval is reported
only when at least 99% of draws are valid (finite, converged, not
separated), otherwise both bounds are NA with
result@diagnostics$interval_status = "bootstrap_incomplete". It models
uncertainty under the specified model; it does not correct misspecification.
Calibration is limited to the documented simulation cells (see NEWS);
it does not establish general nominal coverage.
method = "mbco" returns a Model-Based Constrained Optimization interval
(Tofighi & Kelley, 2020): a likelihood-ratio interval for P_med and for the
indirect effect a*b, obtained by inverting the constrained-likelihood test
rather than by resampling. It is deterministic (no n_boot, no seed) and
supports a Gaussian outcome and mediator, with covariates, and any treatment
contrast x_ref != x_value. For binary or other non-Gaussian models, use the
bootstrap methods.
For single-mediator method = "mbco", converged requires two verified
P_med LR endpoints and no detected re-entry on the scan. It does not
establish coverage calibration. The indirect-effect interval is reported
separately and may be NA on
a degenerate design (e.g. a non-finite delta-method scale for a*b) even when
the P_med interval converges; check ie_ci_lower / ie_ci_upper directly.
Composite null. No mediation is the union of a = 0 and b = 0. At the
double null a = b = 0, the profile LR at P_med = 0.5 is the minimum of the
two branch statistics, so the chi-square(1) interval is conservative:
calibration coverage was 0.997 to 0.998 at the 95% level, and about 0.994
near the origin. @diagnostics$pmed$branch_lr holds the LR statistics for
a = 0 and b = 0. @diagnostics$pmed$near_double_null is TRUE when both
are below the cutoff (neither path individually significant), and print()
then notes that the interval is conservative. Minimum observed coverage
on the tested lattice was 0.9415, not a uniform coverage guarantee.
It stays conservative until one path is about
4 standard errors from zero, so a single null with a weak other path is
conservative too, and the flag can be FALSE there.
PmedResult object
# Toy example: Simple mediation model # Generate data where X affects Y through M set.seed(123) n <- 100 data <- data.frame( X = rnorm(n), C = rnorm(n) ) data$M <- 0.5 * data$X + 0.3 * data$C + rnorm(n, sd = 0.5) data$Y <- 0.4 * data$M + 0.2 * data$X + 0.2 * data$C + rnorm(n, sd = 0.5) # Compute P_med using plugin estimator (fast, no CI) result_plugin <- pmed( Y ~ X + M + C, formula_m = M ~ X + C, data = data, treatment = "X", mediator = "M", method = "plugin" ) print(result_plugin) # With parametric bootstrap for confidence intervals result_boot <- pmed( Y ~ X + M + C, formula_m = M ~ X + C, data = data, treatment = "X", mediator = "M", method = "parametric_bootstrap", n_boot = 200, # Use more (e.g., 1000+) in practice seed = 456 ) print(result_boot)# Toy example: Simple mediation model # Generate data where X affects Y through M set.seed(123) n <- 100 data <- data.frame( X = rnorm(n), C = rnorm(n) ) data$M <- 0.5 * data$X + 0.3 * data$C + rnorm(n, sd = 0.5) data$Y <- 0.4 * data$M + 0.2 * data$X + 0.2 * data$C + rnorm(n, sd = 0.5) # Compute P_med using plugin estimator (fast, no CI) result_plugin <- pmed( Y ~ X + M + C, formula_m = M ~ X + C, data = data, treatment = "X", mediator = "M", method = "plugin" ) print(result_plugin) # With parametric bootstrap for confidence intervals result_boot <- pmed( Y ~ X + M + C, formula_m = M ~ X + C, data = data, treatment = "X", mediator = "M", method = "parametric_bootstrap", n_boot = 200, # Use more (e.g., 1000+) in practice seed = 456 ) print(result_boot)
Computes the joint P_med for k parallel mediators (all moved together from
their control to their treated levels). Gaussian outcome and mediators only.
See vignette("parallel-mediation", package = "probmed") for the estimand.
S7 class for P_med computation results.
PmedResult( estimate = integer(0), ci_lower = integer(0), ci_upper = integer(0), ci_level = integer(0), method = character(0), n_boot = NA_integer_, boot_estimates = numeric(0), ie_estimate = NA_real_, ie_ci_lower = NA_real_, ie_ci_upper = NA_real_, ie_boot_estimates = numeric(0), x_ref = integer(0), x_value = integer(0), source_extract = NULL, converged = logical(0), diagnostics = list(), call = NULL )PmedResult( estimate = integer(0), ci_lower = integer(0), ci_upper = integer(0), ci_level = integer(0), method = character(0), n_boot = NA_integer_, boot_estimates = numeric(0), ie_estimate = NA_real_, ie_ci_lower = NA_real_, ie_ci_upper = NA_real_, ie_boot_estimates = numeric(0), x_ref = integer(0), x_value = integer(0), source_extract = NULL, converged = logical(0), diagnostics = list(), call = NULL )
estimate |
Numeric: P_med point estimate |
ci_lower |
Numeric: lower bound of confidence interval |
ci_upper |
Numeric: upper bound of confidence interval |
ci_level |
Numeric: confidence level (e.g., 0.95) |
method |
Character: inference method used |
n_boot |
Integer: number of bootstrap samples (NA if not bootstrap) |
boot_estimates |
Numeric vector: bootstrap distribution (empty if not bootstrap) |
ie_estimate |
Numeric: Indirect Effect (NIE) point estimate |
ie_ci_lower |
Numeric: lower bound of NIE confidence interval |
ie_ci_upper |
Numeric: upper bound of NIE confidence interval |
ie_boot_estimates |
Numeric vector: bootstrap distribution of NIE |
x_ref |
Numeric: reference treatment value |
x_value |
Numeric: treatment value for contrast |
source_extract |
medfit::MediationData object: source of the estimates |
converged |
Logical: did computation converge |
diagnostics |
List: method-specific numerical interval diagnostics |
call |
Call object: original function call |
S7 class holding a one-row sensitivity analysis for a proportion-mediated
estimand under violation of assumption A4 (no exposure-induced
mediator–outcome confounding). Every estimand in the P_med family shares
the ratio form P_med = indirect / total. A4 failure biases the
cross-world corner mean that enters the indirect (numerator)
component, while the total (denominator) effect is identified from
observable corner means and is robust to A4. The object therefore reports how
P_med moves as the numerator bias b ranges over a stated set, plus the
tipping bias that drives the indirect component (or P_med) through a
threshold.
PmedSensitivity( p_med = integer(0), indirect = integer(0), total = integer(0), bias_grid = integer(0), p_med_bias = integer(0), tipping_indirect = integer(0), tipping_threshold = NA_real_, threshold = NA_real_, evalue = numeric(0), scale = NA_character_, call = NULL )PmedSensitivity( p_med = integer(0), indirect = integer(0), total = integer(0), bias_grid = integer(0), p_med_bias = integer(0), tipping_indirect = integer(0), tipping_threshold = NA_real_, threshold = NA_real_, evalue = numeric(0), scale = NA_character_, call = NULL )
p_med |
Numeric: the point estimate |
indirect |
Numeric: the indirect (mediated) component, the numerator. |
total |
Numeric: the total effect, the denominator. |
bias_grid |
Numeric: the grid of numerator biases |
p_med_bias |
Numeric: |
tipping_indirect |
Numeric: bias |
tipping_threshold |
Numeric: bias |
threshold |
Numeric: the |
evalue |
Numeric: optional E-value summary (length-2: estimate and CI
bound) from the EValue package; |
scale |
Character: label of the |
call |
Call: the originating call. |
Standalone cross-fit DR estimator of P_med^W. Consistent when EITHER the
outcome model OR the propensity model is misspecified, provided the mediator
model is correct (multiply robust; verified in simulation).
pmedW_dr(d, covars = "C", K = 5L, misspec = "none", G = 300L, seed = 1L)pmedW_dr(d, covars = "C", K = 5L, misspec = "none", G = 300L, seed = 1L)
d |
A |
covars |
Character vector of covariate column names. Default |
K |
Integer: cross-fitting folds (default 5). |
misspec |
Character: |
G |
Integer: number of |
seed |
Integer: RNG seed for fold assignment (default 1). |
A named list with NIE_W, NDE_W, and pmedW.
Computes P_med^W for multivariate outcomes (d >= 1) from corner pseudo-samples
via debiased Sinkhorn divergence. Intended for use with pre-computed corner-law draws
(e.g. from g-computation or DR). SE not yet implemented.
pmedW_md(NU11, NU10, NU00, eps = 0.1)pmedW_md(NU11, NU10, NU00, eps = 0.1)
NU11, NU10, NU00
|
Numeric matrices (rows = draws, columns = outcome dimensions). |
eps |
Numeric: Sinkhorn regularization parameter (default 0.1). |
A named list with NIE_W, NDE_W, and pmedW.
m
Applies the coarse-graining operator T_m: within each cluster, units are
randomly partitioned into blocks of size m and each block is replaced by
its mean. The cluster structure is preserved, so every scale keeps the same
clusters and two scales remain comparable cluster-by-cluster.
rg_coarse_grain( data, m, treatment = "A", mediator = "M", outcome = "Y", cluster = "cluster" )rg_coarse_grain( data, m, treatment = "A", mediator = "M", outcome = "Y", cluster = "cluster" )
data |
Data frame with the treatment, mediator, outcome, and cluster columns. |
m |
Positive integer block size. |
treatment, mediator, outcome
|
Character: column names of |
cluster |
Character: column name of the cluster identifier. |
Blocks are complete by construction: within a cluster of size n_s, the
final incomplete block is discarded, so the number of coarse-grained units
is sum(floor(n_s / m)), not n / m. Clusters smaller than m drop out
entirely at that scale. The composition property
T_{m2}(T_{m1}(.)) = T_{m1 * m2}(.) holds exactly on the units a partition
retains.
This function calls sample.int() and is therefore random. It never sets a
seed; seed the caller (or wrap in withr::with_seed()) for reproducibility.
List with Z (numeric matrix of coarse-grained units, columns
A, M, Y) and cluster (their cluster labels).
set.seed(1) d <- data.frame( A = rnorm(200), M = rnorm(200), Y = rnorm(200), sch = rep(1:20, each = 10) ) cg <- rg_coarse_grain(d, m = 5, cluster = "sch") nrow(cg$Z) # 40 = 20 clusters * floor(10 / 5)set.seed(1) d <- data.frame( A = rnorm(200), M = rnorm(200), Y = rnorm(200), sch = rep(1:20, each = 10) ) cg <- rg_coarse_grain(d, m = 5, cluster = "sch") nrow(cg$Z) # 40 = 20 clusters * floor(10 / 5)
Contrasts g at two aggregation scales estimated from the same clusters
by differencing their influence functions cluster-by-cluster. Because the
shared sampling noise cancels, the paired standard error is materially
smaller than the naive one that treats the two scales as independent.
rg_flow_contrast(fit1, fit2, ci_level = 0.95)rg_flow_contrast(fit1, fit2, ci_level = 0.95)
fit1, fit2
|
RgFlowResult objects estimated on the same clusters. The
contrast is |
ci_level |
Numeric: confidence level. |
The gain is largest for two nearby scales, where the influence functions are most similar cluster-by-cluster; the endpoint contrast (unit scale versus fully aggregated) is the worst case for pairing.
List with delta, se_paired, se_naive, gain
(se_naive / se_paired), z, p_value, ci, and n_clusters.
set.seed(1) d <- data.frame( A = rnorm(400), M = rnorm(400), Y = rnorm(400), sch = rep(1:40, each = 10) ) f1 <- g_hat_eif(rg_coarse_grain(d, 1, cluster = "sch")$Z, rg_coarse_grain(d, 1, cluster = "sch")$cluster, m = 1) cg2 <- rg_coarse_grain(d, 2, cluster = "sch") f2 <- g_hat_eif(cg2$Z, cg2$cluster, m = 2) rg_flow_contrast(f1, f2)$se_pairedset.seed(1) d <- data.frame( A = rnorm(400), M = rnorm(400), Y = rnorm(400), sch = rep(1:40, each = 10) ) f1 <- g_hat_eif(rg_coarse_grain(d, 1, cluster = "sch")$Z, rg_coarse_grain(d, 1, cluster = "sch")$cluster, m = 1) cg2 <- rg_coarse_grain(d, 2, cluster = "sch") f2 <- g_hat_eif(cg2$Z, cg2$cluster, m = 2) rg_flow_contrast(f1, f2)$se_paired
S7 class for the scale-indexed proportion mediated
g_l = NIE_l / TE_l estimated on data coarse-grained to aggregation scale
l (block size m). Under the coarse-graining operator T_m – which
replaces each within-cluster block of m units by its mean – g_l flows
with m, and the flow is governed by the mediation beta-function.
This is the classical proportion mediated (a ratio of effects), not
the package's P_med (a probability in [0, 1]); the two are different
quantities and are deliberately named apart. See vignette("probmed") for
P_med.
RgFlowResult( g = integer(0), g_ci = integer(0), g_se = integer(0), a_path = integer(0), b_path = integer(0), TE = integer(0), grad_norm = integer(0), phi_cluster = integer(0), m = integer(0), n_units = integer(0), n_clusters = integer(0), ci_level = integer(0), call = NULL )RgFlowResult( g = integer(0), g_ci = integer(0), g_se = integer(0), a_path = integer(0), b_path = integer(0), TE = integer(0), grad_norm = integer(0), phi_cluster = integer(0), m = integer(0), n_units = integer(0), n_clusters = integer(0), ci_level = integer(0), call = NULL )
g |
Numeric: estimated proportion mediated at this scale. |
g_ci |
Numeric length-2: cluster-robust EIF Wald confidence interval. |
g_se |
Numeric: cluster-robust standard error. |
a_path, b_path, TE
|
Numeric: the |
grad_norm |
Numeric: Euclidean norm of the gradient of |
phi_cluster |
Numeric: per-cluster influence contributions, in cluster
order. Retained so two scales estimated on the same clusters can be
contrasted by |
m |
Numeric: aggregation scale (block size); |
n_units |
Integer: coarse-grained units retained at this scale. |
n_clusters |
Integer: clusters contributing at this scale. |
ci_level |
Numeric: confidence level. |
call |
Call: original call. |
g_l is a smooth functional of five second moments of (A, M, Y), so the
efficient influence function is available in closed form and the plug-in
estimator is already efficient in the linear-Gaussian case (no cross-fitting
required). The sampling unit is the cluster, not the unit: influence
contributions are aggregated to the cluster before their variance is taken.
Debiased Sinkhorn divergence S_eps(X, Y) (Genevay et al.), which approximates
W_2^2(X, Y) as eps -> 0. Use for multivariate (d > 1) Wasserstein P_med^W
via pmedW_md().
sinkhorn_div(X, Y, eps = 0.1, iters = 2000L)sinkhorn_div(X, Y, eps = 0.1, iters = 2000L)
X, Y
|
Numeric matrices (rows = samples, columns = dimensions). |
eps |
Numeric: regularization parameter (default 0.1). |
iters |
Integer: maximum Sinkhorn iterations (default 400). |
Non-negative scalar.
Closed-form map from the four corner means theta(a, a') to the Sobol /
functional-ANOVA variance components and the variance-scale proportion
mediated. Used internally by sobol_pmed() and exported as a helper for
computing the analytic (population) truth from known corner means.
sobol_from_theta(theta, pd = 0.5, pm = 0.5)sobol_from_theta(theta, pd = 0.5, pm = 0.5)
theta |
Numeric length-4: corner means named |
pd, pm
|
Numeric: tilt probabilities for the direct and mediated pathways
(Bernoulli design points); defaults |
A named numeric vector with elements Vd, Vm, Vdm, VT,
Pmed_sobol (= Vm / VT), S1_med (= Vm / VT),
ST_med (= (Vm + Vdm) / VT), and Rint (the interaction remainder).
th <- c("11" = 1.32, "10" = 0.5, "01" = 0.42, "00" = 0) sobol_from_theta(th)th <- c("11" = 1.32, "10" = 0.5, "01" = 0.42, "00" = 0) sobol_from_theta(th)
Estimate the Sobol proportion mediated P_med^{sigma2} = V_med / V_T – the
share of the total interventional outcome variance carried by the mediated
pathway, i.e. the first-order Sobol (functional-ANOVA) sensitivity index of the
outcome with respect to the mediator. The estimator is cross-fitted and built
on the same triply-robust corner pseudo-outcomes as ward_residual() and
incr_pmed(); the standard error uses the delta method through the quadratic
variance components and the ratio identity. See SobolPmedResult for the
variance decomposition and the boundary remedy.
sobol_pmed( object, pd = 0.5, pm = 0.5, covars = "C", K = 5L, seed = 1L, ci_level = 0.95, warn_boundary = TRUE, boundary_test = c("split", "plugin", "none"), procedure = c("B", "A"), reps = 1L, se_method = c("analytic", "bootstrap"), B = 200L, ... )sobol_pmed( object, pd = 0.5, pm = 0.5, covars = "C", K = 5L, seed = 1L, ci_level = 0.95, warn_boundary = TRUE, boundary_test = c("split", "plugin", "none"), procedure = c("B", "A"), reps = 1L, se_method = c("analytic", "bootstrap"), B = 200L, ... )
object |
A |
pd, pm
|
Numeric: tilt probabilities for the direct and mediated pathways
(Bernoulli design points); defaults |
covars |
Character vector of covariate column names. Default |
K |
Integer: number of cross-fitting folds (default |
seed |
Integer: RNG seed for fold assignment and the boundary split
(default |
ci_level |
Numeric: confidence level (default |
warn_boundary |
Logical: warn when the boundary test does not reject and
the reported interval falls back to the Procedure-A upper bound
(default |
boundary_test |
Character: boundary test for |
procedure |
Character: interval procedure for |
reps |
Integer: repeated cross-fitting draws averaged for the point estimate
(default |
se_method |
Character: |
B |
Integer: bootstrap resamples when |
... |
Unused. |
At the null V_med = 0 (equivalently Delta_m = 0) the squared-variance
influence function degenerates and the symmetric Wald interval is non-regular.
A Williamson et al. (2021) sample-split test on the regular contrast Delta_m
decides the boundary. The default reported interval is Procedure B
(procedure = "B"): the image of the regular Delta_m Wald CI under
delta -> (c_m/V_T) delta^2 (continuous mapping, no gate). Its validity inherits
that of the input Delta_m CI, which the analytic se delivers: post-fix
(2026-08-22) the Delta_m se is exact from strong effects to the null and the
Procedure-B interval covers 0.952-0.996 across six cells (see the boundary
coverage section of SobolPmedResult). Procedure A
(procedure = "A") is the legacy gated rule (Wald on rejection, one-sided
[0, Pmed_upper] otherwise), pointwise- but not uniformly-valid across the
near-null transition (Leeb-Potscher pre-test under-coverage). The boundary
machinery recurses through a plain internal fitter (not this generic), so the
point estimate, fold draw, and split are deterministic in seed. See the
boundary coverage section in SobolPmedResult.
A SobolPmedResult object.
Williamson, B. D., Gilbert, P. B., Carone, M., & Simon, N. (2021). Nonparametric variable importance assessment using machine learning techniques. Biometrics, 77(1), 9–22.
set.seed(1) n <- 1500; C <- rnorm(n) A <- rbinom(n, 1, plogis(-0.2 + 0.8 * C)) M <- 0.6 * A + 0.4 * C + rnorm(n) # strong mediation (off-boundary): Wald interval reported Y <- 0.5 * A + 0.7 * M + 0.3 * C + rnorm(n) sobol_pmed(data.frame(A, M, Y, C))set.seed(1) n <- 1500; C <- rnorm(n) A <- rbinom(n, 1, plogis(-0.2 + 0.8 * C)) M <- 0.6 * A + 0.4 * C + rnorm(n) # strong mediation (off-boundary): Wald interval reported Y <- 0.5 * A + 0.7 * M + 0.3 * C + rnorm(n) sobol_pmed(data.frame(A, M, Y, C))
S7 class for the Sobol (functional-ANOVA) proportion mediated
P_med^{sigma2} = V_med / V_T (Paper 3 of the P_med-for-modern-estimands
program). On the variance scale the proportion mediated is the first-order
Sobol sensitivity index of the outcome with respect to the mediated pathway:
the share of the total interventional outcome variance V_T explained by the
mediated component V_med. Unlike the additive split TE = NDE + NIE, the
variance decomposition V_T = V_dir + V_med + V_int always holds (the
interaction variance V_int is a separate non-negative term), so
P_med^{sigma2} is well defined even under treatment-by-mediator interaction.
SobolPmedResult( p_med = integer(0), se = integer(0), ci = integer(0), ci_wald = integer(0), ci_A = NA_real_, ci_B1 = NA_real_, ci_B2 = NA_real_, procedure = "B", boundary = logical(0), Pmed_upper = integer(0), S1_med = integer(0), ST_med = integer(0), Vd = integer(0), Vm = integer(0), Vdm = integer(0), VT = integer(0), Dm = integer(0), se_Dm = integer(0), vmed_split_p = NA_real_, vmed_split_reject = NA_integer_, theta = integer(0), method = character(0), se_method = "analytic", reps = 1L, n = integer(0), ci_level = integer(0), call = NULL )SobolPmedResult( p_med = integer(0), se = integer(0), ci = integer(0), ci_wald = integer(0), ci_A = NA_real_, ci_B1 = NA_real_, ci_B2 = NA_real_, procedure = "B", boundary = logical(0), Pmed_upper = integer(0), S1_med = integer(0), ST_med = integer(0), Vd = integer(0), Vm = integer(0), Vdm = integer(0), VT = integer(0), Dm = integer(0), se_Dm = integer(0), vmed_split_p = NA_real_, vmed_split_reject = NA_integer_, theta = integer(0), method = character(0), se_method = "analytic", reps = 1L, n = integer(0), ci_level = integer(0), call = NULL )
p_med |
Numeric: Sobol proportion mediated |
se |
Numeric: standard error of |
ci |
Numeric length-2: reported confidence interval, selected by
|
ci_wald |
Numeric length-2: ungated symmetric Wald interval (always the Wald interval, even at the boundary). |
ci_A |
Numeric length-2: Procedure-A (gated, pre-test) interval – the Wald
interval off-boundary, or |
ci_B1 |
Numeric length-2: Procedure-B interval from the same-sample
|
ci_B2 |
Numeric length-2: experimental Procedure-B variant from a
single decorrelated split ( |
procedure |
Character: |
boundary |
Logical: |
Pmed_upper |
Numeric: one-sided Procedure-A upper bound for |
S1_med |
Numeric: first-order Sobol index |
ST_med |
Numeric: total Sobol index |
Vd, Vm, Vdm, VT
|
Numeric: direct, mediated, interaction, and total variance components. |
Dm |
Numeric: mediated antecedent contrast |
se_Dm |
Numeric: standard error of |
vmed_split_p |
Numeric: p-value of the sample-split boundary test of
|
vmed_split_reject |
Integer: |
theta |
Numeric length-4: corner means |
method |
Character: estimation method. |
se_method |
Character: |
reps |
Integer: number of repeated cross-fitting fold draws averaged for the
point estimate (default |
n |
Integer: sample size. |
ci_level |
Numeric: confidence level. |
call |
Call: original call. |
Writing the four corner means as theta(a, a') = E[Y(a, M(a'))], the contrasts
are the direct antecedent Delta_d, the mediated antecedent Delta_m, and the
interaction remainder R = theta_11 - theta_10 - theta_01 + theta_00. With
Bernoulli tilt variances c_d = p_d(1 - p_d) and c_m = p_m(1 - p_m),
V_dir = c_d Delta_d^2, V_med = c_m Delta_m^2, V_int = c_d c_m R^2, and
V_T = V_dir + V_med + V_int. The estimator is the cross-fitted one-step
estimator built on the shared triply-robust corner pseudo-outcomes (the same
engine as ward_residual() and incr_pmed()); the delta method propagates the
corner influence functions through the quadratic V-components and the
ratio identity to a standard error for P_med^{sigma2}.
Boundary remedy. Because V_med = c_m Delta_m^2 is a squared quantity, its
influence function phi_{V_med} = 2 c_m Delta_m phi_{Delta_m} degenerates as
Delta_m -> 0, deflating the Wald standard error and collapsing coverage at
the null V_med = 0. The null is equivalent to Delta_m = 0, which has a
regular root-n influence function. Following Williamson et al. (2021,
Biometrics), the boundary is tested with a sample-split statistic that
estimates Delta_m on one half and its standard error on the independent other
half (so z ~ N(0, 1) even at the boundary).
Two interval procedures. The default reported ci is Procedure B: the
image of the regular two-sided Delta_m Wald CI under the map
delta -> (c_m/V_T) delta^2 (ci_B1/ci_B2). Because it is the continuous
image of a uniformly-valid CI for the regular parameter Delta_m, it is
uniformly valid in Delta_m – there is no pre-test, hence no boundary
pathology (subject to the near-null se caveat below). The reported ci is
ci_B1, the same-sample image (centered at the full Delta_m, so it always
brackets p_med); ci_B2 is an experimental single-split variant that did
not improve coverage (see its note). Procedure A (procedure = "A",
ci_A) is the legacy gated rule: Wald on rejection, one-sided
[0, Pmed_upper] otherwise.
Coverage near the boundary (A-15). Two distinct issues.
(1) Procedure A is not uniformly valid. As a pre-test (gating) rule A is only
pointwise valid: near the boundary it routes downward-selected Delta_m_hat to
the contracting one-sided bound, a Leeb-Potscher post-selection effect. Prefer B,
which does not gate. (Measured 2026-08-22: A covers 0.90-0.92 at the
transition, B 0.95-0.97.)
(2) The near-null Delta_m standard error – exact after the corner-weight
fix (2026-08-22). Against exact truth in six cells (n = 800 and 3000; strong,
intermediate and near-null Delta_m; 250 datasets each;
inst/sim/sobol_shipped_check.R, inst/sim/results/sobol_postfix_summary.csv,
docs/specs/FINDINGS-2026-08-22-sobol-pmedw-postfix.md) the analytic
influence-function se of Delta_m matches its sampling SD (ratio 0.99, CV
0.08) from strong effects to the null, with Delta_m Wald coverage 0.95-0.96;
the estimator's empirical SD is within 3% of an oracle's with known nuisances;
and p_med's Procedure-B1 interval covers 0.952-0.996 (conservative under a
strong interaction; 0.952 at the boundary with 98% of draws gated). The only
near-null pathology is the symmetric Wald interval for p_med itself (0.74 at
the boundary – the same for the oracle: the square map, not the se), which
Procedure B avoids by construction.
History. Before the fix this entry reported three "pinned facts": that a
single fold partition was ~80% of Var(Delta_m_hat) near the null, that the
analytic se was ~0.8x its sampling SD so Procedure B covered ~0.85, and that a
refit bootstrap (~1.25x, coverage ~0.97) was the remedy. All three measured the
corner-EIF weight bug fixed in PR #34 (.corner_fit() weighted every row of a
fold by its first row's propensity; NEWS 0.4.0). Post-fix reps > 1 has
nothing left to remove, and the bootstrap se is 0.84-1.05x the sampling SD
(slightly narrower at the transition); both remain options, not remedies.
Practical guidance. The analytic Procedure-B default is valid across the
measured transition; nothing needs switching near the boundary. reps > 1
makes the point estimate reproducible across fold draws (small efficiency
effect). se_method = "bootstrap" is an option costing B (x reps)
refits; post-fix it is not better than the analytic se and slightly
anti-conservative at the transition. V_T is well calibrated throughout. Use
procedure = "A" only to reproduce the legacy gated behavior.
Williamson, B. D., Gilbert, P. B., Carone, M., & Simon, N. (2021). Nonparametric variable importance assessment using machine learning techniques. Biometrics, 77(1), 9–22.
Estimate the interventional proportion mediated together with the gauge
residual W = R / OE that flags non-decomposability (treatment-by-mediator
interaction). Uses a cross-fitted one-step estimator built on the
triply-robust efficient influence functions of the four corner means
theta(a, a') = E[Y(a, M(a'))].
ward_residual( object, covars = "C", K = 5L, ci_level = 0.95, seed = 1L, fieller = TRUE, reps = 1L, se_method = c("analytic", "bootstrap"), B = 200L, a0 = 0, a1 = 1, weak_id_ratio_threshold = 3, oe_snr_threshold = 2, ... )ward_residual( object, covars = "C", K = 5L, ci_level = 0.95, seed = 1L, fieller = TRUE, reps = 1L, se_method = c("analytic", "bootstrap"), B = 200L, a0 = 0, a1 = 1, weak_id_ratio_threshold = 3, oe_snr_threshold = 2, ... )
object |
A |
covars |
Character vector of covariate column names. Default |
K |
Integer: number of cross-fitting folds (default 5). |
ci_level |
Numeric: confidence level (default 0.95). |
seed |
Integer: RNG seed for fold assignment. |
fieller |
Logical: also compute the Fieller confidence set for the ratio
|
reps |
Integer: number of repeated cross-fitting fold draws averaged for
the point estimate (default |
se_method |
Character: History. Every coverage figure this package reported before
2026-08-22 – the 2,000-rep coverage grid
( |
B |
Integer: number of bootstrap resamples when |
a0, a1
|
Reference and comparison exposure levels (defaults |
weak_id_ratio_threshold |
Numeric: the |
oe_snr_threshold |
Numeric: the |
... |
Unused. |
A GaugePmedResult object.
set.seed(1) n <- 800; C <- rnorm(n) A <- rbinom(n, 1, plogis(-0.2 + 0.8 * C)) M <- 0.6 * A + 0.4 * C + rnorm(n) Y <- 0.5 * A + 0.7 * M + 0.8 * A * M + 0.3 * C + rnorm(n) ward_residual(data.frame(A, M, Y, C))set.seed(1) n <- 800; C <- rnorm(n) A <- rbinom(n, 1, plogis(-0.2 + 0.8 * C)) M <- 0.6 * A + 0.4 * C + rnorm(n) Y <- 0.5 * A + 0.7 * M + 0.8 * A * M + 0.3 * C + rnorm(n) ward_residual(data.frame(A, M, Y, C))
Estimate P_med^W = NIE^W / (NIE^W + NDE^W), the transport-scale proportion
mediated, where NIE^W = W_2(nu_11, nu_10) and NDE^W = W_2(nu_10, nu_00) are
the 2-Wasserstein distances between the corner laws nu(a, a').
Three estimation methods are available:
"eif" (default): deterministic parametric g-computation with an
observed-row influence-function SE. Includes both fitted-model and empirical
covariate-distribution uncertainty. Requires correct working models, iid
rows, full-rank fits, and nonzero components; d = 1 only. The legacy name
does not assert nonparametric efficiency.
"bootstrap": row bootstrap refitting the same deterministic estimator.
Used automatically near either zero-component boundary. Boundary coverage
is not guaranteed by this fallback.
"dr": cross-fit DR corner-CDF estimator (multiply robust: consistent when
either the outcome or propensity model is wrong, mediator model correct) with
bootstrap SE. Slowest; d = 1 only.
wasserstein_pmed( object, covars = "C", n_mc = 20000L, binY = FALSE, method = c("eif", "bootstrap", "dr"), n_boot = 200L, K = 5L, G = 300L, ci_level = 0.95, seed = 1L, ... )wasserstein_pmed( object, covars = "C", n_mc = 20000L, binY = FALSE, method = c("eif", "bootstrap", "dr"), n_boot = 200L, K = 5L, G = 300L, ci_level = 0.95, seed = 1L, ... )
object |
A |
covars |
Character vector of covariate column names. Default |
n_mc |
Legacy compatibility argument (default 20000), now ignored. G-computation uses deterministic quadrature instead of simulated corner laws. |
binY |
Logical: if |
method |
Character: |
n_boot |
Integer: bootstrap replications (default 200). Used by
|
K |
Integer: cross-fitting folds for |
G |
Integer: y-grid points for |
ci_level |
Numeric: confidence level (default 0.95). |
seed |
Integer: RNG seed (default 1). |
... |
Unused. |
The g-computation methods assume a homoscedastic Gaussian mediator model
M ~ A + covars and an outcome model Y ~ A * M + covars (homoscedastic
Gaussian, or Bernoulli-logit when binY = TRUE). Covariates are standardized
over the observed rows. Residual variances use mean squared residuals, not
degrees-of-freedom corrections. Gaussian outcomes use 128-point quadrature
over mixture quantiles; binary outcomes use 64-point normal quadrature.
These are numerical approximations, not certified error bounds. Wald
inference is model-based and requires regular, nonzero distance components;
it is not a method for clustered observations or unsupported outcome families.
A WassersteinPmedResult object.
set.seed(1) n <- 800; C <- rnorm(n) A <- rbinom(n, 1, plogis(0.3 * C)) M <- 0.7 * A + 0.5 * C + rnorm(n) Y <- 0.6 * A + 0.8 * M + 0.3 * C + rnorm(n) wasserstein_pmed(data.frame(A, M, Y, C), covars = "C")set.seed(1) n <- 800; C <- rnorm(n) A <- rbinom(n, 1, plogis(0.3 * C)) M <- 0.7 * A + 0.5 * C + rnorm(n) Y <- 0.6 * A + 0.8 * M + 0.3 * C + rnorm(n) wasserstein_pmed(data.frame(A, M, Y, C), covars = "C")
S7 class for the Wasserstein proportion mediated
P_med^W = NIE^W / (NIE^W + NDE^W), where NIE^W = W_2(nu_11, nu_10) and
NDE^W = W_2(nu_10, nu_00) are the 2-Wasserstein distances between the
corner laws nu(a, a') = Law(Y(a, M(a'))). Unlike the additive split
TE = NDE + NIE, the Wasserstein decomposition handles treatment-by-mediator
interaction via the synergy term G^W = NIE^W + NDE^W - TE^W >= 0
(triangle gap).
WassersteinPmedResult( pmedW = integer(0), pmedW_ci = integer(0), pmedW_se = integer(0), NIE_W = integer(0), NDE_W = integer(0), TE_W = integer(0), synergy = integer(0), boundary = FALSE, method = character(0), n = integer(0), ci_level = integer(0), call = NULL )WassersteinPmedResult( pmedW = integer(0), pmedW_ci = integer(0), pmedW_se = integer(0), NIE_W = integer(0), NDE_W = integer(0), TE_W = integer(0), synergy = integer(0), boundary = FALSE, method = character(0), n = integer(0), ci_level = integer(0), call = NULL )
pmedW |
Numeric: transport-scale proportion mediated. |
pmedW_ci |
Numeric length-2: confidence interval (Wald or bootstrap). |
pmedW_se |
Numeric: standard error. |
NIE_W, NDE_W, TE_W
|
Numeric: Wasserstein NIE, NDE, and total effect. |
synergy |
Numeric: triangle gap |
boundary |
Logical: |
method |
Character: estimation method ( |
n |
Integer: sample size. |
ci_level |
Numeric: confidence level. |
call |
Call: original call. |