Package 'medfit'

Description

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,
  cluster = 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

• 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

• Percentile confidence interval bounds (the $\alpha/2$ and $1 - \alpha/2$ quantiles of the bootstrap distribution)

• Bootstrap distribution (for parametric and nonparametric)

• Method used

See Also

BootstrapResult, MediationData, fit_mediation()

Examples

# Fit mediation model
med_data <- fit_mediation(
  formula_y = outcome ~ treatment + mediator1 + covariate1 + covariate2,
  formula_m = mediator1 ~ treatment + covariate1 + covariate2,
  data = mediation_demo,
  treatment = "treatment",
  mediator = "mediator1"
)

# Define indirect effect function
indirect_fn <- function(theta) theta["m_treatment"] * theta["y_mediator1"]

# 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(mediator1 ~ treatment + covariate1 + covariate2,
              data = boot_data)
  fit_y <- lm(outcome ~ treatment + mediator1 + covariate1 + covariate2,
              data = boot_data)
  unname(coef(fit_m)["treatment"] * coef(fit_y)["mediator1"])
}

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



if (requireNamespace("lme4", quietly = TRUE)) {
# Treatment assigned to whole clusters (needs lme4); formulas hold fixed
# effects only, the engine adds the random cluster intercept
set.seed(1)
J <- 30
id <- rep(seq_len(J), each = 6)
cdat <- data.frame(school = factor(id), X = sample(rep(0:1, J / 2))[id])
cdat$M <- 0.5 * cdat$X + rnorm(J, sd = 0.5)[id] + rnorm(J * 6)
cdat$Y <- 0.2 * cdat$X + 0.4 * cdat$M + rnorm(J, sd = 0.5)[id] + rnorm(J * 6)
# Resample whole clusters. A refit that is singular or does not converge is
# dropped, and the warning counts how many were
nie_of <- function(d) {
  unname(nie(fit_mediation(Y ~ X + M, M ~ X, data = d, treatment = "X",
                           mediator = "M", engine = "lmer",
                           cluster = "school")))
}
result_cluster <- bootstrap_mediation(
  nie_of, method = "nonparametric", data = cdat, cluster = "school",
  n_boot = 20, seed = 1
)
print(result_cluster)
}

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
)

Description

S7 class for mediation with a treatment assigned to whole clusters (such as schools) and a mediator and outcome measured on the individuals in them. It holds the paths of a two-level linear mixed model pair: the mediator model M ~ X + ... + (1 | cluster) and the outcome model with the mediator split into its within-cluster deviation and the observed cluster mean. medfit computes the effects; causal interpretation is the user's responsibility.

Usage

ClusterMediationData(a_path, b_within, b_between, c_prime, estimates, vcov,
  kr_df, treatment, mediator, outcome, cluster, n_obs, n_clusters,
  cluster_sizes, parameterization, covariates_centered, se_type, reml,
  converged, sigma_m, sigma_y, tau_m, tau_y, data, source_package)

Arguments

Details

With the within parameterization $Y_{ij} = \theta_0 + c' X_j + b_W (M_{ij} - \bar M_j) + b_B \bar M_j + \dots$ and $a$ the effect of the treatment on the mediator, the natural indirect effect is $a b_B$. It splits into an own-mediator part $a b_W$ and a spillover part $a (b_B - b_W)$.

The split is the large-cluster, cluster-average approximation. With the observed cluster mean, moving only member $i$'s mediator by $a$ also moves the mean by $a / n_j$, so the exact own-mediator effect is $a [b_W + (b_B - b_W) / n_j]$ and the exact spillover is $a (b_B - b_W)(n_j - 1) / n_j$; they still sum to $a b_B$. decompose() warns when the approximation error exceeds half the own effect's standard error.

The estimand assumes randomized treatment, intact clusters, no interference between clusters, linear models without mediator-by-treatment or mediator-by-covariate products, and interference through the observed cluster mean only; all members of a cluster must be in the analysis rows. The own effect also needs no unmeasured lower-level mediator-outcome confounding. The spillover, NIE and NDE need no unmeasured upper-level mediator-outcome confounding either; under it only the sum of the direct effect and $a (b_B - b_W)$ is identified. print() and summary() print this block. The validator ties the alias rows a, c_prime, b_within and b_between of estimates and vcov to the path properties, so an object with inconsistent numbers cannot be built.

kr_df holds the Kenward-Roger degrees of freedom of each path when se_type = "kr", and is empty otherwise.

Value

A ClusterMediationData S7 object.

References

Talloen, W., Moerkerke, B., Loeys, T., De Naeghel, J., Van Keer, H., & Vansteelandt, S. (2016). Estimation of indirect effects in the presence of unmeasured confounding for the mediator-outcome relationship in a multilevel 2-1-1 mediation model. Journal of Educational and Behavioral Statistics, 41(4), 359-391. doi:10.3102/1076998616636855

VanderWeele, T. J. (2010). Direct and indirect effects for neighborhood-based clustered and longitudinal data. Sociological Methods & Research, 38(4), 515-544. doi:10.1177/0049124110366236

Examples

# Hand-built object (the paths and their alias rows must agree)
est <- c(a = 0.5, c_prime = 0.2, b_within = 0.3, b_between = 0.6)
vc <- diag(0.01, 4)
dimnames(vc) <- list(names(est), names(est))
cmd <- ClusterMediationData(
  a_path = 0.5, b_within = 0.3, b_between = 0.6, c_prime = 0.2,
  estimates = est, vcov = vc,
  treatment = "X", mediator = "M", outcome = "Y", cluster = "school",
  n_obs = 400L, n_clusters = 40L, cluster_sizes = rep(10L, 40),
  parameterization = "within", covariates_centered = TRUE, se_type = "model",
  reml = TRUE, converged = TRUE, sigma_m = 1, sigma_y = 1, tau_m = 0.5,
  tau_y = 0.5, source_package = "medfit"
)
cmd@a_path * cmd@b_between  # the NIE, 0.5 * 0.6

Description

Return the components of the total effect for an object whose outcome model has treatment-by-mediator products.

For an InteractionMediationData object this is VanderWeele's (2014) four-way decomposition: controlled direct effect (CDE), reference interaction (INTref), mediated interaction (INTmed), and pure indirect effect (PIE), together with the derived natural direct and indirect effects and the total effect. For a JointMediationData object it is the CDE and the joint natural direct and indirect effects (VanderWeele and Vansteelandt 2014).

Usage

decompose(x, ...)

Arguments

Details

