Package {medfit}

Type:	Package
Title:	Infrastructure for Mediation Model Fitting and Extraction
Version:	0.3.2
Date:	2026-07-20
Description:	Provides S7-based infrastructure for fitting mediation models, extracting path coefficients, and performing bootstrap inference. Designed as a foundation package for the mediation analysis ecosystem, supporting 'probmed', 'RMediation', and 'medrobust' packages. Implements unified interfaces for model fitting across different engines (currently generalized linear models, with future support for mixed models and Bayesian methods), standardized extraction of mediation paths from various model types, and robust bootstrap inference methods. Mediation methods are described in MacKinnon, Lockwood and Williams (2004) <doi:10.1207/s15327906mbr3901_4>, Preacher and Hayes (2008) <doi:10.3758/brm.40.3.879>, Tofighi and MacKinnon (2011) <doi:10.3758/s13428-011-0076-x>, and VanderWeele (2014) <doi:10.1097/EDE.0000000000000121>.
License:	GPL (≥ 3)
URL:	https://data-wise.github.io/medfit/, https://github.com/data-wise/medfit
BugReports:	https://github.com/data-wise/medfit/issues
Depends:	R (≥ 4.1.0)
Imports:	S7 (≥ 0.1.0), stats, methods, checkmate, generics, MASS
Suggests:	lavaan (≥ 0.6-0), sandwich, testthat (≥ 3.0.0), tibble
Encoding:	UTF-8
Language:	en-US
Config/testthat/edition:	3
Config/Needs/website:	knitr, rmarkdown, quarto
Config/roxygen2/version:	8.0.0
NeedsCompilation:	no
Packaged:	2026-07-23 16:49:54 UTC; dt
Author:	Davood Tofighi [aut, cre, cph]
Maintainer:	Davood Tofighi <dtofighi@gmail.com>
Repository:	CRAN
Date/Publication:	2026-07-23 17:00:02 UTC

medfit: Infrastructure for Mediation Model Fitting and Extraction

Description

Provides S7-based infrastructure for fitting mediation models, extracting path coefficients, and performing bootstrap inference. Designed as a foundation package for probmed, RMediation, and medrobust.

Details

Key functions:

• fit_mediation: Fit mediation models

• extract_mediation: Extract from fitted models

• bootstrap_mediation: Bootstrap inference

Key classes:

• MediationData: Mediation model structure

• BootstrapResult: Bootstrap results

Author(s)

Maintainer: Davood Tofighi dtofighi@gmail.com (ORCID) [copyright holder]

Authors:

• Davood Tofighi dtofighi@gmail.com (ORCID) [copyright holder]

See Also

Useful links:

• https://data-wise.github.io/medfit/

• https://github.com/data-wise/medfit

• Report bugs at https://github.com/data-wise/medfit/issues

Classify a multi-mediator structure as serial or parallel

Description

Conservative, backward-compatible inference for structure = "auto". Returns "parallel" only on POSITIVE evidence of a parallel structure (no mediator is regressed on another, and every mediator enters the outcome model); otherwise defaults to "serial" (the historical default for vector mediator). It never errors – malformed inputs fall through to the chosen worker's own validation, which emits specific, directed messages. Users can always set structure explicitly to override.

Usage

.classify_multimediator_structure(med_models, mediators, treatment, model_y)

Arguments

Value

"serial" or "parallel".

Classify a multi-mediator lavaan structure as serial or parallel

Description

Conservative, backward-compatible inference for structure = "auto" on lavaan objects – the SEM analogue of .classify_multimediator_structure() for lm/glm. Returns "parallel" only on POSITIVE evidence (no mediator is regressed on another); otherwise defaults to "serial" (the historical default for vector mediator). It never errors – malformed inputs fall through to the chosen worker's own directed validation. Users can always set structure explicitly to override.

Usage

.classify_multimediator_structure_lavaan(
  object,
  mediators,
  standardized = FALSE
)

Arguments

Details

Detection reads only op == "~" (regression) rows, so residual covariances (⁠~~⁠) among mediators cannot masquerade as serial chain edges. Like the lm/glm classifier, it returns "parallel" only on POSITIVE evidence: no mediator-on-mediator edge AND every mediator enters a single common outcome equation. Anything else (e.g. one mediator missing from the outcome model) falls back to "serial", preserving the historical vector-mediator behavior.

Value

"serial" or "parallel".

Expand a Source Covariance Matrix with Full-Copy Path Aliases

Description

Appends named structural aliases (e.g. a, b, c_prime, d1, ...) to a source variance-covariance matrix, copying the FULL covariance row/column of each alias's source parameter rather than just its diagonal variance. This preserves every covariance the aliased parameter has – both with the original parameters and with the other aliases.

Usage

.expand_vcov_with_aliases(vcov_src, source_idx, aliases_to_add)

Arguments

Details

This is the shared engine behind the alias-vcov contract used by the lm/glm and lavaan extract_mediation() methods (simple and serial). Factoring it here keeps the two extractors from drifting: each computes its own source_idx mapping (the lavaan path tries labels then variable names; the lm path maps to the prefixed coefficient names) and then hands the mechanical expansion to this single routine.

Value

A symmetric numeric matrix of dimension nrow(vcov_src) + length(aliases_to_add), with the original block intact, each alias row/column populated from its source, and the alias-to-alias intersections filled from the corresponding source-to-source covariances.

Extract Treatment-Mediator Interaction Structure from a lavaan Model

Description

Internal worker for the four-way (VanderWeele 2014) branch of extract_mediation() on lavaan objects. The SEM analogue of .extract_interaction_mediation_lm(): it returns an InteractionMediationData object for continuous Y and M with binary treatment and reference level m_star.

Usage

.extract_interaction_mediation_lavaan(
  object,
  treatment,
  mediator,
  int_term,
  outcome = NULL,
  m_star = 0,
  standardized = FALSE
)

Arguments

