| 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 <dtofighi@gmail.com> |
| License: | MIT + file LICENSE |
| Version: | 0.3.0 |
| Built: | 2026-10-14 19:05:16 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
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), 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), 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 |
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 * gamma_med - med * gamma_dir) / tot^2
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.
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))
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. |
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) ## Not run: # View bootstrap distribution summary(result) plot(result) ## End(Not run)# 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) ## Not run: # View bootstrap distribution summary(result) plot(result) ## End(Not run)
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" |
n_boot |
Number of bootstrap samples (default: 1000) |
ci_level |
Confidence level (default: 0.95) |
seed |
Random seed for reproducibility |
... |
Additional arguments |
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 method = "mbco", the converged flag reflects the P_med interval
only. 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.
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), 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), 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 |
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. |
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 – so near the boundary use reps > 1 and
se_method = "bootstrap" for a valid interval (the analytic Delta_m se is
anti-conservative there; see the boundary coverage section). 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 – mechanism and remedies (A-15). Two distinct issues; the second is now decomposed (reproducible fixed-seed variance decomposition) into three findings, each with a concrete remedy.
(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.
(2) The near-null Delta_m standard error. Three pinned facts:
Fold-split Monte-Carlo variance dominates. A single cross-fitting
fold partition contributes ~80% of Var(Delta_m_hat) in the near-null regime
(~45% at ordinary effect sizes), persistent across n = 2000 to 8000. This
is an algorithmic nuisance, not sampling information: it inflates the
estimator's variance and makes the single-split point estimate seed-dependent.
Remedy: reps > 1 (repeated cross-fitting – the corner influence matrix
is averaged over reps independent fold draws), which removes it and yields a
reproducible, lower-variance point estimate (mean DML aggregation).
Delta_m_hat is approximately normal at the null. Recomputing coverage
with the oracle Monte-Carlo SD gives ~0.95, so the Wald shape is correct for
the regular contrast Delta_m: the under-coverage is a wrong interval width,
not a wrong shape (there is no boundary non-normality for Delta_m itself).
A residual analytic-se bias remains. Even after fold noise is removed,
the analytic influence-function se for Delta_m is ~0.8x its true sampling SD
near the null (structural; persistent in n); it is well calibrated (ratio
~1) at ordinary effect sizes. So the default analytic near-null interval –
and its Procedure-B image – covers only ~0.85. Remedy: se_method = "bootstrap" (resample rows, refit the cross-fit estimator), which recovers a
valid but mildly conservative se (~1.25x the true SD, coverage ~0.97 near
null) and is calibrated off-boundary.
Practical guidance. Off the boundary the analytic default is calibrated and
fast. Near the boundary (the split test does not reject), pass reps > 1 for a
reproducible point and se_method = "bootstrap" for a valid interval (conservative
there – it trades width for guaranteed coverage). V_T is well calibrated
throughout. The full coverage grid across the Delta_m transition is validated by
a separate large simulation. Use procedure = "A" only to reproduce the legacy
gated behaviour.
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, ... )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, ... )
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: |
B |
Integer: number of bootstrap resamples when |
a0, a1
|
Reference and comparison exposure levels (defaults |
... |
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))