For a 0/1 treatment, outcome model $Y = \theta_0 + \theta_1 X + \theta_2 M + \theta_3 X M + \dots$, mediator model $M = \beta_0 + \beta_1 X + \dots$, and reference mediator level $m^*$:

$CDE = \theta_1 + \theta_3 m^*$

$INTref = \theta_3 (E[M \mid X = 0, \bar c] - m^*)$

$INTmed = \theta_3 \beta_1$

$PIE = \theta_2 \beta_1$

with $NDE = CDE + INTref$, $NIE = INTmed + PIE$, and $TE = NDE + NIE$. $E[M \mid X = 0, \bar c]$ is the mediator model's prediction at no treatment and the covariate means.

Value

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

References

VanderWeele, T. J. (2014). A unification of mediation and interaction: A 4-way decomposition. Epidemiology, 25(5), 749–761.

VanderWeele, T. J., & Vansteelandt, S. (2014). Mediation analysis with multiple mediators. Epidemiologic Methods, 2(1), 95–115. doi:10.1515/em-2012-0010

See Also

nie(), nde(), te()

Examples

fit_m <- lm(mediator1 ~ treatment + covariate1 + covariate2,
            data = mediation_demo)
fit_y <- lm(outcome ~ treatment * mediator1 + covariate1 + covariate2,
            data = mediation_demo)
med_int <- extract_mediation(fit_m, model_y = fit_y,
                             treatment = "treatment", mediator = "mediator1")
decompose(med_int)

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, and a native SEMFit).

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 two linear mixed models for a treatment assigned to whole clusters, as a ClusterMediationData: the mediator model is object, the outcome model is model_y, and cluster names the cluster variable. Both models must be fitted to the same rows; the cluster means must be computed on those rows, and a pair of fits that kept different individuals is an error (best effort: it compares the data row names the fits kept, so it cannot see a mismatch between frames whose row names were both reset)

• SEMFit: Extract from a model fitted by fit_sem() (the experimental native maximum-likelihood engine), with the same arguments and structures as the lavaan method. A treatment-by-mediator product written x:m (or m:x) in the outcome equation gives the four-way decomposition, stored under the canonical name m:x; without it decomposition = "four_way" is refused by name. ⁠interaction =⁠ takes that product written with :.

• brmsfit: Extract from Bayesian models (future)

Note: OpenMx extraction is planned for a future release.

The returned objects share one interface (nie(), nde(), te(), confint(), bootstrap_mediation(), ...) for use by other medfit functions and dependent packages (probmed, RMediation, medrobust).

The class returned depends on the structure. A single mediator gives a MediationData object, or an InteractionMediationData object when the outcome model has a treatment-by-mediator term. A mediator vector of length two or more gives a SerialMediationData or ParallelMediationData object; for lm/glm fits whose outcome model has a treatment-by-mediator term written with : or *, it gives a JointMediationData object with the joint natural effects of the mediators. Any other product term in a multi-mediator fit errors.

Arguments for lm and glm models

With object the fitted mediator model (the first mediator's model for several mediators), the lm/glm method takes:

• model_y: the fitted outcome model. With several mediators include every mediator, not only the last: omitting one that also affects the outcome biases the others' coefficients.

• treatment: name of the treatment variable.

• mediator: name of the mediator, or an ordered character vector of two or more mediator names.

• mediator_models: list of the fitted models for mediators 2 to k, in order; required whenever mediator has length two or more.

• structure: "auto" (default), "serial", or "parallel"; "auto" classifies the mediators from the mediator models.

• decomposition: "auto" (default) uses the four-way decomposition when a single mediator's outcome model has a treatment-by-mediator term; "four_way" requires that term; "two_way" ignores it and returns a MediationData object.

• m_star: the reference mediator level for the controlled direct effect (default 0; for JointMediationData a scalar or a vector named by the interacting mediators).

• vcov_fun: function returning each model's covariance matrix (default stats::vcov(); for example sandwich::vcovHC). Not supported for JointMediationData.

• outcome, data: optional; detected from the models when omitted.

The lavaan method takes a fitted lavaan model and the variable names; see ?extract_mediation_lavaan for its arguments.

Covariance of the estimates

For lm/glm fits, each equation's block of ⁠@vcov⁠ is that model's own covariance. The blocks for different equations are stored as zero. This is exact for simple mediation, serial chains, and the four-way model, whose later equations contain every regressor of the earlier ones, but for parallel mediators it leaves out the covariance between their a paths. JointMediationData stores the full stacked least-squares covariance, and a lavaan fit estimates all equations jointly, so for the same data the intervals can differ between engines.

Coefficient names

The simple single-mediator MediationData from lm/glm fits names three coefficients ⁠m_<treatment>⁠, ⁠y_<mediator>⁠ and ⁠y_<treatment>⁠ (m_X, y_M and y_X for variables X and M), and code that looks coefficients up by name, such as the parametric bootstrap in probmed, relies on them. The lavaan method adds the same three names to ⁠@estimates⁠ and ⁠@vcov⁠ for an observed, continuous treatment, mediator and outcome. Each is a second name for the a, b or c_prime path, with the same estimate and covariances, appended after the existing rows. A lavaan fit with a latent or ordered treatment, mediator or outcome, and any extraction with standardized = TRUE, gets none of the three and works with the plugin estimators only. Serial, parallel and four-way objects are not affected.

Value

An S7 object whose class depends on the mediation structure (see Details): MediationData, InteractionMediationData, SerialMediationData, ParallelMediationData, JointMediationData, or ClusterMediationData. Each carries the path coefficients, the full parameter vector and its variance-covariance matrix, residual standard deviations (for Gaussian models), variable names, and the original data when available.

See Also

MediationData, fit_mediation(), bootstrap_mediation()

Examples

# Extract the mediation structure from fitted lm models, using the
# simulated mediation_demo data bundled with medfit
fit_m <- lm(mediator1 ~ treatment + covariate1 + covariate2,
            data = mediation_demo)
fit_y <- lm(outcome ~ treatment + mediator1 + covariate1 + covariate2,
            data = mediation_demo)
med_data <- extract_mediation(fit_m, model_y = fit_y,
                              treatment = "treatment", mediator = "mediator1")



if (requireNamespace("lme4", quietly = TRUE)) {
# Treatment assigned to whole clusters (needs lme4); formulas hold fixed
# effects only, the engine adds the random cluster intercept
set.seed(1)
J <- 30
id <- rep(seq_len(J), each = 6)
cdat <- data.frame(school = factor(id), X = sample(rep(0:1, J / 2))[id])
cdat$M <- 0.5 * cdat$X + rnorm(J, sd = 0.5)[id] + rnorm(J * 6)
cdat$Y <- 0.2 * cdat$X + 0.4 * cdat$M + rnorm(J, sd = 0.5)[id] + rnorm(J * 6)
# The outcome model carries the mediator as its within-cluster deviation
# plus the observed cluster mean
cdat$M_bar <- ave(cdat$M, cdat$school)
cdat$M_w <- cdat$M - cdat$M_bar
fit_m <- lme4::lmer(M ~ X + (1 | school), data = cdat)
fit_y <- lme4::lmer(Y ~ X + M_w + M_bar + (1 | school), data = cdat)
cluster_med <- extract_mediation(fit_m, model_y = fit_y, treatment = "X",
                                 mediator = "M", cluster = "school")
nie(cluster_med)
}

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.

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", "kr"),
  engine_args = list(),
  m_star = 0,
  cluster = NULL,
  model = NULL,
  ...
)