Details

Because lavaan fits one joint system, the expanded vcov preserves the FULL covariance among the paths – including cov(beta1, theta3) and cov(beta0, theta3) – unlike the block-diagonal lm/glm engine, so the delta-method standard errors reflect the joint estimation. The mediator intercept beta0 (needed for INTref) is read from the ~1 row, so the model must be fit with meanstructure = TRUE.

Value

An InteractionMediationData object.

Extract Treatment-Mediator Interaction Structure from lm/glm Models

Description

Internal worker for the four-way (VanderWeele 2014) branch of the lm/glm extract_mediation() method. Invoked by .extract_mediation_lm_impl() when a single mediator's outcome model carries an X:M term. Builds an InteractionMediationData object for continuous Y and M with binary treatment (0 -> 1) and reference mediator level m_star.

Usage

.extract_interaction_mediation_lm(
  model_m,
  model_y,
  treatment,
  mediator,
  int_term,
  outcome = NULL,
  data = NULL,
  m_star = 0
)

Arguments

Details

With mediator model M = \beta_0 + \beta_1 X + \beta_2^\top C and outcome model Y = \theta_0 + \theta_1 X + \theta_2 M + \theta_3 XM + \dots the components are CDE = \theta_1 + \theta_3 m^*, INTref = \theta_3 (E[M\mid X=0] - m^*), INTmed = \theta_3 \beta_1, PIE = \theta_2 \beta_1, where E[M\mid X=0] evaluates covariates at their sample means. The combined vcov is block-diagonal across the two separately-fitted equations and named with the aliases a, b, c_prime, theta3, b0.

Value

An InteractionMediationData object.

Internal Implementation for lm/glm Extraction

Description

Internal Implementation for lm/glm Extraction

Usage

.extract_mediation_lm_impl(
  model_m,
  model_y,
  treatment,
  mediator,
  mediator_models = NULL,
  outcome = NULL,
  data = NULL,
  structure = c("auto", "serial", "parallel"),
  decomposition = c("auto", "four_way", "two_way"),
  m_star = 0,
  vcov_fun = stats::vcov
)

Arguments

Value

MediationData object, or SerialMediationData when mediator is a vector of length >= 2

Extract Parallel Mediation Structure from a lavaan Model

Description

Internal worker for the parallel branch of extract_mediation() on lavaan objects (X -> M_j -> Y for k independent mediators). It is the SEM analogue of .extract_parallel_mediation_lm() and returns a ParallelMediationData object. Total indirect effect = ⁠sum_j a_j * b_j⁠.

Usage

.extract_parallel_mediation_lavaan(
  object,
  treatment,
  mediators,
  outcome = NULL,
  standardized = FALSE,
  ...
)

Arguments

Details

Paths are located in the lavaan parameter table by variable name:

• a_j: M_j ~ X

• b_j: Y ~ M_j

• ⁠c'⁠ : Y ~ X (defaults to 0 with a warning if absent – full mediation)

