probmed

probmed provides a robust framework for computing $P_{med}$, a scale-free probabilistic effect size for causal mediation analysis.

What is $P_{med}$?

$P_{med}$ offers an intuitive alternative to traditional mediation effect sizes. It represents the probability that the outcome for an individual in the treatment group is higher than the outcome for an individual in the control group, specifically due to the indirect path through the mediator.

Unlike the “Proportion Mediated” or simple indirect effect coefficients, $P_{med}$: * Is Scale-Free: Interpretable regardless of the units of measurement for $X$, $M$, or $Y$. * Is Bounded: Always falls between 0 and 1 (or 0 and 0.5 for some definitions), making it comparable across studies. * Handles Non-Linearity: Can be extended to Generalized Linear Models (GLMs) where coefficients are not directly comparable.

Features

Mediationverse Ecosystem

probmed is part of the mediationverse ecosystem for mediation analysis in R:

Package Purpose Role
medfit Model fitting, extraction, bootstrap Foundation
probmed (this) Probabilistic effect size (P_med) Application
RMediation Confidence intervals (DOP, MBCO) Application
medrobust Sensitivity analysis Application
medsim Simulation infrastructure Support

See Ecosystem Coordination for version compatibility and development guidelines.

Installation

Install from the Data-Wise r-universe — no compiler needed, and dependencies (e.g. medfit) resolve automatically:

install.packages("probmed", repos = "https://data-wise.r-universe.dev")

Or install the development version from GitHub:

# install.packages("devtools")
devtools::install_github("data-wise/probmed")

Quick Example

Here is a basic example showing how to simulate data and compute $P_{med}$ using a parametric bootstrap approach.

1. Simulate Data

First, we generate a dataset with a continuous treatment ($X$), mediator ($M$), and outcome ($Y$), along with a covariate ($C$).

library(probmed)

set.seed(123)
n <- 300
data <- data.frame(X = rnorm(n), C = rnorm(n))

# Mediator model: M ~ X + C
data$M <- 0.5 * data$X + 0.3 * data$C + rnorm(n)

# Outcome model: Y ~ M + X + C
data$Y <- 0.4 * data$M + 0.2 * data$X + 0.2 * data$C + rnorm(n)

2. Compute $P_{med}$

We use the pmed() function to estimate the effect size. We specify the models for the outcome and mediator, and select the bootstrap method.

result <- pmed(
  Y ~ X + M + C,             # Outcome model formula
  formula_m = M ~ X + C,     # Mediator model formula
  data = data,
  treatment = "X",
  mediator = "M",
  method = "parametric_bootstrap",
  n_boot = 1000,
  seed = 42
)

3. View Results

The results object provides a clear summary of both $P_{med}$ and the Indirect Effect with confidence intervals.

print(result)
#>
#> P_med: Probability of Mediated Shift
#> ====================================
#>
#> Estimate: 0.563
#> 95% CI: [0.520, 0.605]
#>
#> Indirect Effect (Product of Coefficients):
#> Estimate: 0.198
#> 95% CI: [0.134, 0.268]
#>
#> Inference: parametric_bootstrap
#> Bootstrap samples: 1000
#>
#> Treatment contrast: X = 1 vs. X* = 0
#>
#> Interpretation:
#>   P(Y(1, M(1)) > Y(1, M(0))) = 0.563
#>   P that the mediator shift (M(0) -> M(1)) leaves a random individual better off, holding X = 1.

summary(result)

Interpretation:

Workflow Options

probmed offers three workflow approaches to suit your needs:

1. Formula Interface (Recommended for Most Users)

Directly specify your models using R formulas:

result <- pmed(
  Y ~ X + M + C,             # Outcome model
  formula_m = M ~ X + C,     # Mediator model
  data = data,
  treatment = "X",
  mediator = "M",
  method = "parametric_bootstrap"
)

2. Integration with lavaan (For SEM Users)

Extract from fitted Structural Equation Models with full support for FIML, robust estimators, and standardized estimates:

library(lavaan)

# Define SEM model
model <- '
  M ~ a*X + C
  Y ~ b*M + cp*X + C
'

# Fit with FIML for missing data
fit <- sem(model, data = data, missing = "fiml")

# Extract and compute P_med
extract <- extract_mediation(fit, treatment = "X", mediator = "M")
result <- pmed(extract, method = "parametric_bootstrap", n_boot = 5000)
print(result)

3. Integration with mediation Package

Seamlessly work with objects from the mediation package:

library(mediation)

# Fit models
model_m <- lm(M ~ X + C, data = data)
model_y <- lm(Y ~ M + X + C, data = data)

# Run mediate
med_out <- mediate(model_m, model_y, treat = "X", mediator = "M", boot = TRUE)

# Compute P_med from mediate object
result <- pmed(extract_mediation(med_out))
print(result)

Advanced Features

Generalized Linear Models (GLMs)

probmed supports GLMs for non-normal outcomes. For example, with a binary outcome:

# Binary outcome example
data_bin <- data
data_bin$Y_binary <- rbinom(n, 1, plogis(0.4 * data$M + 0.2 * data$X))

result_bin <- pmed(
  Y_binary ~ X + M + C,
  formula_m = M ~ X + C,
  data = data_bin,
  family_y = binomial(),    # Logistic regression for binary outcome
  treatment = "X",
  mediator = "M",
  method = "parametric_bootstrap",
  n_boot = 1000
)
print(result_bin)

Bootstrap Methods Comparison

When to use each method:

# Compare methods
result_plugin <- pmed(Y ~ X + M, formula_m = M ~ X, data = data,
                      treatment = "X", mediator = "M", method = "plugin")

result_param <- pmed(Y ~ X + M, formula_m = M ~ X, data = data,
                     treatment = "X", mediator = "M",
                     method = "parametric_bootstrap", n_boot = 1000)

result_nonparam <- pmed(Y ~ X + M, formula_m = M ~ X, data = data,
                        treatment = "X", mediator = "M",
                        method = "nonparametric_bootstrap", n_boot = 1000)

References

Contributing

We welcome contributions! Please see the Issue Tracker for potential features or bug reports.