Arguments

Details

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.)

• With a non-identity link the path coefficients are on the link scale, where the product $a b$ is not a natural indirect effect on the outcome scale; interpret it with care or use engine = "regmedint" for closed-form effects of a logistic outcome

• A treatment-by-mediator term in formula_y requires Gaussian identity-link models (the four-way formulas are linear); other families error

• 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

regmedint Engine

When engine = "regmedint":

• Delegates to regmedint::regmedint() (suggested package) for closed-form natural (in)direct effects, with or without a treatment-mediator interaction

• Returns MediationData or InteractionMediationData depending on whether formula_y contains a treatment-by-mediator interaction term; use engine_args to override the derived regmedint arguments

Reference Mediator Level

When the fit yields an InteractionMediationData, m_star fixes the level $m^*$ at which the controlled direct effect is read off: $CDE = \theta_1 + \theta_3 m^*$ and $INTref = \theta_3 (E[M \mid X = 0] - m^*)$. The two shift in compensating directions, so nde(), nie(), te(), and pm() are invariant to m_star; only the CDE/INTref split moves.

The engines reach the same result by different routes. engine = "glm" applies m_star at extraction time, after the coefficients are fit; engine = "regmedint" passes it to regmedint::regmedint() as m_cde, where it is consumed by that package's closed-form estimator at fitting time. Supplying m_star for a fit with no treatment-by-mediator term is an error, not a silent no-op.

Value

A MediationData object containing the fitted mediation structure, or an InteractionMediationData object when formula_y contains a treatment-by-mediator interaction term.

References

VanderWeele, T. J. (2014). A unification of mediation and interaction: A 4-way decomposition. Epidemiology, 25(5), 749–761.

See Also

MediationData, extract_mediation(), bootstrap_mediation()

Examples

# mediation_demo is simulated data bundled with medfit; its covariates
# confound the mediator-outcome relation, so both models adjust for them
med_data <- fit_mediation(
  formula_y = outcome ~ treatment + mediator1 + covariate1 + covariate2,
  formula_m = mediator1 ~ treatment + covariate1 + covariate2,
  data = mediation_demo,
  treatment = "treatment",
  mediator = "mediator1"
)
print(med_data)


# Binary outcome (takes longer to fit): dichotomize the outcome at its median
demo <- mediation_demo
demo$outcome_bin <- as.integer(demo$outcome > stats::median(demo$outcome))
med_data_bin <- fit_mediation(
  formula_y = outcome_bin ~ treatment + mediator1 + covariate1 + covariate2,
  formula_m = mediator1 ~ treatment + covariate1 + covariate2,
  data = demo,
  treatment = "treatment",
  mediator = "mediator1",
  family_y = binomial()
)



if (requireNamespace("lme4", quietly = TRUE)) {
# Treatment assigned to whole clusters (needs lme4); formulas hold fixed
# effects only, the engine adds the random cluster intercept
set.seed(1)
J <- 30
id <- rep(seq_len(J), each = 6)
cdat <- data.frame(school = factor(id), X = sample(rep(0:1, J / 2))[id])
cdat$M <- 0.5 * cdat$X + rnorm(J, sd = 0.5)[id] + rnorm(J * 6)
cdat$Y <- 0.2 * cdat$X + 0.4 * cdat$M + rnorm(J, sd = 0.5)[id] + rnorm(J * 6)
cluster_fit <- fit_mediation(
  Y ~ X + M, M ~ X, data = cdat,
  treatment = "X", mediator = "M",
  engine = "lmer", cluster = "school"
)
nie(cluster_fit)       # a * b_between
decompose(cluster_fit) # own-mediator and spillover parts
}

Description

Fits a model written in lavaan-style syntax to raw data with medfit's native maximum-likelihood engine, and returns a SEMFit object. The engine handles observed and latent variables, user labels, fixed values, starts and bounds, equality and inequality constraints (==, <, >, linear or nonlinear), and defined parameters (⁠:=⁠) with delta-method standard errors.

Usage

fit_sem(
  model,
  data,
  information = c("observed", "expected"),
  n_starts = NULL,
  control = list(),
  weights = NULL,
  se_type = c("model", "sandwich"),
  meanstructure = NULL,
  fixed_x = FALSE,
  sampling_weights = lifecycle::deprecated()
)

Arguments

Value

A SEMFit object. With a product term whose factor is a regression target (x:m with m regressed on x), glance() reports the chi-square, df and p-value as NA with a note, and anova() and mbco_sem() refuse; see SEMFit-methods.

See Also

SEMFit, fit_mediation()

Examples

fit <- fit_sem(
  "mediator1 ~ a*treatment\noutcome ~ b*mediator1 + treatment\nab := a*b",
  data = mediation_demo
)
fit@defined

# A treatment-by-mediator product: the four-way decomposition
fit4 <- fit_sem(
  "mediator1 ~ treatment\noutcome_int ~ treatment + mediator1 + treatment:mediator1",
  data = mediation_demo, se_type = "sandwich"
)
decompose(extract_mediation(fit4, treatment = "treatment", mediator = "mediator1",
                            outcome = "outcome_int"))

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^*$, $INTref = \theta_3 (E[M \mid X = 0, \bar c] - 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$).

$E[M \mid X = 0, \bar c] = \beta_0 + \gamma^\top \bar c$ is the mediator model's prediction at no treatment and the covariate means $\bar c$. The effects are for the treatment contrast 0 to 1, and the lm/glm extractor requires Gaussian identity-link models. From fit_mediation(engine = "regmedint"), ⁠@estimates⁠ and ⁠@vcov⁠ also carry regmedint's own delta-method covariance of the four components and the derived effects.