Unlike the lm/glm engine – where the M_j come from separate regressions so cov(a_j, b_j) = 0 – lavaan estimates the whole system jointly, so the expanded vcov preserves the FULL off-diagonal structure (including cov(a_j, b_j) and ⁠cov(a_j, a_j')⁠). Downstream SEs therefore reflect the true joint covariance; tests must not hardcode any of these to zero.

Value

A ParallelMediationData object.

Extract Parallel Mediation Structure from lm/glm Models

Description

Internal worker for parallel mediation (X -> M_j -> Y, independent mediators). Mirrors .extract_serial_mediation_lm() but the mediator models are NOT chained: mediator_models[[j - 1]] is the model for mediators[j] regressed on the treatment (and covariates), in mediator-index order.

Usage

.extract_parallel_mediation_lm(
  object,
  mediator_models,
  model_y,
  treatment,
  mediators,
  outcome = NULL,
  data = NULL
)

Arguments

Value

A ParallelMediationData object.

Extract Serial Mediation Structure from a lavaan Model

Description

Internal worker for the serial branch of extract_mediation() on lavaan objects. It is invoked by extract_mediation_lavaan() when mediator is a character vector of length >= 2, and returns a SerialMediationData object describing the chain X -> M1 -> M2 -> ... -> Mk -> Y.

Usage

.extract_serial_mediation_lavaan(
  object,
  treatment,
  mediators,
  outcome = NULL,
  standardized = FALSE,
  ...
)

Arguments

Details

Paths are located in the lavaan parameter table by variable name:

• a : M1 ~ X

• d_i: ⁠M_{i+1} ~ M_i⁠ for ⁠i = 1 .. k-1⁠ (the k - 1 inter-mediator paths)

• b : Y ~ Mk

• ⁠c'⁠ : Y ~ X (defaults to 0 with a warning if absent – full mediation)

As in the simple-mediation extractor, named structural aliases (a, d1, ..., ⁠d{k-1}⁠, b, c_prime) are appended to estimates and the variance-covariance matrix is expanded so that the FULL covariance row/column of each source parameter is preserved. This lets downstream code recover the true joint covariance of the chain (including off-diagonals) via, for example, vcov[c("a", "d1", "b"), c("a", "d1", "b")] – which is required for serial indirect-effect standard errors.

Value

A SerialMediationData object.

Extract Serial Mediation Structure from lm/glm Models

Description

Internal worker for the serial branch of the lm/glm extract_mediation() method. Invoked by .extract_mediation_lm_impl() when mediator is a character vector of length >= 2. It assembles a SerialMediationData object for the chain X -> M1 -> M2 -> ... -> Mk -> Y from k + 1 separately fitted regressions.

Usage

.extract_serial_mediation_lm(
  object,
  mediator_models,
  model_y,
  treatment,
  mediators,
  outcome = NULL,
  data = NULL
)

Arguments

Details

Path resolution: a = coefficient of treatment in object; d_i = coefficient of mediators[i] in mediator_models[[i]] (the predecessor mediator, read regardless of any additional covariates in that equation); b = coefficient of mediators[k] in model_y; ⁠c'⁠ = coefficient of treatment in model_y (0 with a warning if absent).

The combined vcov is block-diagonal across the separately-fitted equations (so cov(a, d_i) = cov(d_i, b) = 0) but preserves the within-model_y covariance, so ⁠cov(b, c')⁠ is non-zero. See the extract_mediation lm method docs for the lm-vs-lavaan covariance divergence this implies.

Value

A SerialMediationData object.

Extract Residual Standard Deviation from Model

Description

Extract Residual Standard Deviation from Model

Usage

.extract_sigma(model)

Arguments

Value

Numeric scalar or NULL

Locate a treatment-by-mediator interaction term in an outcome model

Description

Returns the coefficient name of the X:M product term in model_y, trying both orderings (treatment:mediator and mediator:treatment, since the formula order determines which lm()/glm() emits). Returns NA_character_ when no interaction term is present.

Usage

.find_interaction_term(model_y, treatment, mediator)

Locate a treatment-by-mediator interaction term in a lavaan outcome equation

Description

In lavaan the interaction enters as a product predictor of the outcome (a data column, e.g. XM). This returns its coefficient name, preferring an explicit interaction argument and otherwise trying treatment:mediator / mediator:treatment. Returns NA_character_ when none is found.

Usage

.find_interaction_term_lavaan(
  object,
  treatment,
  mediator,
  interaction = NULL,
  standardized = FALSE
)

Extract Response Variable Name from Model

Description

Extract Response Variable Name from Model

Usage

.get_response_var(model)

Arguments

Value

Character string: response variable name

Register lavaan Method for extract_mediation

Description

This function is called from .onLoad() to register the S7 method for lavaan objects when the lavaan package is available.

Usage

.register_lavaan_method()

BootstrapResult S7 Class

Description

S7 class containing results from bootstrap inference, including point estimates, confidence intervals, and bootstrap distribution.

Usage

BootstrapResult(estimate, ci_lower, ci_upper, ci_level, boot_estimates,
  n_boot, method, call)

Arguments

Details

This class standardizes bootstrap inference results across different bootstrap methods (parametric, nonparametric, plugin).

The class includes validation to ensure consistency between method type and required fields.

Value

A BootstrapResult S7 object

Examples

# Parametric bootstrap result
result <- BootstrapResult(
  estimate = 0.15,
  ci_lower = 0.10,
  ci_upper = 0.20,
  ci_level = 0.95,
  boot_estimates = rnorm(1000, 0.15, 0.02),
  n_boot = 1000L,
  method = "parametric",
  call = NULL
)

InteractionMediationData: Mediation with Treatment-Mediator Interaction

Description

S7 class for simple mediation with a treatment-by-mediator interaction (X \rightarrow M \rightarrow Y where the outcome model contains an X \times M term). It carries VanderWeele's (2014) four-way decomposition of the total effect into controlled direct effect (CDE), reference interaction (INTref), mediated interaction (INTmed), and pure indirect effect (PIE):

TE = CDE + INTref + INTmed + PIE

with NDE = CDE + INTref and NIE = INTmed + PIE. medfit computes the decomposition; causal interpretation is the user's responsibility (it requires the four no-unmeasured-confounding assumptions of VanderWeele 2014).

Usage

InteractionMediationData(a_path, b_path, c_prime, interaction, cde,
  int_ref, int_med, pie, nde, nie, total_effect, m_star, estimates, vcov,
  sigma_m, sigma_y, treatment, mediator, outcome, mediator_predictors,
  outcome_predictors, data, n_obs, converged, source_package)

Arguments

Details

Path coefficients follow the outcome model Y = \theta_0 + \theta_1 X + \theta_2 M + \theta_3 XM + \dots and mediator model M = \beta_0 + \beta_1 X + \dots: a_path = \beta_1, b_path = \theta_2, c_prime = \theta_1, interaction = \theta_3. With reference level m_star (m^*) the components are CDE = \theta_1 + \theta_3 m^*, INTmed = \theta_3 \beta_1, and PIE = \theta_2 \beta_1. When \theta_3 = 0 the decomposition collapses to standard simple mediation (CDE = NDE = \theta_1; INTref = INTmed = 0; NIE = PIE = \theta_2\beta_1).

Value

An InteractionMediationData S7 object.

Examples

# Hand-built object (theta3 = 0.2 interaction, m* = 0)
imd <- InteractionMediationData(
  a_path = 0.5, b_path = 0.3, c_prime = 0.1, interaction = 0.2,
  cde = 0.1, int_ref = 0.04, int_med = 0.10, pie = 0.15,
  nde = 0.14, nie = 0.25, total_effect = 0.39, m_star = 0,
  estimates = c(a = 0.5, b = 0.3, c_prime = 0.1, theta3 = 0.2),
  vcov = diag(0.01, 4),
  treatment = "X", mediator = "M", outcome = "Y",
  mediator_predictors = "X", outcome_predictors = c("X", "M", "X:M"),
  n_obs = 200L, converged = TRUE, source_package = "medfit"
)

nie(imd)       # INTmed + PIE = 0.25
decompose(imd) # all four components + derived effects

MediationData S7 Class

Description

S7 class containing standardized mediation model structure, including path coefficients, parameter estimates, variance-covariance matrix, and metadata.

Usage

MediationData(a_path, b_path, c_prime, estimates, vcov, sigma_m, sigma_y,
  family_m, family_y, treatment, mediator, outcome, mediator_predictors,
  outcome_predictors, data, n_obs, converged, source_package)

Arguments

Details

This class provides a unified container for mediation model information extracted from various model types (lm, glm, lavaan, etc.). It ensures consistency across the mediation analysis ecosystem.

The class includes comprehensive validation to ensure data integrity.

Value

A MediationData S7 object

Examples

# Create a MediationData object
med_data <- MediationData(
  a_path = 0.5,
  b_path = 0.3,
  c_prime = 0.2,
  estimates = c(0.5, 0.3, 0.2),
  vcov = diag(3) * 0.01,
  sigma_m = 1.0,
  sigma_y = 1.2,
  treatment = "X",
  mediator = "M",
  outcome = "Y",
  mediator_predictors = "X",
  outcome_predictors = c("X", "M"),
  data = NULL,
  n_obs = 100L,
  converged = TRUE,
  source_package = "stats"
)

ParallelMediationData: Parallel (Multiple-Mediator) Mediation Structure

Description

S7 class for parallel mediation, where a treatment affects an outcome through two or more independent mediators operating in parallel (X \rightarrow M_j \rightarrow Y for j = 1, \dots, k). The total indirect effect is the sum of the per-mediator products, \sum_{j=1}^{k} a_j b_j. This complements MediationData (simple) and SerialMediationData (serial chains).

Usage

ParallelMediationData(a_paths, b_paths, c_prime, estimates, vcov,
  sigma_mediators, sigma_y, treatment, mediators, outcome,
  mediator_predictors, outcome_predictors, data, n_obs, converged,
  source_package)

Arguments

Value

A ParallelMediationData S7 object.

Examples

pmd <- ParallelMediationData(
  a_paths = c(0.5, 0.4),
  b_paths = c(0.6, 0.3),
  c_prime = 0.2,
  estimates = c(0.5, 0.4, 0.6, 0.3, 0.2),
  vcov = diag(0.01, 5),
  treatment = "X",
  mediators = c("M1", "M2"),
  outcome = "Y",
  mediator_predictors = list("X", "X"),
  outcome_predictors = c("X", "M1", "M2"),
  n_obs = 200L,
  converged = TRUE,
  source_package = "medfit"
)

nie(pmd)   # total indirect effect: sum(a_j * b_j) = 0.42
paths(pmd) # a1, b1, a2, b2, c_prime

SerialMediationData S7 Class

Description

S7 class for serial mediation models where the effect flows through multiple mediators in sequence: X -> M1 -> M2 -> ... -> Mk -> Y.

This class supports serial mediation chains of any length, from simple two-mediator models (product-of-three: a * d * b) to complex chains with many mediators (product-of-k).

Usage

SerialMediationData(a_path, d_path, b_path, c_prime, estimates, vcov,
  sigma_mediators, sigma_y, treatment, mediators, outcome,
  mediator_predictors, outcome_predictors, data, n_obs, converged,
  source_package)

Arguments

Details

Serial Mediation Structure

Serial mediation models the indirect effect flowing through a sequence of mediators. The total indirect effect is the product of all path coefficients:

• 2 mediators (product-of-three): Indirect = a * d * b

• 3 mediators (product-of-four): Indirect = a * d21 * d32 * b

• k mediators (product-of-k+1): Indirect = a * d21 * d32 * ... * d(k,k-1) * b

Path Notation

• a: Treatment -> First mediator (X -> M1)

• d21: First -> Second mediator (M1 -> M2)

• d32: Second -> Third mediator (M2 -> M3)

• dji: Previous mediator -> Current mediator

• b: Last mediator -> Outcome (Mk -> Y)

• ⁠c'⁠: Direct effect (X -> Y, controlling for all mediators)

Extensibility

This class is designed to handle serial chains of any length:

• Minimal case: 2 mediators (length(d_path) = 1)

• No upper limit on chain length

• Validator ensures consistency between mediators and paths

Value

A SerialMediationData S7 object

Examples

# Two-mediator serial mediation (X -> M1 -> M2 -> Y)
# Product-of-three: a * d * b
serial_data <- SerialMediationData(
  a_path = 0.5,       # X -> M1
  d_path = 0.4,       # M1 -> M2 (scalar for 2 mediators)
  b_path = 0.3,       # M2 -> Y
  c_prime = 0.1,      # X -> Y (direct)
  estimates = c(0.5, 0.4, 0.3, 0.1),
  vcov = diag(4) * 0.01,
  sigma_mediators = c(1.0, 1.1),  # SD for M1, M2 models
  sigma_y = 1.2,
  treatment = "X",
  mediators = c("M1", "M2"),
  outcome = "Y",
  mediator_predictors = list(
    c("X"),           # M1 ~ X
    c("X", "M1")      # M2 ~ X + M1
  ),
  outcome_predictors = c("X", "M1", "M2"),  # Y ~ X + M1 + M2
  data = NULL,
  n_obs = 100L,
  converged = TRUE,
  source_package = "lavaan"
)

# Three-mediator serial mediation (X -> M1 -> M2 -> M3 -> Y)
# Product-of-four: a * d21 * d32 * b
serial_data_3 <- SerialMediationData(
  a_path = 0.5,           # X -> M1
  d_path = c(0.4, 0.35),  # M1 -> M2, M2 -> M3 (vector for 3 mediators)
  b_path = 0.3,           # M3 -> Y
  c_prime = 0.1,
  estimates = c(0.5, 0.4, 0.35, 0.3, 0.1),
  vcov = diag(5) * 0.01,
  sigma_mediators = c(1.0, 1.1, 1.05),  # SD for M1, M2, M3 models
  sigma_y = 1.2,
  treatment = "X",
  mediators = c("M1", "M2", "M3"),
  outcome = "Y",
  mediator_predictors = list(
    c("X"),              # M1 ~ X
    c("X", "M1"),        # M2 ~ X + M1
    c("X", "M1", "M2")   # M3 ~ X + M1 + M2
  ),
  outcome_predictors = c("X", "M1", "M2", "M3"),
  data = NULL,
  n_obs = 100L,
  converged = TRUE,
  source_package = "lavaan"
)

Perform Bootstrap Inference for Mediation Statistics

Description

Conduct bootstrap inference to compute confidence intervals for mediation statistics. Supports parametric, nonparametric, and plugin methods.

Conduct bootstrap inference to compute confidence intervals for mediation statistics. Supports parametric, nonparametric, and plugin methods.

Usage

bootstrap_mediation(
  statistic_fn,
  method = c("parametric", "nonparametric", "plugin"),
  mediation_data = NULL,
  data = NULL,
  n_boot = 1000L,
  ci_level = 0.95,
  parallel = FALSE,
  ncores = NULL,
  seed = NULL,
  ...
)

bootstrap_mediation(
  statistic_fn,
  method = c("parametric", "nonparametric", "plugin"),
  mediation_data = NULL,
  data = NULL,
  n_boot = 1000L,
  ci_level = 0.95,
  parallel = FALSE,
  ncores = NULL,
  seed = NULL,
  ...
)

Arguments

Details

Bootstrap Methods

Parametric Bootstrap (method = "parametric"):

• Samples parameter vectors from N(\hat{\theta}, \hat{\Sigma})

• Fast and efficient

• Assumes asymptotic normality of parameters

• Recommended for most applications with n > 50

Nonparametric Bootstrap (method = "nonparametric"):

• Resamples observations with replacement

• Refits models for each bootstrap sample

• More robust, no normality assumption

• Computationally intensive

• Use when normality is questionable or n is small

Plugin Estimator (method = "plugin"):

• Computes point estimate only

• No confidence interval

• Fastest method

• Use for quick checks or when CI not needed

Parallel Processing

Set parallel = TRUE to use multiple cores:

• Automatically detects available cores

• Falls back to sequential if parallel fails

• Seed handling ensures reproducibility

Reproducibility

Always set a seed for reproducible results:

bootstrap_mediation(..., seed = 12345)

Bootstrap Methods

Parametric Bootstrap (method = "parametric"):

• Samples parameter vectors from N(\hat{\theta}, \hat{\Sigma})

• Fast and efficient

• Assumes asymptotic normality of parameters

• Recommended for most applications with n > 50

• Requires mediation_data argument

Nonparametric Bootstrap (method = "nonparametric"):

• Resamples observations with replacement

• Refits models for each bootstrap sample

• More robust, no normality assumption

• Computationally intensive

• Use when normality is questionable or n is small

• Requires data argument

Plugin Estimator (method = "plugin"):

• Computes point estimate only

• No confidence interval

• Fastest method

• Use for quick checks or when CI not needed

• Requires mediation_data argument

Statistic Function

The statistic_fn should be a function that:

• For parametric/plugin: Takes a named numeric vector of parameters

• For nonparametric: Takes a data frame

• Returns a single numeric value

Common statistic functions for indirect effect:

# Using parameter names from MediationData
indirect_fn <- function(theta) {
  theta["m_X"] * theta["y_M"]
}

Parallel Processing

Set parallel = TRUE to use multiple cores:

• Uses parallel::mclapply() on Unix systems

• Falls back to sequential on Windows

• Automatically detects available cores

Reproducibility

Always set a seed for reproducible results:

bootstrap_mediation(..., seed = 12345)

Value

A BootstrapResult object containing:

• Point estimate

• Confidence interval bounds

• Bootstrap distribution (for parametric and nonparametric)

• Method used

A BootstrapResult object containing:

• Point estimate

• Confidence interval bounds

• Bootstrap distribution (for parametric and nonparametric)

• Method used

See Also

BootstrapResult, MediationData, extract_mediation()

BootstrapResult, MediationData, fit_mediation()

Examples

## Not run: 
# Parametric bootstrap for indirect effect
result <- bootstrap_mediation(
  statistic_fn = function(theta) theta["a"] * theta["b"],
  method = "parametric",
  mediation_data = med_data,
  n_boot = 5000,
  ci_level = 0.95,
  seed = 12345
)

# Nonparametric bootstrap with parallel processing
result <- bootstrap_mediation(
  statistic_fn = function(data) {
    # Refit models and compute statistic
    # ...
  },
  method = "nonparametric",
  data = mydata,
  n_boot = 5000,
  parallel = TRUE,
  seed = 12345
)

# Plugin estimator (no CI)
result <- bootstrap_mediation(
  statistic_fn = function(theta) theta["a"] * theta["b"],
  method = "plugin",
  mediation_data = med_data
)

## End(Not run)

# Generate example data
set.seed(123)
n <- 100
mydata <- data.frame(X = rnorm(n))
mydata$M <- 0.5 * mydata$X + rnorm(n)
mydata$Y <- 0.3 * mydata$X + 0.4 * mydata$M + rnorm(n)

# Fit mediation model
med_data <- fit_mediation(
  formula_y = Y ~ X + M,
  formula_m = M ~ X,
  data = mydata,
  treatment = "X",
  mediator = "M"
)

# Define indirect effect function
indirect_fn <- function(theta) theta["m_X"] * theta["y_M"]

# Plugin estimator (point estimate only, fastest)
result_plugin <- bootstrap_mediation(
  statistic_fn = indirect_fn,
  method = "plugin",
  mediation_data = med_data
)
print(result_plugin)


# Parametric bootstrap (recommended for most applications)
result <- bootstrap_mediation(
  statistic_fn = indirect_fn,
  method = "parametric",
  mediation_data = med_data,
  n_boot = 1000,
  ci_level = 0.95,
  seed = 12345
)
print(result)

# Nonparametric bootstrap (slower but more robust)
refit_fn <- function(boot_data) {
  fit_m <- lm(M ~ X, data = boot_data)
  fit_y <- lm(Y ~ X + M, data = boot_data)
  unname(coef(fit_m)["X"] * coef(fit_y)["M"])
}

result_np <- bootstrap_mediation(
  statistic_fn = refit_fn,
  method = "nonparametric",
  data = mydata,
  n_boot = 500,
  seed = 12345
)
print(result_np)

Four-Way Decomposition of a Mediation Effect

Description

Return VanderWeele's (2014) four-way decomposition of the total effect for an InteractionMediationData object: controlled direct effect (CDE), reference interaction (INTref), mediated interaction (INTmed), and pure indirect effect (PIE), together with the derived natural direct/indirect and total effects.

Usage

decompose(x, ...)

Arguments

Value

A named numeric vector: c(cde, int_ref, int_med, pie, nde, nie, total).

See Also

nie(), nde(), te()

Extract Mediation Structure from Fitted Models

Description

Generic function to extract mediation structure (a, b, c' paths and variance-covariance matrices) from fitted models. This function provides a unified interface for extracting mediation information from various model types (lm, glm, lavaan, lmer, brms, etc.).

Usage

extract_mediation(object, ...)

Arguments

Details

The extract_mediation() generic provides methods for different model types:

• lm/glm: Extract from linear and generalized linear models

• lavaan: Extract from structural equation models

• lmerMod: Extract from mixed-effects models (future)

• brmsfit: Extract from Bayesian models (future)

Note: OpenMx extraction is planned for a future release.

All methods return a standardized MediationData object that can be used with other medfit functions and dependent packages (probmed, RMediation, medrobust).

Value

A MediationData object containing:

• Path coefficients (a, b, c')

• Full parameter vector and variance-covariance matrix

• Residual variances (for Gaussian models)

• Variable names and metadata

• Original data (if available)

See Also

MediationData, fit_mediation(), bootstrap_mediation()

Examples

# Simulate data with a single mediator (X -> M -> Y)
set.seed(123)
n <- 200
X <- rnorm(n)
M <- 0.5 * X + rnorm(n)
Y <- 0.3 * M + 0.2 * X + rnorm(n)
dat <- data.frame(X = X, M = M, Y = Y)

# Extract the mediation structure from fitted lm models
fit_m <- lm(M ~ X, data = dat)
fit_y <- lm(Y ~ X + M, data = dat)
med_data <- extract_mediation(fit_m, model_y = fit_y,
                              treatment = "X", mediator = "M")

Extract Mediation Structure from lavaan Model

Description

Internal function for extracting mediation structure from lavaan models. This function is registered as an S7 method in .onLoad() when lavaan is available.

Usage

extract_mediation_lavaan(
  object,
  treatment,
  mediator,
  outcome = NULL,
  a_label = "a",
  b_label = "b",
  cp_label = "cp",
  standardized = FALSE,
  structure = c("auto", "serial", "parallel"),
  decomposition = c("auto", "four_way", "two_way"),
  interaction = NULL,
  m_star = 0,
  ...
)

Arguments

Details

This method extracts mediation structure from a fitted lavaan SEM model. The lavaan model should specify labeled paths for the mediation structure.

Typical lavaan Model Specification

model <- "
  # Mediator model
  M ~ a*X

  # Outcome model
  Y ~ b*M + cp*X

  # Indirect and total effects (optional)
  indirect := a*b
  total := cp + a*b
"

Path Labels

By default, the function looks for paths labeled:

• a: Treatment -> Mediator path

• b: Mediator -> Outcome path

• cp: Treatment -> Outcome (direct effect) path

You can customize these labels using the a_label, b_label, and cp_label arguments.

Alternative: Unlabeled Paths

If paths are not labeled, the function will attempt to identify them by variable names. This requires specifying treatment, mediator, and outcome arguments.

Value

A MediationData object; a SerialMediationData object when mediator is a length >= 2 vector resolving to a serial chain; or a ParallelMediationData object when it resolves to parallel mediation.

Examples

if (requireNamespace("lavaan", quietly = TRUE)) {
  # Simulate a simple mediation data set (X -> M -> Y)
  set.seed(123)
  n <- 200
  X <- rnorm(n)
  M <- 0.5 * X + rnorm(n)
  Y <- 0.3 * M + 0.2 * X + rnorm(n)
  dat <- data.frame(X = X, M = M, Y = Y)

  # Fit a labeled lavaan mediation model
  model <- "
    M ~ a*X
    Y ~ b*M + cp*X
  "
  fit <- lavaan::sem(model, data = dat)

  # Extract the mediation structure (dispatches to the lavaan method)
  med_data <- extract_mediation(
    fit,
    treatment = "X",
    mediator = "M",
    outcome = "Y"
  )
}

Fit Mediation Models

Description

Fit mediation models using a specified modeling engine. This function provides a convenient formula-based interface for fitting both the mediator and outcome models simultaneously.

Fit mediation models using a specified modeling engine. This function provides a convenient formula-based interface for fitting both the mediator and outcome models simultaneously.

Usage

fit_mediation(
  formula_y,
  formula_m,
  data,
  treatment,
  mediator,
  engine = "glm",
  family_y = stats::gaussian(),
  family_m = stats::gaussian(),
  weights = NULL,
  se_type = c("model", "sandwich"),
  ...
)

fit_mediation(
  formula_y,
  formula_m,
  data,
  treatment,
  mediator,
  engine = "glm",
  family_y = stats::gaussian(),
  family_m = stats::gaussian(),
  weights = NULL,
  se_type = c("model", "sandwich"),
  ...
)

Arguments

Details

The fit_mediation() function fits both the mediator model and outcome model using the specified engine, then extracts the mediation structure using extract_mediation().

Supported Engines

GLM (engine = "glm"):

• Fits models using stats::glm()

• Supports all GLM families (gaussian, binomial, poisson, etc.)

• For Gaussian models, extracts residual variances

Future Engines:

• "lmer": Mixed-effects models via lme4

• "brms": Bayesian models via brms

Model Specification

The formulas should follow standard R formula syntax:

• formula_m: Mediator model (e.g., M ~ X + C1 + C2)

• formula_y: Outcome model (e.g., Y ~ X + M + C1 + C2)

The mediator must appear in formula_y, and the treatment must appear in both formulas.

Model Specification

The function fits two models:

1. Mediator model: formula_m (e.g., M ~ X + C1 + C2)

2. Outcome model: formula_y (e.g., Y ~ X + M + C1 + C2)

The treatment variable must appear in both formulas. The mediator variable must appear in the outcome formula but NOT in the mediator formula (as it is the response).

GLM Engine

When engine = "glm" (default):

• Models are fit using stats::glm()

• Supports all GLM families (gaussian, binomial, poisson, etc.)

• For Gaussian models, residual standard deviations are extracted

• Non-Gaussian outcomes have sigma_y = NULL

Common Family Specifications

• gaussian(): Continuous outcomes (default)

• binomial(): Binary outcomes

• poisson(): Count outcomes

• Gamma(): Positive continuous outcomes

Value

A MediationData object containing the fitted mediation structure

A MediationData object containing the fitted mediation structure

See Also

MediationData, extract_mediation(), bootstrap_mediation()

Examples

## Not run: 
# Fit Gaussian mediation model
med_data <- fit_mediation(
  formula_y = Y ~ X + M + C,
  formula_m = M ~ X + C,
  data = mydata,
  treatment = "X",
  mediator = "M",
  engine = "glm"
)

# Fit with binary outcome
med_data <- fit_mediation(
  formula_y = Y ~ X + M + C,
  formula_m = M ~ X + C,
  data = mydata,
  treatment = "X",
  mediator = "M",
  engine = "glm",
  family_y = binomial()
)

## End(Not run)

# Generate example data
set.seed(123)
n <- 100
mydata <- data.frame(
  X = rnorm(n),
  C = rnorm(n)
)
mydata$M <- 0.5 * mydata$X + 0.2 * mydata$C + rnorm(n)
mydata$Y <- 0.3 * mydata$X + 0.4 * mydata$M + 0.1 * mydata$C + rnorm(n)

# Simple mediation with continuous variables
med_data <- fit_mediation(
  formula_y = Y ~ X + M,
  formula_m = M ~ X,
  data = mydata,
  treatment = "X",
  mediator = "M"
)
print(med_data)

# With covariates
med_data_cov <- fit_mediation(
  formula_y = Y ~ X + M + C,
  formula_m = M ~ X + C,
  data = mydata,
  treatment = "X",
  mediator = "M"
)


# Binary outcome (takes longer to fit)
mydata$Y_bin <- rbinom(n, 1, plogis(0.3 * mydata$X + 0.4 * mydata$M))
med_data_bin <- fit_mediation(
  formula_y = Y_bin ~ X + M,
  formula_m = M ~ X,
  data = mydata,
  treatment = "X",
  mediator = "M",
  family_y = binomial()
)

Simple Mediation Analysis

Description

A simplified entry point for mediation analysis. Specify the data and variable names, and get results with minimal configuration.

This is the recommended starting point for most mediation analyses. For more control over model specifications, use fit_mediation() directly.

Usage

med(
  data,
  treatment,
  mediator,
  outcome,
  covariates = NULL,
  boot = FALSE,
  n_boot = 1000L,
  seed = NULL,
  ...
)

Arguments

Details

med() is designed to be the simplest way to run a mediation analysis. It constructs the model formulas automatically from variable names.

Default Behavior

• Fits Gaussian (continuous) mediator and outcome models

• No covariates unless specified

• No bootstrap unless requested (use boot = TRUE)

Accessing Results

After running med(), use:

• nie(result): Natural indirect effect

• nde(result): Natural direct effect

• te(result): Total effect

• pm(result): Proportion mediated

• quick(result): One-line summary

• summary(result): Detailed summary

Value

A MediationData object with mediation results

See Also

fit_mediation() for full control, quick() for instant summary, nie(), nde(), te(), pm() for extracting effects

Examples

# Generate example data
set.seed(123)
n <- 200
mydata <- data.frame(
  treatment = rnorm(n),
  covariate = rnorm(n)
)
mydata$mediator <- 0.5 * mydata$treatment + 0.2 * mydata$covariate + rnorm(n)
mydata$outcome <- 0.3 * mydata$treatment + 0.4 * mydata$mediator +
                  0.1 * mydata$covariate + rnorm(n)

# Simple mediation (no covariates)
result <- med(
  data = mydata,
  treatment = "treatment",
  mediator = "mediator",
  outcome = "outcome"
)
print(result)

# With covariates
result_cov <- med(
  data = mydata,
  treatment = "treatment",
  mediator = "mediator",
  outcome = "outcome",
  covariates = "covariate"
)

# Quick summary
quick(result)


# With bootstrap CI (slower)
result_boot <- med(
  data = mydata,
  treatment = "treatment",
  mediator = "mediator",
  outcome = "outcome",
  boot = TRUE,
  n_boot = 1000,
  seed = 42
)

Extract Natural Direct Effect (NDE)

Description

Extract the natural direct effect from a mediation analysis result. The NDE represents the effect of treatment on outcome that does NOT operate through the mediator.

Usage

nde(x, ...)

Arguments

Details

For both simple and serial mediation:

NDE = c'

where c' is the direct effect coefficient.

Value

A numeric value with optional attributes for confidence intervals

See Also

nie(), te(), pm(), paths()

Examples

# Generate example data
set.seed(123)
n <- 100
mydata <- data.frame(X = rnorm(n))
mydata$M <- 0.5 * mydata$X + rnorm(n)
mydata$Y <- 0.3 * mydata$X + 0.4 * mydata$M + rnorm(n)

med_data <- fit_mediation(
  formula_y = Y ~ X + M,
  formula_m = M ~ X,
  data = mydata,
  treatment = "X",
  mediator = "M"
)

nde(med_data)

Extract Natural Indirect Effect (NIE)

Description

Extract the natural indirect effect from a mediation analysis result. The NIE represents the effect of treatment on outcome that operates through the mediator.

Usage

nie(x, ...)

Arguments

Details

For simple mediation (MediationData):

NIE = a \times b

For serial mediation (SerialMediationData):

NIE = a \times d_{21} \times d_{32} \times \ldots \times b

Value

A numeric value (or named vector for SerialMediationData) with optional attributes for confidence intervals if available

See Also

nde(), te(), pm(), paths()

Examples

# Generate example data
set.seed(123)
n <- 100
mydata <- data.frame(X = rnorm(n))
mydata$M <- 0.5 * mydata$X + rnorm(n)
mydata$Y <- 0.3 * mydata$X + 0.4 * mydata$M + rnorm(n)

med_data <- fit_mediation(
  formula_y = Y ~ X + M,
  formula_m = M ~ X,
  data = mydata,
  treatment = "X",
  mediator = "M"
)

nie(med_data)

Extract All Path Coefficients

Description

Extract all path coefficients from a mediation analysis result.

Usage

paths(x, ...)

Arguments

Details

For simple mediation (MediationData):

• a: Treatment -> Mediator (X -> M)

• b: Mediator -> Outcome (M -> Y | X)

• c_prime: Direct effect (X -> Y | M)

For serial mediation (SerialMediationData):

• a: Treatment -> First mediator

• d21, d32, ...: Mediator-to-mediator paths

• b: Last mediator -> Outcome

• c_prime: Direct effect

Value

A named numeric vector of path coefficients

See Also

nie(), nde(), te(), pm()

Examples

# Generate example data
set.seed(123)
n <- 100
mydata <- data.frame(X = rnorm(n))
mydata$M <- 0.5 * mydata$X + rnorm(n)
mydata$Y <- 0.3 * mydata$X + 0.4 * mydata$M + rnorm(n)

med_data <- fit_mediation(
  formula_y = Y ~ X + M,
  formula_m = M ~ X,
  data = mydata,
  treatment = "X",
  mediator = "M"
)

paths(med_data)

Extract Proportion Mediated (PM)

Description

Extract the proportion of the total effect that is mediated (operates through the mediator).

Usage

pm(x, ...)

Arguments

Details

PM = \frac{NIE}{TE} = \frac{NIE}{NIE + NDE}

The proportion mediated can be:

• Between 0 and 1: Normal mediation

• Greater than 1: Suppression (direct and indirect effects have opposite signs)

• Negative: Inconsistent mediation

Value

A numeric value between 0 and 1 (or negative/greater than 1 in cases of suppression effects)

See Also

nie(), nde(), te(), paths()

Examples

# Generate example data
set.seed(123)
n <- 100
mydata <- data.frame(X = rnorm(n))
mydata$M <- 0.5 * mydata$X + rnorm(n)
mydata$Y <- 0.3 * mydata$X + 0.4 * mydata$M + rnorm(n)

med_data <- fit_mediation(
  formula_y = Y ~ X + M,
  formula_m = M ~ X,
  data = mydata,
  treatment = "X",
  mediator = "M"
)

pm(med_data)

Print Method for mediation_effect

Description

Print Method for mediation_effect

Usage

## S3 method for class 'mediation_effect'
print(x, ...)

Arguments

Value

Invisibly returns x (the mediation_effect object). Called for its side effect of printing a formatted effect summary to the console.

Print Summary for BootstrapResult

Description

Print Summary for BootstrapResult

Usage

## S3 method for class 'summary.BootstrapResult'
print(x, ...)

Arguments

Value

Invisibly returns x (the summary.BootstrapResult object). Called for its side effect of printing the formatted summary to the console.

Print Summary for MediationData

Description

Print Summary for MediationData

Usage

## S3 method for class 'summary.MediationData'
print(x, ...)

Arguments

Value

Invisibly returns x (the summary.MediationData object). Called for its side effect of printing the formatted summary to the console.

Print Summary for SerialMediationData

Description

Print Summary for SerialMediationData

Usage

## S3 method for class 'summary.SerialMediationData'
print(x, ...)

Arguments

Value

Invisibly returns x (the summary.SerialMediationData object). Called for its side effect of printing the formatted summary to the console.

Quick Summary of Mediation Results

Description

Print a one-line summary of mediation results, perfect for quick checks or ADHD-friendly workflows.

Usage

quick(x, digits = 3, ...)

Arguments

Details

Prints a compact one-line summary showing:

• NIE (Natural Indirect Effect) with CI if available

• NDE (Natural Direct Effect)

• Proportion Mediated (PM)

If bootstrap results are available (from med(..., boot = TRUE)), confidence intervals are shown for NIE.

Value

Invisibly returns x

See Also

med(), nie(), nde(), pm()

Examples

# Generate example data
set.seed(123)
n <- 100
mydata <- data.frame(X = rnorm(n))
mydata$M <- 0.5 * mydata$X + rnorm(n)
mydata$Y <- 0.3 * mydata$X + 0.4 * mydata$M + rnorm(n)

result <- med(
  data = mydata,
  treatment = "X",
  mediator = "M",
  outcome = "Y"
)

# One-line summary
quick(result)

Extract Total Effect (TE)

Description

Extract the total effect from a mediation analysis result. The TE is the sum of the indirect and direct effects.

Usage

te(x, ...)

Arguments

Details

TE = NIE + NDE

Value

A numeric value with optional attributes for confidence intervals

See Also

nie(), nde(), pm(), paths()

Examples

# Generate example data
set.seed(123)
n <- 100
mydata <- data.frame(X = rnorm(n))
mydata$M <- 0.5 * mydata$X + rnorm(n)
mydata$Y <- 0.3 * mydata$X + 0.4 * mydata$M + rnorm(n)

med_data <- fit_mediation(
  formula_y = Y ~ X + M,
  formula_m = M ~ X,
  data = mydata,
  treatment = "X",
  mediator = "M"
)

te(med_data)

# Verify: TE = NIE + NDE
nie(med_data) + nde(med_data)