Value

An InteractionMediationData S7 object.

References

VanderWeele, T. J. (2014). A unification of mediation and interaction: A 4-way decomposition. Epidemiology, 25(5), 749–761.

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

Description

Recomputes the joint CDE, NDE, NIE and total effect of a JointMediationData object from a named parameter vector, holding the object's structure, reference levels (m_star) and sample covariate means fixed. With the default estimates = object@estimates it returns the stored effects. Its main use is as the statistic of a parametric bootstrap: the draws bootstrap_mediation() passes to statistic_fn are named like ⁠@estimates⁠, and the NDE needs the prefixed intercept and covariate rows (m1_, ..., y_) and the covariate means, not only the path aliases.

Usage

joint_effects(object, estimates = object@estimates)

Arguments

Value

Named numeric vector: cde, nde, nie, te.

Examples

d <- mediation_demo
fit <- extract_mediation(
  lm(mediator1 ~ treatment + covariate1 + covariate2, d),
  model_y = lm(outcome_int ~ treatment * mediator1 + mediator2 +
                 covariate1 + covariate2, d),
  treatment = "treatment", mediator = c("mediator1", "mediator2"),
  mediator_models = list(lm(mediator2 ~ treatment + mediator1 +
                              covariate1 + covariate2, d))
)
joint_effects(fit)

# Parametric bootstrap of the joint NIE
boot <- bootstrap_mediation(
  function(theta) joint_effects(fit, theta)[["nie"]],
  method = "parametric", mediation_data = fit, n_boot = 500, seed = 1
)
boot@ci_lower
boot@ci_upper

Description

S7 class for the joint natural effects of two or more mediators, serial or parallel, when the outcome model may contain treatment-by-mediator product terms (VanderWeele and Vansteelandt 2014). The mediators are treated as one block: the natural indirect effect (NIE) runs through every mediator and every path among them, and the natural direct effect (NDE) is the remainder. medfit computes the effects; causal interpretation is the user's responsibility.

Usage

JointMediationData(structure, mediators, treatment, outcome, interactions,
  a_total, b_paths, theta3, c_prime, cde, nde, nie, total_effect, m_star,
  estimates, vcov, data, n_obs, converged, source_package)

Arguments

Details

With outcome model $Y = \theta_1 X + \sum_i (\theta_{2i} + \theta_{3i} X) M_i + \dots$ and $\beta^*_{1i}$ the total (propagated) effect of the treatment on mediator $i$ (for a serial chain this includes the paths through earlier mediators), the unit-contrast effects are

$NIE = \sum_i (\theta_{2i} + \theta_{3i}) \beta^*_{1i}$

$NDE = \theta_1 + \sum_i \theta_{3i} \mu^*_{0i}$

$CDE = \theta_1 + \sum_i \theta_{3i} m^*_i$

and $TE = NDE + NIE$. $\theta_{3i}$ is zero for a mediator without a product term, and $\mu^*_{0i}$ is the mean of mediator $i$ under no treatment at the covariate means (propagated down a serial chain like $\beta^*_{1i}$). The validator enforces the NIE and CDE identities, so an object with inconsistent numbers cannot be built.

What the joint NIE is, and is not

The NIE is the effect through the mediators as a block: every path from the treatment through any mediator, including paths among the mediators. It is not split into per-mediator or per-path pieces, and none of its parts should be reported as the effect "through M1". For a serial chain it therefore differs from the default nie() of a SerialMediationData (a * d * b, the effect through the full chain only): with no product term, the joint NIE of M1 -> M2 is $a_1 b_1 + (a_2 + d a_1) b_2$, not $a_1 d b_2$. The effects use the unit contrast of a 0/1 treatment (0 to 1).

Assumptions

The joint effects are identified under no-unmeasured-confounding assumptions stated for the whole mediator vector: none for the treatment and outcome, none for the mediators and the outcome, none for the treatment and the mediators, and no mediator-outcome confounder affected by the treatment outside the mediator vector (VanderWeele and Vansteelandt 2014). In a serial chain the earlier mediators are affected by the treatment and confound the later ones; that is allowed precisely because they are part of the vector. Mediator-outcome and treatment-mediator confounders must be controlled for every mediator.

medfit's requirements

• Every model (each mediator model and the outcome model) must carry the same covariates. This is a medfit limitation that keeps the serial-chain algebra and the covariance exact, not a requirement of the method; differing sets error, naming the terms.

• The models must be Gaussian with the identity link, unweighted, with an intercept, fit to the same rows; every mediator must appear in the outcome model; mediator lists the mediators in causal order.

• A product must be written in the outcome formula with : or *. A product precomputed as a data column (for example XM <- X * M1 added to the data and then used as Y ~ X + M1 + M2 + XM) cannot be recognized from the formula: medfit treats it as an ordinary covariate and the effects ignore the interaction, with no error. Write X * M1 in the formula instead.

Covariates and standard errors

The NDE depends on the covariates; it is evaluated at their sample means. The effects are linear in the covariates, so this equals the sample average of the per-observation effects. Standard errors use the delta method with analytic gradients over a stacked-OLS covariance that includes the correlation between parallel mediator equations. They treat the covariate means as fixed, so they are conditional on the observed covariates and slightly understate the uncertainty of a population-average effect. For a parametric bootstrap use joint_effects() as the statistic.

Value

A JointMediationData S7 object.

References

VanderWeele, T. J., & Vansteelandt, S. (2014). Mediation analysis with multiple mediators. Epidemiologic Methods, 2(1), 95–115. doi:10.1515/em-2012-0010

Examples

# Hand-built parallel object with an X x M2 product (m2* = 0)
jmd <- JointMediationData(
  structure = "parallel", mediators = c("M1", "M2"),
  treatment = "X", outcome = "Y", interactions = "M2",
  a_total = c(M1 = 0.5, M2 = 0.3), b_paths = c(M1 = 0.3, M2 = 0.4),
  theta3 = c(M2 = 0.2), c_prime = 0.1,
  cde = 0.1, nde = 0.12, nie = 0.33, total_effect = 0.45,
  m_star = c(M2 = 0),
  estimates = c(a1 = 0.5, a2 = 0.3, b1 = 0.3, b2 = 0.4,
                theta3_M2 = 0.2, c_prime = 0.1),
  vcov = diag(0.01, 6),
  n_obs = 200L, converged = TRUE, source_package = "medfit"
)
jmd@nie  # (0.3 + 0) * 0.5 + (0.4 + 0.2) * 0.3 = 0.33

Description

Tests the null hypothesis that an indirect effect, the product of two or more labeled paths ($a b$, or a serial $a d b$), is zero in a model fitted by fit_sem(). The product is zero exactly when one factor is, so the model is refitted once per path with that path set to zero; each is a linear equality constraint, so each refit is a single solve. The statistic is the smallest of the likelihood-ratio statistics, referred to a chi-square distribution with one degree of freedom (the minimum-of-k LRTs; see the Methods and Formulas article). The nonlinear constraint $a b = 0$ is never solved, because a nonlinear solver can stop at the branch with the larger statistic.

Usage

mbco_sem(
  object,
  a,
  b,
  n_starts = NULL,
  control = NULL,
  fit = lifecycle::deprecated(),
  paths = NULL
)

Arguments

Details

A fit made with weights is refused: the weighted discrepancy is not a likelihood, so its difference is not chi-square. So is a fit that did not converge, since its discrepancy is the reference both refits are measured against. A refit that fails the acceptance gate of the optimizer is an error rather than a statistic.

Inequality constraints in the model matter only when they put the null on the boundary of the parameter space. If a bound or an inequality that involves the tested parameter is active at its refit (a > 0 with the null a = 0), the function warns that the statistic is a mixture of chi-squares, not chi-square(1), and that the p-value is not valid. A bound the null sits inside (a < 0.2), an inactive inequality and one on another parameter are not mentioned. A null the constraints rule out (a > 0.3 with the null a = 0) is an error that says so.

Only a product of labeled paths is tested; a sum of products (several parallel indirect effects together) is not. If two of the paths are tied by an equality (a == b), their refits impose the same restriction and their diff_ll values coincide. The test is conservative: under a null where only one factor is zero its size is at most the nominal level, and with more than one zero factor it is smaller still.

Value

An object of class c("mbco_sem", "htest"): statistic (the smallest diffLL), parameter (df = 1), p.value, method, data.name, n_obs, and tests, a data frame with one row per refit (constraint, the minimized discrepancy f, diff_ll equal to n_obs * (f - fit@f), and converged, always TRUE because a refit that does not converge is an error). It prints with the standard htest layout.

See Also

fit_sem(), SEMFit

Examples

fit <- fit_sem(
  "mediator1 ~ a*treatment\noutcome ~ b*mediator1 + treatment",
  data = mediation_demo
)
mbco_sem(fit, "a", "b")

# a serial indirect effect, the product of three paths
serial <- fit_sem(
  paste0("mediator1 ~ a*treatment\nmediator2 ~ d*mediator1 + treatment\n",
         "outcome ~ b*mediator2 + treatment"),
  data = mediation_demo
)
mbco_sem(serial, paths = c("a", "d", "b"))

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

# mediation_demo is simulated data bundled with medfit; its covariates
# confound the mediator-outcome relation, so adjust for them
result <- med(
  data = mediation_demo,
  treatment = "treatment",
  mediator = "mediator1",
  outcome = "outcome",
  covariates = c("covariate1", "covariate2")
)
print(result)

# Quick summary
quick(result)


# With bootstrap CI (slower)
result_boot <- med(
  data = mediation_demo,
  treatment = "treatment",
  mediator = "mediator1",
  outcome = "outcome",
  covariates = c("covariate1", "covariate2"),
  boot = TRUE,
  n_boot = 1000,
  seed = 42
)

Description

The model syntax read by fit_sem() and by fit_mediation(engine = "native"). It is a lavaan-style core plus two extensions, each marked "medfit extension" below. One statement per line (or separated by ⁠;⁠).

Core statements

Constraints may be nonlinear in the labels (cp == a*b, a*b == 0.1, a*b > 0.1, ab > 0 with ab := a*b), solved with the Jacobian recomputed at every point; a model with one tries every start and keeps the best accepted solution, warning when the starts disagree. A product of labels set to zero (a*b == 0, b1*b2*b3 == 0) is instead solved exactly: one linear fit per factor (a == 0, b == 0), keeping the smallest discrepancy. A sum of products (a1*b1 + a2*b2 == 0) has no exact reduction, so it is best-effort: the fit is the best of the starts tried, a disagreement warning means a better optimum may exist, and n_starts raises the number of starts. Labels used in a constraint or a definition must be defined in the model. A constraint with no free parameter is an error. Exogenous variables are given free variances and covariances, as lavaan does with fixed.x = FALSE; fixed_x = TRUE (or a vector of names) in fit_sem() fixes the moments of the selected ones at the sample values instead, as fixed.x = TRUE does. Unlike lavaan, which keeps an explicit x ~~ x or x ~ 1 free under fixed.x = TRUE and lets the degrees of freedom move, medfit stops with an error when the model has a statement on a selected variable's moments. Any intercept statement (y ~ 1, or 1 as a term: y ~ 1 + x) turns the mean structure on and frees every observed intercept the model does not write, as lavaan does; y ~ 0*1 then restricts the mean of y to 0 (one degree of freedom), and meanstructure = TRUE in fit_sem() turns it on without any statement, and a product term (x:m, see Product terms) always turns it on. A latent variable's mean is fixed at 0; a free latent ~1 is an error. A repeated 1 for one variable is an error. A variable cannot appear on both sides of ~ or ⁠=~⁠ (y ~ y), and a path, a (co)variance or an intercept takes one row: a repeated y ~ x is an error, not a second row that competes for the label.

Product terms

x:m is the product of two observed variables, written with no spaces, on the right side of a regression (y ~ x + m + x:m). As in lavaan, the product is an observed exogenous variable: the engine builds the column x * m from data after listwise deletion, with a free variance, free covariances with the other exogenous variables and (the mean structure is on) a free mean. The likelihood does not change. x:m and m:x are one term, stored and printed under the canonical name, the two names sorted in the C locale (m:x); the coefficient is named y~m:x, or by its label (b*x:m labels the whole product). extract_mediation() finds the treatment-by-mediator product through the same name and returns the four-way decomposition; centering the mediator moves the split between the controlled direct effect and the reference interaction, not the natural effects.

A column in data named m:x or x:m is used only if it equals x * m on the rows kept; otherwise the fit stops. A product term turns the mean structure on, and meanstructure = FALSE is then an error. A coefficient equals the equation-wise least-squares estimate inside the class below, and the sandwich covariance (se_type = "sandwich") is the standard-error choice, because a product is not normal; the engine warns once when it is not used. fixed_x includes the product columns.

Rules 8 and 9 are checked on the parsed table, so they need no data. Rule 8 rejects a regressor built from the outcome (it is correlated with the equation's error); rule 9 names what is not yet supported with products, and the same model works without the product. The covariance that completion would add between two pure outcomes is fixed at 0 in a product model, with one note. Products of observed variables only: a product with a latent variable, x:x and three-way products are not supported.

The global chi-square is not valid when a factor of a product is itself a regression target (the model treats the product as exogenous although it is a function of that variable) or when the exogenous (co)variance block is restricted. glance() then reports chi-square, df and p-value as NA with a note, anova() and mbco_sem() refuse, and expected information is refused for an endogenous product; logLik(), AIC and BIC stay defined. A product of exogenous variables with the block free keeps the ordinary test.

medfit extensions

• Comma shorthand (medfit extension). A comma-separated left side expands to one statement per name: ⁠y1, y2 ~ x⁠ means y1 ~ x and y2 ~ x. Commas on the right side are not supported; use +. lavaan reads ⁠y1, y2 ~ x⁠ differently: it silently keeps only y2. A model copied between the two packages can therefore change meaning without an error.

• CONSTRAINT(...) (medfit extension). A constraint line may be written CONSTRAINT(a == b), equivalent to the bare a == b. The keyword is case-sensitive and starts the statement, one pair of parentheses encloses the whole constraint, it holds exactly one comparison, and no text follows the closing parenthesis.

Expressions

Both sides of a constraint and the right side of a definition are arithmetic expressions in the labels, with ⁠+ - * / ^⁠, parentheses, finite numbers, and the functions exp, log, log10, sqrt, abs, min, max, sin, cos, tan, pnorm and qnorm. Anything else is rejected with the construct named.

See Also

fit_sem(), SEMFit

Examples

# Labeled paths and a defined parameter
fit_sem(
  "mediator1 ~ a*treatment\noutcome ~ b*mediator1 + cp*treatment\nab := a*b",
  data = mediation_demo
)@defined

# Comma shorthand: one left-hand side per name
coef(fit_sem(
  "mediator1, mediator2 ~ treatment\noutcome ~ mediator1 + mediator2 + treatment",
  data = mediation_demo
))[1:3]

# An equality constraint, written with the CONSTRAINT() keyword
coef(fit_sem(
  "mediator1 ~ a*treatment\noutcome ~ b*mediator1 + treatment\nCONSTRAINT(a == b)",
  data = mediation_demo
))[c("a", "b")]

Description

A simulated dataset for demonstrating medfit's simple, serial, parallel, and treatment-by-mediator interaction workflows with one running example. The data are simulated: they describe no real study, population, or finding.

Usage

mediation_demo

Format

A data frame with 400 rows and 8 variables:

treatment

Integer, 0/1. Randomized treatment assignment.

mediator1

Numeric. First mediator; affected by treatment.

mediator2

Numeric. Serial mediator; affected by treatment and mediator1.

mediator3

Numeric. Parallel mediator; affected by treatment only.

covariate1

Numeric. Continuous confounder of the mediator-outcome relations.

covariate2

Integer, 0/1. Binary confounder of the mediator-outcome relations.

outcome

Numeric. Outcome with no product terms; use it for the simple, serial, and parallel models.

outcome_int

Numeric. outcome plus a treatment by mediator1 interaction; use it only for models with that interaction: the four-way decomposition with a single mediator, or joint effects with several mediators (JointMediationData).

Details

The variables follow a fictional workplace-training story, used only so the examples read concretely: randomized assignment to a training program (treatment) raises skill confidence (mediator1), which builds task mastery (mediator2); the assignment also triggers supervisor check-ins (mediator3); outcome is job performance.

All errors are independent standard normal. The generating equations are

$M_1 = 0.5X + 0.3C_1 + 0.3C_2 + e_1$

$M_2 = 0.2X + 0.5M_1 + 0.2C_1 + e_2$

$M_3 = 0.5X + e_3$

$Y = 0.2X + 0.4M_1 + 0.3M_2 + 0.3M_3 + 0.3C_1 + 0.2C_2 + e_4$

and outcome_int adds $0.5 X M_1$.

Because the covariates confound the mediator-outcome relations, every model should adjust for both covariate1 and covariate2. A model that omits a downstream mediator estimates reduced-form coefficients rather than the structural ones above: for example, the simple model outcome ~ treatment + mediator1 + covariates has limiting coefficients 0.55 for mediator1 and 0.41 for treatment. For a serial model, include mediator1 in the outcome model as well; the serial indirect effect $a \times d \times b$ is then the effect through the chain mediator1 to mediator2 only.

Source

Simulated for package demonstration; not drawn from or representing any real study. The generating script is data-raw/mediation_demo.R in the package source repository.

Examples

data(mediation_demo)
str(mediation_demo)

# Simple mediation, adjusting for both covariates
fit <- fit_mediation(
  formula_y = outcome ~ treatment + mediator1 + covariate1 + covariate2,
  formula_m = mediator1 ~ treatment + covariate1 + covariate2,
  data = mediation_demo,
  treatment = "treatment",
  mediator = "mediator1"
)
nie(fit)

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"
)

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(s).

Usage

nde(x, ...)

Arguments

Details

Without a treatment-by-mediator product (MediationData, SerialMediationData, ParallelMediationData) the NDE is the direct-path coefficient, $NDE = c'$. With a product it also depends on the mediator's mean under no treatment:

• InteractionMediationData: $NDE = CDE + INTref = \theta_1 + \theta_3 E[M \mid X = 0, \bar c]$.

• JointMediationData: $NDE = \theta_1 + \sum_i \theta_{3i} \mu^*_{0i}$, with $\mu^*_{0i}$ the mean of mediator $i$ under no treatment at the covariate means.

Value

A numeric scalar of class mediation_effect carrying a type attribute.

See Also

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

Examples

med_data <- fit_mediation(
  formula_y = outcome ~ treatment + mediator1 + covariate1 + covariate2,
  formula_m = mediator1 ~ treatment + covariate1 + covariate2,
  data = mediation_demo,
  treatment = "treatment",
  mediator = "mediator1"
)

nde(med_data)

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(s).

Usage

nie(x, ...)

Arguments

Details

The effects are for a unit contrast of the treatment. By class:

• MediationData: $NIE = a b$.

• SerialMediationData: with type = "chain" (the default), the effect through the full chain only, $NIE = a \, d_1 \cdots d_{k-1} \, b$; with type = "total", the total indirect effect, the sum over every treatment-to-outcome path through at least one mediator (including paths that skip a mediator, such as X -> M2 -> Y), which equals te(x) - nde(x).

• ParallelMediationData: $NIE = \sum_j a_j b_j$.

• InteractionMediationData: $NIE = INTmed + PIE = (\theta_2 + \theta_3) \beta_1$.

• JointMediationData: the joint NIE through all the mediators, $\sum_i (\theta_{2i} + \theta_{3i}) \beta^*_{1i}$.

The "Methods and Formulas" article on the package website gives the full formulas.

Value

A numeric scalar of class mediation_effect carrying a type attribute. For intervals use confint() or bootstrap_mediation().

See Also

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

Examples

med_data <- fit_mediation(
  formula_y = outcome ~ treatment + mediator1 + covariate1 + covariate2,
  formula_m = mediator1 ~ treatment + covariate1 + covariate2,
  data = mediation_demo,
  treatment = "treatment",
  mediator = "mediator1"
)

nie(med_data)

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

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

• d (two mediators) or d21, d32, ...: Mediator-to-mediator paths

• b: Last mediator -> Outcome

• c_prime: Direct effect

For parallel mediation (ParallelMediationData): a1, b1, a2, b2, ..., c_prime.

For InteractionMediationData: a, b, c_prime, and theta3 (the treatment-by-mediator coefficient).

For JointMediationData: the raw coefficients a1..aK, dij, b1..bK, ⁠theta3_<mediator>⁠, and c_prime. The a paths here are the raw coefficients, not the propagated ⁠a*⁠ values used by the effects.

Value

A named numeric vector of path coefficients

See Also

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

Examples

med_data <- fit_mediation(
  formula_y = outcome ~ treatment + mediator1 + covariate1 + covariate2,
  formula_m = mediator1 ~ treatment + covariate1 + covariate2,
  data = mediation_demo,
  treatment = "treatment",
  mediator = "mediator1"
)

paths(med_data)

Description

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

Usage

pm(x, ...)

Arguments

Details

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

For serial mediation (SerialMediationData) the numerator is the total indirect effect, nie(x, type = "total"), and the denominator the full total effect from te().

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

The ratio has no delta-method standard error in medfit; bootstrap it with bootstrap_mediation().

Value

A numeric scalar of class mediation_effect, usually between 0 and 1 (negative or greater than 1 in cases of suppression effects), or NA with a warning when the total effect is numerically zero.

See Also

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

Examples

med_data <- fit_mediation(
  formula_y = outcome ~ treatment + mediator1 + covariate1 + covariate2,
  formula_m = mediator1 ~ treatment + covariate1 + covariate2,
  data = mediation_demo,
  treatment = "treatment",
  mediator = "mediator1"
)

pm(med_data)

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

result <- med(
  data = mediation_demo,
  treatment = "treatment",
  mediator = "mediator1",
  outcome = "outcome",
  covariates = c("covariate1", "covariate2")
)

# One-line summary
quick(result)

Description

S7 class holding a structural equation model fitted by the native maximum-likelihood engine (see fit_sem()): estimates, their covariance, the parameter table, defined parameters (⁠:=⁠) with delta-method standard errors, and the diagnostics of the optimizer.

Usage

SEMFit(theta, vcov, table, defined, f, n_obs, n_dropped, df, information,
  converged, data, model, diagnostics, internals, call)

Arguments

Value

A SEMFit S7 object

See Also

fit_sem()

Examples

fit <- fit_sem(
  "mediator1 ~ a*treatment\noutcome ~ b*mediator1 + treatment\nab := a*b",
  data = mediation_demo
)
fit@table[fit@table$op == "~", c("lhs", "rhs", "est", "se")]
fit@defined
fit@converged

Description

Accessors and display for a SEMFit returned by fit_sem().

Arguments

Details

• coef() returns the free-parameter estimates, named as lavaan names its coefficients: the label when the row has one, else ⁠lhs op rhs⁠ without spaces ("M~X", "X~~X", "eta=~m2"). vcov(), tidy() and the accessors extract_mediation() reads use the same names.

• vcov() returns their covariance, projected onto the active constraints and bounds; rows for a pinned parameter are zero.

• nobs() returns the number of rows used after listwise deletion.

• logLik() returns the log likelihood with the means free; its df attribute counts the free parameters plus one mean per observed variable that is not fixed. With fixed_x it is conditional on the fixed variables (lavaan's value), so AIC() and BIC() of fits with different fixed sets or data stop instead of ranking them; the same holds for anova() and mbco_sem(). A fixed_x fit also lists its fixed moments after the free parameters in tidy() (fixed_x column TRUE, no standard error) and in summary(); coef() and vcov() leave them out. print() and summary() say what the log likelihood is conditional on.

• Model test with a product term. The global chi-square is not valid when a factor of a product term (x:m) is itself a regression target, because the model then treats the product as exogenous although it is a function of m; nor when the variances and covariances of the exogenous variables (every product included) are restricted (fixed, labeled, bounded, or left out because a ⁠~~⁠ statement names one of them; a fixed_x moment is not a restriction). fit_sem() then warns once, glance() reports chisq, df and p_value as NA (with the reason in attr(, "note")), summary() prints the reason, and anova() and mbco_sem() stop naming the product. The coefficients equal equation-wise least squares and the sandwich standard errors are unaffected, but only with observed information: fit_sem(information = "expected") stops for a product with an endogenous factor, because its standard errors (the sandwich built on them included) are wrong there. fit_sem() also warns once for any product term fitted without se_type = "sandwich". logLik(), AIC() and BIC() stay defined.

• anova() tests nested fits by the likelihood-ratio statistic $n (F_r - F_f)$ on the difference in model degrees of freedom: anova(fit_full, fit_restricted), or several fits, each tested against the one with fewer degrees of freedom. It compares fits that were already made and never refits. The fits must use the same observed variables and rows, and neither may have sampling weights or have failed to converge. Nesting is checked by parameter-table inclusion: the restricted fit may free no parameter the other leaves out and must have more degrees of freedom. That is necessary for nesting, not sufficient; whether two models are nested in substance is the analyst's responsibility. A restricted fit with a smaller discrepancy than the model it restricts is an error (one fit is not at its global optimum). When a bound or an inequality on a tested parameter is active at the restricted fit, the null lies on the boundary of the parameter space, the statistic is a mixture of chi-squares, and the function warns that the p-value is not valid, as mbco_sem() does. Returns an object of class "anova" with columns Df, Chisq (n times the discrepancy), ⁠Chisq diff⁠, ⁠Df diff⁠ and ⁠Pr(>Chisq)⁠. Name rows by naming the arguments after the first.

• confint() returns Wald (normal) intervals for the free parameters and the defined parameters, confint(fit, parm, level = 0.95); parm takes names or positions. A parameter pinned by a constraint or sitting on a bound has a zero standard error and a point interval.

• tidy() returns one row per free parameter, then the moments fixed_x fixed (a fixed_x column, only when something is fixed), then per defined parameter: term (the name coef() uses), op (the operator of its row in the parameter table, ⁠:=⁠ for a defined parameter), estimate, std_error, statistic (the z value, NA when the standard error is zero), p_value (two-sided), and with conf_int = TRUE the interval confint() returns (conf_low, conf_high, at conf_level). type is "all" (default), "parameters" or "defined".

• glance() returns one row: nobs, npar (free parameters after the equality constraints), df, chisq (n times the minimized discrepancy, the test of the model against the saturated one), p_value (NA when df is 0), logLik, AIC and BIC (with the means free, as logLik() defines them) and fixed_x, the names the fit holds at their sample moments, comma-joined ("" when none). A fit with sampling weights has no likelihood, so its chisq, p_value, logLik, AIC and BIC are NA. There is no converged column: fit_sem() returns a fit only for a solution that passed the acceptance gate and errors otherwise, so ⁠@converged⁠ is TRUE for every SEMFit and would carry no information.

• summary() collects the estimates, standard errors and z statistics, the defined parameters, and the diagnostics: rows dropped, information type, degrees of freedom, retries, stationarity, active bounds and constraints, improper-solution flags and the multiplier sign check.

Value

coef() a named numeric vector; vcov() a matrix; nobs() an integer; logLik() an object of class "logLik"; summary() an object of class "summary_SEMFit"; print() invisibly returns its argument.

See Also

fit_sem(), SEMFit

Examples

fit <- fit_sem(
  "mediator1 ~ a*treatment\noutcome ~ b*mediator1 + treatment\nab := a*b",
  data = mediation_demo
)
coef(fit)
vcov(fit)[c("a", "b"), c("a", "b")]
nobs(fit)
summary(fit)
generics::tidy(fit, conf_int = TRUE)
generics::glance(fit)

# a likelihood-ratio test of b = 0 against the fitted model
fit0 <- fit_sem(paste0(fit@model, "\nb == 0"), data = mediation_demo)
anova(fit, fit0)

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 indirect effect through the full chain is the product of the chain's 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

This is the effect through the chain only, and what nie() returns by default. Paths that skip a mediator (X -> M2, M1 -> Y, ...) are not part of it; nie(x, type = "total") adds them, and te() is the full total effect over every path, so pm() is the total indirect share of it. With treatment-by-mediator products in the outcome model, extract_mediation() returns a JointMediationData object with the effect through all the mediators together.

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"
)

Description

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

Usage

te(x, ...)

Arguments

Details

$TE = NIE + NDE$

For a SerialMediationData object the total effect is the sum over every directed X-to-Y path in the fitted models: the direct path, the full chain, and every path that skips a mediator (for two mediators, $c' + a_1 b_1 + a_2 b_2 + a_1 d_{21} b_2$). With the same covariates in every equation and linear models this equals the treatment coefficient of the outcome regressed on the treatment and covariates alone. A path missing from its model counts as zero; when a path is in a model but its coefficient was not recorded (a hand-built object), te() returns NA with a warning. For glm fits with a non-identity link the sum of path products is on the linear-predictor scale, as for MediationData.

Value

A numeric scalar of class mediation_effect carrying a type attribute.

See Also

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

Examples

med_data <- fit_mediation(
  formula_y = outcome ~ treatment + mediator1 + covariate1 + covariate2,
  formula_m = mediator1 ~ treatment + covariate1 + covariate2,
  data = mediation_demo,
  treatment = "treatment",
  mediator = "mediator1"
)

te(med_data)

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

Description

tidy() converts a mediation data object or a BootstrapResult into a tidy tibble, one row per path coefficient or effect. glance() returns a one-row summary. The base generics stats::coef(), stats::vcov(), stats::confint(), and stats::nobs() also have methods for every mediation class; see Details.

Usage

tidy_s7_object(x, ...)

glance_s7_object(x, ...)

Arguments

Details

Path standard errors are the square roots of the diagonal of ⁠@vcov⁠. Effect standard errors (NIE, NDE, TE, and for interaction objects the four-way components) use the delta method over the full ⁠@vcov⁠, the same computation as confint(parm = "effects"), so tidy(conf_int = TRUE) reproduces confint() exactly. Intervals are normal approximations ($\hat{\theta} \pm z \, SE$); the sampling distribution of a product of coefficients is skewed, so for inference on indirect effects prefer bootstrap_mediation(). tidy() raises no warning about this; confint() does.

For a serial chain fitted as separate lm/glm regressions, ⁠@vcov⁠ has zero covariances between equations, and the effect standard errors inherit that. The proportion mediated (reported by glance()) has no standard error: it is a ratio whose delta-method standard error is unstable when the total effect is near zero, so bootstrap it instead.

Base methods

• coef(object, type = "paths"): the path coefficients; type = "effects" gives the effects, "all" both, and for interaction objects "components" gives the four-way components.

• vcov(object): the stored ⁠@vcov⁠, covering every entry of ⁠@estimates⁠.

• confint(object, parm = "paths", level = 0.95): normal intervals for the paths, or with parm = "effects" for the effects using the delta-method standard errors above (interaction objects also accept parm = "components"). With parm = "effects" it warns that the normal approximation may be inaccurate for an indirect effect.

• nobs(object): the number of observations.

Value

tidy(): a tibble (a data frame if tibble is not installed) with columns term, estimate, std_error, and, when conf_int = TRUE, conf_low and conf_high.

glance(): a one-row tibble with nie, nde, te, pm, nobs, and converged; interaction objects add interaction and m_star, and joint objects add cde, structure, n_mediators, interactions, and m_star.

See Also

bootstrap_mediation(), MediationData

Examples

med_data <- fit_mediation(
  formula_y = outcome ~ treatment + mediator1 + covariate1 + covariate2,
  formula_m = mediator1 ~ treatment + covariate1 + covariate2,
  data = mediation_demo,
  treatment = "treatment",
  mediator = "mediator1"
)
# tidy() is the generic from the generics package (also re-exported by broom)
generics::tidy(med_data)
generics::tidy(med_data, type = "effects", conf_int = TRUE)