Package 'iconic'

Description

Reference-named entry point: applies UNADJ, DIRECT, COCA, IV2SLS, and PGC (matrix bridge) mediation estimators to each tested outcome feature and returns tidy per-feature results with significance flags for both NDE and NIE. When the dataset includes a mediator instrument (Gm), the 2-stage MR estimator fit_iv2sls_mediation2 (method "IV2SLS2") is also run. When the dataset includes path-specific negative controls (W1, W2;), the two-stage proximal estimator fit_pgc_mediation2 is run as "PGC2" (without Gm) and "PGC2Gm" (with Gm, when Gm is also present).

Usage

analyze_mediation_robust(
  iteration_data,
  test_features = NULL,
  alpha = 0.05,
  n_cores = 1,
  se_method = c("delta", "bootstrap", "composite"),
  n_boot = 500
)

Arguments

Details

Negative controls are always summarised over the full W panel for COCA, so restricting test_features does not change its control summary. The matrix-bridge PGC uses the full W panel.

A scalar-bridge variant (fit_pgc_scalar_mediation()) is exported for standalone use but is not included in the default pipeline.

Value

Data frame: feature, method, NDE, NDE_se, NDE_p, NIE, NIE_se, NIE_p, NDE_significant, NIE_significant.

Examples

dat <- generate_toy_data(n = 300, mo_confounding = 0.8, seed = 1)
analyze_mediation_robust(dat)

Parallel version of analyze_methods_robust()

Description

Distributes the per-feature analysis across workers via parallel::mclapply (falling back to sequential on Windows / one core).

Usage

analyze_methods_parallel(
  iteration_data,
  test_features = NULL,
  alpha = 0.05,
  debug = FALSE,
  n_cores = 1
)

Arguments

Value

Data frame: feature, method, beta, se, pvalue, significant.

Examples

dat <- run_single_iteration(NULL, n_synthetic_samples = 100,
  n_features = 5, n_confounders = 1, seed = 1)
analyze_methods_parallel(dat, n_cores = 1)

Description

Reference-named entry point: applies UNADJ, DIRECT, COCA, IV2SLS, and PGC (matrix bridge) to each tested outcome feature and returns tidy per-feature results with a significance flag. Negative controls are summarised over the full W panel for COCA, so restricting test_features does not change its control summary. The matrix-bridge PGC uses the full W panel (or a validity-screened subset if provided).

Usage

analyze_methods_robust(
  iteration_data,
  test_features = NULL,
  alpha = 0.05,
  debug = FALSE,
  n_cores = 1
)

Arguments

Details

A scalar-bridge variant (fit_pgc_scalar()) is exported for standalone use but is not included in the default pipeline.

Value

Data frame: feature, method, beta, se, pvalue, significant.

Examples

dat <- run_single_iteration(NULL, n_synthetic_samples = 100,
  n_features = 5, n_confounders = 1, seed = 1)
analyze_methods_robust(dat)

Description

Predicts genetically regulated expression by applying published FUSION per-gene weight files (.wgt.RDat) to a genotype dosage matrix. For each gene, the model with the best cross-validation performance (cv.performance["rsq", ]) is used. The resulting predicted-expression matrix is germline-determined and — after PCA via build_w_pcs() — serves as a negative-control panel that captures confounding structure without direct exposure mediation.

Usage

apply_fusion_weights(
  genotypes,
  pos,
  weights_dir = NULL,
  weights = NULL,
  min_snps = 2,
  n_cores = 1
)

Arguments

Value

A list with:

predicted

numeric matrix, genes x samples: predicted expression.

info

data.frame per predicted gene: gene, n_snps, best_method, rsq (cross-validation R2 of the selected model).

Examples

# Simulate two FUSION-style weight objects in memory
dos <- matrix(rbinom(10 * 50, 2, 0.3), nrow = 10,
              dimnames = list(paste0("rs", 1:10), paste0("S", 1:50)))
mk_wgt <- function(snps, w) {
  list(wgt.matrix = matrix(w, ncol = 1, dimnames = list(NULL, "enet")),
       snps = data.frame(V2 = snps),
       cv.performance = matrix(0.2, nrow = 1,
                               dimnames = list("rsq", "enet")))
}
wlist <- list(GeneA = mk_wgt(paste0("rs", 1:5), rep(0.1, 5)),
              GeneB = mk_wgt(paste0("rs", 6:10), rep(-0.2, 5)))
pos <- data.frame(ID = c("GeneA", "GeneB"), WGT = c("a.wgt.RDat", "b.wgt.RDat"))
out <- apply_fusion_weights(dos, pos = pos, weights = wlist)
dim(out$predicted)

Description

S3 generic bridging other data interfaces to the estimation interface. Methods:

Usage

as_iconic_data(input, ...)

## Default S3 method:
as_iconic_data(input, ...)

## S3 method for class 'SummarizedExperiment'
as_iconic_data(
  input,
  assay = 1,
  mediator_assay = NULL,
  exposure,
  instrument = NULL,
  mediator_instrument = NULL,
  negative_controls = NULL,
  covariates = NULL,
  surv_time = NULL,
  surv_event = NULL,
  bin_outcome = NULL,
  outcome_type = c("continuous", "survival", "binary"),
  ...
)

Arguments

Details

• default: a list returned by load_real_input_data(), an existing iconic_data object (returned as-is), or the exposure vector X with named arguments matching iconic_data() (Y, M, G, Gm, W, W1, W2, covariates, feature_names, mediator_names), which delegates to iconic_data().

• SummarizedExperiment: extracts the outcome panel from an assay and sample-level fields (exposure, instruments, negative controls, covariates) from colData. Requires the SummarizedExperiment package (listed under Suggests).

Value

An iconic_data S3 object.

Examples

# From a load_real_input_data() result
input <- load_real_input_data(example = TRUE)
data <- as_iconic_data(input)
print(data)

# From named components (delegates to iconic_data())
data <- as_iconic_data(rnorm(100), Y = matrix(rnorm(100*10), 10, 100),
G = rnorm(100), Gm = rnorm(100),
W = matrix(rnorm(100*10), 10, 100))
if (requireNamespace("SummarizedExperiment", quietly = TRUE)) {
  se <- SummarizedExperiment::SummarizedExperiment(
    assays = list(expr = matrix(rnorm(20 * 60), 20, 60,
                                dimnames = list(paste0("gene", 1:20),
                                                paste0("S", 1:60)))),
    colData = S4Vectors::DataFrame(
      bmi = rnorm(60), prs = rnorm(60),
      nc1 = rnorm(60), nc2 = rnorm(60), age = rnorm(60))
  )
  data <- as_iconic_data(se, assay = "expr", exposure = "bmi",
                         instrument = "prs",
                         negative_controls = c("nc1", "nc2"),
                         covariates = "age")
  print(data)
}

Description

Applies the logit transform $\log_2(\beta / (1 - \beta))$ after clipping betas away from 0 and 1. M-values are approximately homoscedastic and are the recommended scale for linear modeling and PCA of methylation data.

Usage

beta_to_m(beta, clip = 1e-04, drop_nonfinite = TRUE)

Arguments

Value

A numeric matrix of M-values, probes x samples.

Examples

b <- matrix(runif(200, 0.01, 0.99), nrow = 20,
            dimnames = list(paste0("cg", 1:20), paste0("S", 1:10)))
m <- beta_to_m(b)
range(m)

Description

Trains one cross-validated elastic-net model per gene predicting expression from cis SNPs, after residualizing both expression and dosages on covariates (Frisch-Waugh-Lovell). Genes are kept when the model is non-degenerate (at least one non-zero SNP weight), the signed out-of-fold R2 exceeds cv_r2_min, and the one-sided out-of-fold correlation p-value is below cv_p_max. The in-sample predicted expression of passing genes is the composite mediator instrument Gm consumed by iconic_data().

Usage

build_mediator_instruments(
  expression,
  genotypes,
  gene_pos,
  snp_pos,
  covariates = NULL,
  cis_snps = NULL,
  alpha = 0.3,
  nfolds = 5,
  cis_dist = 1e+06,
  cv_r2_min = 0.005,
  cv_p_max = 0.1,
  n_cores = 1,
  seed = NULL,
  verbose = TRUE
)

Arguments

Details

Requires the glmnet package (listed under Suggests).

Value

A list with:

Gm

numeric matrix, genes x samples: in-sample predicted expression for genes passing the quality filter. Pass to iconic_data() as Gm (it is transposed internally).

qc

data.frame of per-gene quality metrics: gene, n_cis_snps, n_nonzero, cv_r2, cor_oof, cv_p, lambda_min, pass.

weights

named list (all trained genes) of data.frames with snp and weight (non-zero coefficients at lambda.min).

Examples

if (requireNamespace("glmnet", quietly = TRUE)) {
  set.seed(42)
  n <- 150; p <- 80
  dos <- matrix(rbinom(p * n, 2, 0.3), nrow = p,
                dimnames = list(paste0("1:", 1:p, ":A:G"),
                                paste0("S", 1:n)))
  expr <- matrix(rnorm(4 * n), nrow = 4,
                 dimnames = list(paste0("Gene", 1:4), paste0("S", 1:n)))
  expr[1, ] <- expr[1, ] + 0.9 * scale(dos[5, ]) + 0.9 * scale(dos[15, ])
  gp <- data.frame(gene = paste0("Gene", 1:4), chr = "1",
                   tss = c(5, 25, 45, 65))
  sp <- data.frame(snp = rownames(dos), chr = "1", pos = 1:p)
  fit <- build_mediator_instruments(expr, dos, gp, sp, seed = 1)
  fit$qc
}

Description

Runs the LDpred2-auto Bayesian shrinkage pipeline (Privé et al.) via the bigsnpr package: match summary statistics to an LD reference, estimate SNP heritability by LD score regression, run a grid of auto chains, filter chains by convergence, and average the per-chain effects. Optionally scores a target genotype set.

Usage

build_prs_ldpred2(
  sumstats,
  ld_ref,
  genotypes = NULL,
  hapmap3 = NULL,
  max_variants = NULL,
  n_chains = 30,
  p_range = c(1e-04, 0.2),
  burn_in = 500,
  num_iter = 500,
  shrink_corr = 0.95,
  allow_jump_sign = FALSE,
  use_mle = FALSE,
  chain_filter_quantile = 0.95,
  impute_mean = TRUE,
  ncores = 1,
  verbose = TRUE
)

Arguments

Details

This function requires the bigsnpr package (listed under Suggests) and a PLINK/bed-format LD reference panel converted with bigsnpr::snp_readBed().

Value

A list with:

beta

data.frame of LDpred2 effect sizes (chr, pos, a0, a1, beta) in LD-reference orientation.

score

numeric vector of per-sample polygenic scores for genotypes (scaled to unit variance), or NULL when genotypes was not supplied.

h2

LD-score-regression SNP heritability estimate used to initialize the chains.

chains_kept

integer indices of chains passing the convergence filter.

qc

list with variant counts at each stage.

Examples

if (requireNamespace("bigsnpr", quietly = TRUE)) {
  # Tiny fake genotype panel standing in for an LD reference
  # (snp_fake genotypes are all-missing mocks; fill with random dosages)
  fake <- bigsnpr::snp_fake(100, 500)
  fake$genotypes[] <- rbinom(100 * 500, 2, 0.3)
  fake$map$chromosome <- 1L
  fake$map$physical.pos <- sort(sample(1:1e6, 500))
  ss <- data.frame(
    chr = fake$map$chromosome, pos = fake$map$physical.pos,
    a0 = fake$map$allele2, a1 = fake$map$allele1,
    beta = rnorm(500, 0, 0.05), beta_se = 0.02, n_eff = 20000
  )
  prs <- build_prs_ldpred2(ss, ld_ref = fake, n_chains = 3,
                           burn_in = 50, num_iter = 50, verbose = FALSE)
  str(prs$beta)
}

Description

Computes the top n_pcs principal components of a features x samples matrix (typically residualized with residualize_matrix()) and returns the sample scores as a W panel (PCs x samples) in the orientation expected by iconic_data(). Uses irlba for fast truncated PCA when available, falling back to stats::prcomp().

Usage

build_w_pcs(x, n_pcs = 20, scale_features = FALSE, prefix = "PC")

Arguments

Value

A list with:

W

numeric matrix, PCs x samples, suitable as the W argument of iconic_data().

variance_explained

numeric vector: percent of total variance explained by each returned PC.

Examples

x <- matrix(rnorm(500 * 60), nrow = 500,
            dimnames = list(paste0("f", 1:500), paste0("S", 1:60)))
w <- build_w_pcs(x, n_pcs = 5)
dim(w$W)
w$variance_explained

Description

For each gene, tests every cis SNP (within cis_dist of the gene start, optionally restricted to a supplied candidate set) for association with expression, after residualizing both expression and dosages on covariates (Frisch-Waugh-Lovell). A lightweight, pure-R alternative to MatrixEQTL for preparing mediator-instrument candidate sets.

Usage

call_cis_eqtls(
  expression,
  genotypes,
  gene_pos,
  snp_pos,
  covariates = NULL,
  cis_snps = NULL,
  cis_dist = 1e+06,
  min_maf = 0.01,
  fdr = 0.05,
  n_cores = 1
)

Arguments

Value

A list with:

all

data.frame of all tested gene-SNP pairs: gene, snp, beta, se, t, p.

best

data.frame with the top SNP per gene (gene, snp, beta, p, q (BH within gene), pass).

Examples

n <- 80
dos <- matrix(rbinom(200 * n, 2, 0.3), nrow = 200,
              dimnames = list(paste0("1:", 1:200, ":A:G"), NULL))
expr <- matrix(rnorm(20 * n), nrow = 20,
               dimnames = list(paste0("Gene", 1:20), NULL))
expr[1, ] <- expr[1, ] + 0.5 * dos[5, ]   # one true cis-eQTL
gp <- data.frame(gene = paste0("Gene", 1:20), chr = "1",
                 tss = seq(1, 191, by = 10))
sp <- data.frame(snp = rownames(dos), chr = "1", pos = 1:200)
hits <- call_cis_eqtls(expr, dos, gp, sp)
head(hits$best)

Description

Regresses the exposure X on the instrument G (plus optional covariates) and reports the partial F statistic and partial R2 for the instrument — the weak-instrument diagnostic. The partial F (not the overall model F) is the relevant quantity when covariates such as ancestry PCs explain much of the exposure.

Usage

check_instrument_strength(G, X, covariates = NULL, min_f = 10)

Arguments

Value

A list with F (partial F statistic), df1, df2, pvalue, partial_r2, n, and weak (TRUE when F < min_f).

Examples

n <- 300
G <- rnorm(n)
X <- 0.3 * G + rnorm(n)
pcs <- matrix(rnorm(n * 3), n, 3, dimnames = list(NULL, paste0("PC", 1:3)))
check_instrument_strength(G, X, covariates = pcs)

Description

Verifies the package is installed, libtorch is present, and a basic tensor op succeeds. Used to gate the GAN path; a FALSE result triggers the multivariate-normal fallback.

Usage

check_torch_setup()

Value

TRUE if torch is usable, otherwise FALSE.

Examples

check_torch_setup()

Description

Computes the Huang (2019) JT-comp composite p-value for testing H0: alpha * beta = 0 against H1: alpha * beta != 0.

Usage

composite_p_value(a, b, var_a = 1, var_b = 1)

Arguments

Details

The composite null decomposes into three cases (alpha=beta=0, alpha!=0/beta=0, alpha=0/beta!=0). The Wald/Sobel test is conservative under the first case because the product of two independent normals follows a normal product distribution, not a normal. This function computes a closed-form p-value that accounts for all three cases without estimating their proportions.

Value

Numeric scalar in $[0, 1]$.

References

Huang, Y.-T. (2019). Genome-wide analyses of sparse mediation effects under composite null hypotheses. Annals of Applied Statistics, 13(1), 60-84.

Examples

composite_p_value(0, 0) # null: p = 1
composite_p_value(2, 2) # strong signal: small p
composite_p_value(2, 0) # H0(2): p = 1 (no mediation)

Description

Implements the Correlated Outcome Control Approach (COCA). Fits w ~ y + X + covars and recovers the causal effect as $\hat\beta = -\hat\beta_X / \hat\beta_Y$. Standard errors are obtained via the delta method.

Usage

fit_coca(y, X, w, covars = NULL, ratio_cap = 10, se_cap = 5)

Arguments

Details

The negative-control W should be an outcome that shares the same unmeasured confounders as Y but has no direct causal path from X.

Value

Named list: beta, se, pvalue. Returns list(beta=NA, se=NA, pvalue=NA) if estimation is unstable (near-zero $\hat\beta_Y$ or extreme ratio).

Examples

set.seed(1)
dat <- generate_toy_data(n = 200, seed = 1)
fit_coca(dat$Y[, 1], dat$X, rowMeans(dat$W))

Description

COCA (Correlated Outcome Control Approach) regresses the negative control W on the outcome Y (W ~ y + X) and recovers the causal effect as a ratio $-\hat\beta_X / \hat\beta_Y$. That identification argument assumes a linear structural outcome model: the W-Y regression coefficient must be proportional to the confounder-Outcome association on the same scale as the causal effect. With a binary outcome the structural model is nonlinear (logistic), so the linear COCA ratio recovers neither the causal log-odds ratio nor the risk difference. COCA is therefore unsupported for binary outcomes and always returns list(beta=NA, se=NA, pvalue=NA).

Usage

fit_coca_bin(y, X, w, covars = NULL, ...)

Arguments

Value

list(beta = NA, se = NA, pvalue = NA) with an informative "reason" attribute.

Examples

set.seed(1)
dat <- generate_toy_data(n = 200, outcome_type = "binary", seed = 1)
fit_coca_bin(dat$y_bin, dat$X, dat$W[, 1])
# $beta [1] NA (COCA unsupported for binary outcomes)

Description

Uses the negative-control outcome W to calibrate both the mediator and outcome regressions, extending fit_coca to the mediation setting.

Usage

fit_coca_mediation(y, X, M, w, covars = NULL, ratio_cap = 10)

Arguments

Details

Stage 1: W ~ M + X -> calibrated alpha_M = -beta_X / beta_M. Stage 2: W ~ Y + X + M -> calibrated NDE = -beta_X / beta_Y, calibrated beta_M = -beta_M / beta_Y. NIE = alpha_M * beta_M (both calibrated).

Value

Named list: NDE, NDE_se, NDE_p, NIE, NIE_se, NIE_p.

Examples

set.seed(1)
dat <- generate_toy_data(n = 200, mo_confounding = 0.8, seed = 1)
fit_coca_mediation(dat$Y[, 1], dat$X, dat$M, rowMeans(dat$W))

Description

COCA mediation regresses the negative control W on the outcome Y in both stages (W ~ M + X, W ~ y + X + M) and recovers the effect as a ratio of regression coefficients. That identification argument assumes a linear structural outcome model; with a binary (nonlinear, logistic) outcome the ratio recovers neither the causal log-odds ratio nor the risk difference. COCA mediation is therefore unsupported for binary outcomes and always returns NA.

Usage

fit_coca_mediation_bin(y, X, M, w, covars = NULL, ...)

Arguments

Value

All-NA list with an informative "reason" attribute.

Examples

set.seed(1)
dat <- generate_toy_data(n = 200, outcome_type = "binary", seed = 1)
fit_coca_mediation_bin(dat$y_bin, dat$X, dat$M, dat$W[, 1])
# all NA (COCA mediation unsupported for binary outcomes)

Description

COCA mediation regresses the negative control W on the outcome Y in both stages (W ~ M + X, W ~ y + X + M), placing the outcome on the right-hand side. This is structurally impossible when the outcome is a survival::Surv object. COCA mediation is therefore unsupported for survival outcomes and always returns NA.

Usage

fit_coca_mediation_surv(time, event, X, M, w, covars = NULL, ...)

Arguments

Value

All-NA list with an informative "reason" attribute.

Examples

set.seed(1)
dat <- generate_toy_data(n = 200, outcome_type = "survival", seed = 1)
fit_coca_mediation_surv(dat$surv_time, dat$surv_event, dat$X, dat$M, dat$W[, 1])
# all NA (COCA mediation unsupported for survival outcomes)

Description

COCA (Correlated Outcome Control Approach) regresses the negative control W on the outcome Y (W ~ y + X) and recovers the causal effect as a ratio $-\hat\beta_X / \hat\beta_Y$. This places the outcome on the right-hand side of the regression, which is structurally impossible when the outcome is a survival::Surv object (time-to-event). COCA is therefore unsupported for survival outcomes and always returns list(beta=NA, se=NA, pvalue=NA).

Usage

fit_coca_surv(time, event, X, w, covars = NULL, ...)

Arguments

Value

list(beta = NA, se = NA, pvalue = NA) with an informative "reason" attribute.

Examples

set.seed(1)
dat <- generate_toy_data(n = 200, outcome_type = "survival", seed = 1)
fit_coca_surv(dat$surv_time, dat$surv_event, dat$X, dat$W[, 1])
# $beta [1] NA (COCA unsupported for survival outcomes)

Description

Regresses Y on X plus the genetic instrument G, the negative-control W, and any additional covariates. This is a "naive" adjustment that uses whatever observables are available but does NOT correct for unmeasured confounding via a ratio or IV approach.

Usage

fit_direct(y, X, g, w, covars = NULL)

Arguments

Value

Named list: beta, se, pvalue.

Examples

set.seed(1)
dat <- generate_toy_data(n = 200, seed = 1)
fit_direct(dat$Y[, 1], dat$X, dat$G[, 1], dat$W[, 1])

Description

Regresses the binary outcome on X plus the genetic instrument G, the negative-control panel W, and covariates. Naive adjustment; does not correct for unmeasured confounding via a ratio or IV approach.

Usage

fit_direct_bin(
  y,
  X,
  g,
  w,
  covars = NULL,
  effect_scale = c("logor", "riskdiff")
)

Arguments

Value

Named list: beta, se, pvalue.

Examples

set.seed(1)
dat <- generate_toy_data(n = 200, outcome_type = "binary", seed = 1)
fit_direct_bin(dat$y_bin, dat$X, dat$G[, 1], dat$W[, 1])

Description

Adjusts for the genetic instrument G and negative-control W in both the mediator and outcome regressions. Like fit_direct, this is a naive adjustment that does not correct for unmeasured confounding via a ratio or IV approach.

Usage

fit_direct_mediation(y, X, M, g, w, covars = NULL)

Arguments

Value

Named list: NDE, NDE_se, NDE_p, NIE, NIE_se, NIE_p.

Examples

set.seed(1)
dat <- generate_toy_data(n = 200, mo_confounding = 0.8, seed = 1)
fit_direct_mediation(dat$Y[, 1], dat$X, dat$M, dat$G[, 1], dat$W[, 1])

Description

Adjusts for the instrument G and negative-control W in both the mediator (OLS) and outcome (logistic / LPM) stages. Naive adjustment.

Usage

fit_direct_mediation_bin(
  y,
  X,
  M,
  g,
  w,
  covars = NULL,
  effect_scale = c("logor", "riskdiff")
)

Arguments

Value

Named list (same fields as fit_unadj_mediation_bin).

Examples

set.seed(1)
dat <- generate_toy_data(n = 200, outcome_type = "binary",
mo_confounding = 0.8, seed = 1)
fit_direct_mediation_bin(dat$y_bin, dat$X, dat$M,
dat$G[, 1], dat$W[, 1])

Description

Adjusts for the instrument G and negative-control W in both the mediator (OLS) and outcome (Cox / RMST) stages. Naive adjustment.

Usage

fit_direct_mediation_surv(
  time,
  event,
  X,
  M,
  g,
  w,
  covars = NULL,
  effect_scale = c("loghr", "rmst"),
  tau = NULL
)

Arguments

Value

Named list (same fields as fit_unadj_mediation_surv).

Examples

set.seed(1)
dat <- generate_toy_data(n = 200, outcome_type = "survival",
mo_confounding = 0.8, seed = 1)
fit_direct_mediation_surv(dat$surv_time, dat$surv_event, dat$X, dat$M,
dat$G[, 1], dat$W[, 1])

Description

Regresses the survival outcome on X plus the genetic instrument G, the negative-control panel W, and covariates. Naive adjustment; does not correct for unmeasured confounding via a ratio or IV approach.

Usage

fit_direct_surv(
  time,
  event,
  X,
  g,
  w,
  covars = NULL,
  effect_scale = c("loghr", "rmst"),
  tau = NULL
)

Arguments

Value

Named list: beta, se, pvalue.

Examples

set.seed(1)
dat <- generate_toy_data(n = 200, outcome_type = "survival", seed = 1)
fit_direct_surv(dat$surv_time, dat$surv_event, dat$X, dat$G[, 1], dat$W[, 1])

Description

Uses the genetic instrument G to instrument for the exposure X, controlling for the negative-control W and any additional covariates. Requires AER.

Usage

fit_iv2sls(y, X, g, w, covars = NULL, min_f = 10)

Arguments

Details

A weak-instrument check is applied using the partial F-statistic for the excluded instrument G (testing G conditional on W and covariates in the first stage), following Stock & Yogo (2005). If the partial F is below min_f, the function returns NA.

Value

Named list: beta, se, pvalue.

Examples

set.seed(1)
dat <- generate_toy_data(n = 200, seed = 1)
fit_iv2sls(dat$Y[, 1], dat$X, dat$G[, 1], dat$W[, 1])

Description

Two-stage predictor substitution (2SPS): the first stage regresses X on the instrument G (plus W and covariates) via OLS, producing fitted $\hat X$; the second stage regresses the binary outcome on $\hat X$ (plus W and covariates) via logistic regression (log-OR) or a linear probability model (risk difference).

Usage

fit_iv2sls_bin(
  y,
  X,
  g,
  w,
  covars = NULL,
  min_f = 10,
  effect_scale = c("logor", "riskdiff")
)

Arguments

Details

A weak-instrument check (partial F for the excluded instrument G, Stock & Yogo 2005) is applied to the OLS first stage. If the partial F is below min_f, the function returns NA.

Value

Named list: beta, se, pvalue.

Examples

set.seed(1)
dat <- generate_toy_data(n = 200, outcome_type = "binary", seed = 1)
fit_iv2sls_bin(dat$y_bin, dat$X, dat$G[, 1], dat$W[, 1])

Description

Uses the genetic instrument G to purge U1 from X, then estimates the mediator and outcome regressions with the cleaned exposure.

Usage

fit_iv2sls_mediation(y, X, M, g, w, covars = NULL, min_f = 10)

Arguments

Details

Strategy:

1. X ~ G + W -> X_hat (purge U1 from X).

2. M ~ X_hat -> alpha_M (clean effect of X on M).

3. Y ~ X_hat + M + W -> NDE = beta_X, beta_M (OLS).

The IV cleans X of U1 confounding, but M remains endogenous via U1 -> M. With a single instrument, natural effects are not fully identified (Rudolph et al., 2024); NDE and NIE are approximations whose bias from M-O confounding is the key finding the simulation demonstrates.

Value

Named list: NDE, NDE_se, NDE_p, NIE, NIE_se, NIE_p.

References

Rudolph, K. E., et al. (2024). Natural direct and indirect effects with an instrumental variable. Biometrics.

Examples

set.seed(1)
dat <- generate_toy_data(n = 300, mo_confounding = 0.8, seed = 1)
fit_iv2sls_mediation(dat$Y[, 1], dat$X, dat$M, dat$G[, 1], dat$W[, 1])

Description

Two-stage predictor substitution with a single instrument (G for X, no mediator instrument):

1. OLS: X ~ g + W + covars -> X_hat (purge U1 from X).

2. OLS: M ~ X_hat + covars -> alpha_M.

3. Logistic / LPM: y ~ X_hat + M + W + covars -> NDE (coef on X_hat), beta_M (coef on M).

NIE = alpha_M * beta_M. Weak-instrument gate (partial F for G) applies to the OLS first stage. Unlike fit_iv2sls_mediation2_bin, this estimator does not instrument the mediator and is NOT point-identified under M-O confounding — it is included for parity with the continuous fit_iv2sls_mediation.

Usage

fit_iv2sls_mediation_bin(
  y,
  X,
  M,
  g,
  w,
  covars = NULL,
  min_f = 10,
  effect_scale = c("logor", "riskdiff")
)

Arguments

Value

Named list (same fields as fit_unadj_mediation_bin). Returns all-NA if the first-stage partial F for G is below min_f.

Examples

set.seed(1)
dat <- generate_toy_data(n = 500, outcome_type = "binary",
mo_confounding = 0.8, seed = 1)
fit_iv2sls_mediation_bin(dat$y_bin, dat$X, dat$M, dat$G[, 1], dat$W[, 1])

Description

Two-stage predictor substitution with a single instrument (G for X, no mediator instrument):

1. OLS: X ~ g + W + covars -> X_hat (purge U1 from X).

2. OLS: M ~ X_hat + covars -> alpha_M.

3. Cox / RMST: Surv(t,e) ~ X_hat + M + W + covars -> NDE (coef on X_hat), beta_M (coef on M).

NIE = alpha_M * beta_M. Weak-instrument gate (partial F for G) applies to the OLS first stage. Unlike fit_iv2sls_mediation2_surv, this estimator does not instrument the mediator and is NOT point-identified under M-O confounding — it is included for parity with the continuous fit_iv2sls_mediation.

Usage

fit_iv2sls_mediation_surv(
  time,
  event,
  X,
  M,
  g,
  w,
  covars = NULL,
  min_f = 10,
  effect_scale = c("loghr", "rmst"),
  tau = NULL
)

Arguments

Value

Named list (same fields as fit_unadj_mediation_surv).

Examples

set.seed(1)
dat <- generate_toy_data(n = 200, outcome_type = "survival", seed = 1)
fit_iv2sls_mediation_surv(dat$surv_time, dat$surv_event, dat$X, dat$M,
dat$G[, 1], dat$W[, 1])

Description

Uses TWO genetic instruments: G for the exposure X and Gm for the mediator M. This is the key extension: by instrumenting both endogenous variables, NDE and NIE become point-identified even under mediator-outcome (M-O) confounding, resolving the identification failure that limits the single-instrument estimators.

Usage

fit_iv2sls_mediation2(
  y,
  X,
  M,
  g,
  gm,
  covars = NULL,
  min_f = 10,
  W1 = NULL,
  W2 = NULL,
  w = NULL
)

Arguments

Details

The motivating design uses two distinct genetic instruments: a mediator instrument (e.g. a cis-eQTL) instruments the mediator (M), while a separate exposure instrument (e.g. a polygenic score) instruments the exposure (X). The mediator set must be restricted to features for which a mediator instrument is available.

Strategy (sequential 2SLS, three OLS stages), with optional path-specific negative-control (NC) augmentation:

1. X ~ G (+ W1) + covars -> X_hat (purge U1 from X). Weak-IV check: partial F for G >= min_f.

2. M ~ X_hat + Gm (+ W2) + covars -> M_hat, alpha_M = coef on X_hat. Weak-IV check: partial F for Gm >= min_f.

3. Y ~ X_hat + M_hat (+ W2) + covars -> NDE = beta_X_hat, beta_M = coef on M_hat.

NIE = alpha_M * beta_M (delta-method SE).

Path-specific NC augmentation (optional). W1 proxies the exposure-mediator confounder (U1, the X->M path) and is added to stage 1 only; W2 proxies the mediator-outcome confounder (U2, the M->Y path) and is added to stages 2 and 3. Either panel may be omitted: with W1 = NULL stage 1 is unaugmented, with W2 = NULL stages 2–3 are unaugmented, and with both NULL the estimator reduces to plain two-instrument 2-stage MR. Conditioning on a pooled panel (W1 + W2) / 2 in all three stages is a collider under multi-confounder designs (the pooled panel is a common child of the independent confounders U1 and U2) and is therefore not supported: if W1 and W2 are identical they are treated as absent (pure MR). Pass the two path-specific panels separately to obtain the coverage-improves-accuracy behaviour.

Path-specific augmentation is only defined when W1 and W2 proxy distinct confounders. When the two panels are distinct-noise proxies of the same latent composite (the single-confounder / shared-loading design), their column spaces are near-collinear and stage-1 W1 would inject the shared M->Y confounder into X_hat. The estimator detects this (leading-PC correlation between the panels) and falls back to plain two-instrument 2-stage MR, exactly as for an identical pooled panel. Genuinely distinct panels are unaffected.

Including X_hat in the M first-stage (stage 2) is essential: M structurally depends on X, so the X -> M path must be captured for alpha_M and the NIE to be correctly estimated. Gm provides the exogenous variation that identifies the M -> Y effect net of U1 confounding.

The sequential 2SLS implementation matches the pattern of fit_iv2sls_mediation (three separate OLS stages). Standard errors in the outcome stage do not account for the two-stage estimation of M_hat, consistent with the existing estimators; the simulation benchmarks bias, not SE accuracy. A unit test cross-validates this estimator's NDE and beta_M against AER::ivreg on the just-identified system (Y ~ X + M | G + Gm, pure-MR form).

Value

Named list: NDE, NDE_se, NDE_p, NIE, NIE_se, NIE_p. Returns all-NA if either first-stage partial F is below min_f.

References

Rudolph, K. E., et al. (2024). Natural direct and indirect effects with an instrumental variable. Biometrics.

Examples

set.seed(1)
dat <- generate_toy_data(n = 500, mo_confounding = 0.8,
phi = 0.8, lambda_XM = c(1, 0), lambda_MY = c(0, 1),
omega_1 = 0.7, omega_2 = 0.7, seed = 1)
fit_iv2sls_mediation2(dat$Y[, 1], dat$X, dat$M, dat$G[, 1], dat$Gm,
W1 = dat$W1, W2 = dat$W2)

Description

Two-stage Mendelian-randomization mediation with a binary outcome stage (2SPS), with optional path-specific negative-control (NC) augmentation:

1. OLS: X ~ g (+ W1) + covars -> X_hat (purge U1 from X).

2. OLS: M ~ X_hat + gm (+ W2) + covars -> M_hat, alpha_M.

3. Logistic / LPM: y ~ X_hat + M_hat (+ W2) + covars -> NDE (coef on X_hat), beta_M (coef on M_hat).

NIE = alpha_M * beta_M. Weak-instrument gates (partial F for G and Gm) apply to the OLS first stages.

Usage

fit_iv2sls_mediation2_bin(
  y,
  X,
  M,
  g,
  gm,
  covars = NULL,
  min_f = 10,
  W1 = NULL,
  W2 = NULL,
  w = NULL,
  effect_scale = c("logor", "riskdiff")
)

Arguments

Details

Path-specific NC augmentation (optional). W1 proxies the exposure-mediator confounder (U1, X->M path) and is added to stage 1 only; W2 proxies the mediator-outcome confounder (U2, M->Y path) and is added to stages 2 and 3. Either panel may be omitted; with both NULL the estimator reduces to plain two-instrument 2-stage MR. Conditioning on a pooled panel in all three stages is a collider under multi-confounder designs and is not supported: identical W1/W2 are treated as absent (pure MR). When W1 and W2 are distinct-noise proxies of the same latent composite (single-confounder design), their column spaces are near-collinear and the estimator likewise falls back to plain two-instrument 2-stage MR.

Value

Named list (same fields as fit_unadj_mediation_bin). Returns all-NA if either first-stage partial F is below min_f.

References

Rudolph, K. E., et al. (2024). Natural direct and indirect effects with an instrumental variable. Biometrics.

Examples

set.seed(1)
dat <- generate_toy_data(n = 500, outcome_type = "binary",
mo_confounding = 0.8, phi = 0.8, lambda_XM = c(1, 0),
lambda_MY = c(0, 1), omega_1 = 0.7, omega_2 = 0.7, seed = 1)
fit_iv2sls_mediation2_bin(dat$y_bin, dat$X, dat$M,
dat$G[, 1], dat$Gm, W1 = dat$W1, W2 = dat$W2)

Description

Two-stage Mendelian-randomization mediation with a survival outcome stage (2SPS), with optional path-specific negative-control (NC) augmentation:

1. OLS: X ~ g (+ W1) + covars -> X_hat (purge U1 from X).

2. OLS: M ~ X_hat + gm (+ W2) + covars -> M_hat, alpha_M.

3. Cox / RMST: Surv(t,e) ~ X_hat + M_hat (+ W2) + covars -> NDE (coef on X_hat), beta_M (coef on M_hat).

NIE = alpha_M * beta_M. Weak-instrument gates (partial F for G and Gm) apply to the OLS first stages.

Usage

fit_iv2sls_mediation2_surv(
  time,
  event,
  X,
  M,
  g,
  gm,
  covars = NULL,
  min_f = 10,
  W1 = NULL,
  W2 = NULL,
  w = NULL,
  effect_scale = c("loghr", "rmst"),
  tau = NULL
)

Arguments

Details

Path-specific NC augmentation (optional). W1 proxies the exposure-mediator confounder (U1, X->M path) and is added to stage 1 only; W2 proxies the mediator-outcome confounder (U2, M->Y path) and is added to stages 2 and 3. Either panel may be omitted; with both NULL the estimator reduces to plain two-instrument 2-stage MR. Conditioning on a pooled panel in all three stages is a collider under multi-confounder designs and is not supported: identical W1/W2 are treated as absent (pure MR). When W1 and W2 are distinct-noise proxies of the same latent composite (single-confounder design), their column spaces are near-collinear and the estimator likewise falls back to plain two-instrument 2-stage MR.

Value

Named list (same fields as fit_unadj_mediation_surv). Returns all-NA if either first-stage partial F is below min_f.

References

Rudolph, K. E., et al. (2024). Natural direct and indirect effects with an instrumental variable. Biometrics.

Examples

set.seed(1)
dat <- generate_toy_data(n = 500, outcome_type = "survival",
mo_confounding = 0.8, phi = 0.8, lambda_XM = c(1, 0),
lambda_MY = c(0, 1), omega_1 = 0.7, omega_2 = 0.7, seed = 1)
fit_iv2sls_mediation2_surv(dat$surv_time, dat$surv_event, dat$X, dat$M,
dat$G[, 1], dat$Gm, W1 = dat$W1, W2 = dat$W2)

Description

Two-stage predictor substitution (2SPS): the first stage regresses X on the instrument G (plus W and covariates) via OLS, producing fitted $\hat X$; the second stage regresses the survival outcome on $\hat X$ (plus W and covariates) via Cox (log-HR) or RMST pseudo-observation OLS.

Usage

fit_iv2sls_surv(
  time,
  event,
  X,
  g,
  w,
  covars = NULL,
  min_f = 10,
  effect_scale = c("loghr", "rmst"),
  tau = NULL
)

Arguments

Details

A weak-instrument check (partial F for the excluded instrument G, Stock & Yogo 2005) is applied to the OLS first stage. If the partial F is below min_f, the function returns NA.

Value

Named list: beta, se, pvalue.

Examples

set.seed(1)
dat <- generate_toy_data(n = 200, outcome_type = "survival", seed = 1)
fit_iv2sls_surv(dat$surv_time, dat$surv_event, dat$X, dat$G[, 1], dat$W[, 1])

Description

A three-step bridge-function estimator:

1. Residualises X on G to isolate the U-driven component X_resid.

2. Regresses X_resid on the FULL W matrix to construct W_hat, a proxy for unmeasured confounding. This step requires ncol(W) >= k (the proximal completeness condition): if W has fewer valid columns than confounders, the bridge cannot span the confounder subspace and the estimator is under-identified.

3. Fits Y ~ X + W_hat to absorb confounding bias.

Usage

fit_pgc(y, X, g, W, covars = NULL)

Arguments

Details

Unlike the scalar version (fit_pgc_scalar), which collapses W to rowMeans(W) and is algebraically equivalent to IV/2SLS, the matrix bridge preserves the dimensional structure of W and is the estimator for which the completeness condition is binding.

Value

Named list: beta, se, pvalue.

Examples

set.seed(1)
dat <- generate_toy_data(n = 200, seed = 1)
fit_pgc(dat$Y[, 1], dat$X, dat$G[, 1], dat$W)

Description

Three-step bridge-function estimator with a binary outcome stage:

1. Residualise X on G -> X_resid (OLS).

2. Bridge X_resid on the FULL W matrix -> W_hat (OLS).

3. Regress the binary outcome on X + W_hat via logistic regression (log-OR) or a linear probability model (risk difference).

Usage

fit_pgc_bin(y, X, g, W, covars = NULL, effect_scale = c("logor", "riskdiff"))

Arguments

Value

Named list: beta, se, pvalue.

Examples

set.seed(1)
dat <- generate_toy_data(n = 200, outcome_type = "binary", seed = 1)
fit_pgc_bin(dat$y_bin, dat$X, dat$G[, 1], dat$W)

Description

Extends fit_pgc (matrix bridge) to the mediation setting by constructing a confounding proxy $\hat W$ from the full W matrix and including it in both the mediator and outcome regressions.

Usage

fit_pgc_mediation(y, X, M, g, W, covars = NULL)

Arguments

Details

Steps:

1. Residualise X on G -> X_resid (U-driven residual).

2. Bridge X_resid on the FULL W matrix -> W_hat (proxy for U).

3. M ~ X + W_hat -> alpha_M (adjusted for confounding proxy).

4. Y ~ X + M + W_hat -> NDE = beta_X, beta_M (adjusted).

The matrix bridge requires ncol(W) >= k (proximal completeness) for the bridge to span the confounder subspace.

Value

Named list: NDE, NDE_se, NDE_p, NIE, NIE_se, NIE_p.

Examples

set.seed(1)
dat <- generate_toy_data(n = 200, mo_confounding = 0.8, seed = 1)
fit_pgc_mediation(dat$Y[, 1], dat$X, dat$M, dat$G[, 1], dat$W)

Description

Bridge-function-adjusted natural direct and indirect effects with a binary outcome stage, using a single negative-control panel W for both the X->M and M->Y confounding paths. This is the binary analogue of fit_pgc_mediation (continuous outcome):

1. OLS: residualise X on G -> X_resid.

2. OLS: bridge X_resid on the FULL W matrix -> W_hat (proxy for U).

3. OLS: M ~ X + W_hat + covars -> alpha_M.

4. Logistic / LPM: y ~ X + M + W_hat + covars -> NDE (coef on X), beta_M (coef on M).

NIE = alpha_M * beta_M. Unlike fit_pgc_mediation2_bin, which uses path-specific W1/W2 bridges, this estimator uses a single combined W panel and is appropriate when separate conf_XM / conf_MY confounders are not assumed.

Usage

fit_pgc_mediation_bin(
  y,
  X,
  M,
  g,
  W,
  covars = NULL,
  min_f = 10,
  effect_scale = c("logor", "riskdiff")
)

Arguments

Value

Named list (same fields as fit_unadj_mediation_bin).

Examples

set.seed(1)
dat <- generate_toy_data(n = 500, outcome_type = "binary",
mo_confounding = 0.8, seed = 1)
fit_pgc_mediation_bin(dat$y_bin, dat$X, dat$M, dat$G[, 1], dat$W)

Description

Bridge-function-adjusted natural direct and indirect effects with a survival outcome stage, using a single negative-control panel W for both the X->M and M->Y confounding paths. This is the survival analogue of fit_pgc_mediation (continuous outcome):

1. OLS: residualise X on G -> X_resid.

2. OLS: bridge X_resid on the FULL W matrix -> W_hat (proxy for U).

3. OLS: M ~ X + W_hat + covars -> alpha_M.

4. Cox / RMST: Surv(t,e) ~ X + M + W_hat + covars -> NDE (coef on X), beta_M (coef on M).

NIE = alpha_M * beta_M. Unlike fit_pgc_mediation2_surv, which uses path-specific W1/W2 bridges, this estimator uses a single combined W panel and is appropriate when separate conf_XM / conf_MY confounders are not assumed.

Usage

fit_pgc_mediation_surv(
  time,
  event,
  X,
  M,
  g,
  W,
  covars = NULL,
  min_f = 10,
  effect_scale = c("loghr", "rmst"),
  tau = NULL
)

Arguments

Value

Named list (same fields as fit_unadj_mediation_surv). Returns all-NA if the first-stage partial F for G is below min_f.

Examples

set.seed(1)
dat <- generate_toy_data(n = 200, outcome_type = "survival", seed = 1)
fit_pgc_mediation_surv(dat$surv_time, dat$surv_event, dat$X, dat$M,
dat$G[, 1], dat$W)

Description

Extends the proximal-inference approach to the two-linked-DAG mediation setting. Uses path-specific negative controls — W1 capturing the X->M confounder conf_XM and W2 capturing the M->Y confounder conf_MY — to purge confounding at both stages via bridge functions.

Usage

fit_pgc_mediation2(y, X, M, g, W1, W2, gm = NULL, covars = NULL, min_f = 10)

Arguments

Details

When gm is NULL (no mediator instrument), Stage 2 residualises M on the cleaned exposure X_hat to isolate the conf_MY-driven component, then bridges on W2. This is pure negative-control identification at both stages — no mediator instrument is required.

When gm is supplied (mediator instrument present but its exogeneity may be in doubt), Stage 2 uses Gm to help isolate conf_MY's effect on M before bridging on W2. The bridge W_hat_M does the confounding removal, so the estimator is robust to Gm being correlated with conf_MY — the residual bias from Gm-U correlation is smaller than IV2SLS2's.

Strategy (three stages):

1. Bridge for X (purge conf_XM from X): X_resid = residuals(X ~ G1 + C); W_hat_X = bridge(X_resid ~ W1) (fitted values, proxy for conf_XM); X_hat = fitted(X ~ G1 + W_hat_X + C). Weak-IV check: partial F for G1 >= min_f.

2. Bridge for M (purge conf_MY from M): M_resid = residuals(M ~ X_hat + C) (gm = NULL), or M_resid = residuals(M ~ Gm + C) (gm supplied); W_hat_M = bridge(M_resid ~ W2) (fitted values, proxy for conf_MY); M_hat = fitted(M ~ X_hat + W_hat_M + C) (gm = NULL), or M_hat = fitted(M ~ X_hat + Gm + W_hat_M + C) (gm supplied); alpha_M = coefficient on X_hat.

3. Outcome: Y ~ X_hat + M_hat + W_hat_X + W_hat_M + C; NDE = coefficient on X_hat, beta_M = coefficient on M_hat.

NIE = alpha_M * beta_M (delta-method SE).

The key advantage over fit_iv2sls_mediation2: PGC-2 does not require instrument exogeneity. When the mediator instrument Gm is correlated with the confounder conf_MY (rho_G2 > 0), IV2SLS2 is biased but PGC-2's bridge absorbs conf_MY regardless of the instrument violation. The tipping-point simulation maps where PGC-2 bias crosses below IV2SLS2 bias.

Value

Named list: NDE, NDE_se, NDE_p, NIE, NIE_se, NIE_p. Returns all-NA if the first-stage partial F for G1 is below min_f.

References

Miao, W., Geng, Z., & Tchetgen Tchetgen, E. (2018). Identifying causal effects with proxy variables of an unmeasured confounder. Biometrika, 105(4), 987-993.

Examples

set.seed(1)
dat <- generate_toy_data(n = 500, mo_confounding = 0.8,
rho_G2 = 0.3, lambda_XM = c(1, 0), lambda_MY = c(0, 1),
omega_1 = 0.7, omega_2 = 0.7, seed = 1)
# Without mediator instrument (pure NC identification)
fit_pgc_mediation2(dat$Y[, 1], dat$X, dat$M, dat$G1, dat$W1, dat$W2)
# With (possibly imperfect) mediator instrument
fit_pgc_mediation2(dat$Y[, 1], dat$X, dat$M, dat$G1, dat$W1, dat$W2,
gm = dat$Gm)

Description

Two-stage proximal mediation with path-specific negative controls and a binary outcome stage. Stages 1-2 (X bridge on W1, M bridge on W2, X_hat / M_hat construction) are identical to fit_pgc_mediation2 and remain OLS. Only stage 3 switches to logistic regression (log-OR) or a linear probability model (risk difference): y ~ X_hat + M_hat + W_hat_X + W_hat_M + covars.

Usage

fit_pgc_mediation2_bin(
  y,
  X,
  M,
  g,
  W1,
  W2,
  gm = NULL,
  covars = NULL,
  min_f = 10,
  effect_scale = c("logor", "riskdiff")
)

Arguments

Details

When gm = NULL (PGC2), stage 2 uses pure NC identification. When gm is supplied (PGC2Gm), the mediator instrument helps isolate conf_MY before bridging W2.

Value

Named list (same fields as fit_unadj_mediation_bin). Returns all-NA if the first-stage partial F for G1 is below min_f.

Examples

set.seed(1)
dat <- generate_toy_data(n = 500, outcome_type = "binary",
mo_confounding = 0.8, rho_G2 = 0.3,
lambda_XM = c(1, 0), lambda_MY = c(0, 1), seed = 1)
fit_pgc_mediation2_bin(dat$y_bin, dat$X, dat$M,
dat$G[, 1], dat$W1, dat$W2, gm = dat$Gm)

Description

Two-stage proximal mediation with path-specific negative controls and a survival outcome stage. Stages 1-2 (X bridge on W1, M bridge on W2, X_hat / M_hat construction) are identical to fit_pgc_mediation2 and remain OLS. Only stage 3 switches to Cox (log-HR) or RMST pseudo-observation OLS: Surv(t,e) ~ X_hat + M_hat + W_hat_X + W_hat_M + covars.

Usage

fit_pgc_mediation2_surv(
  time,
  event,
  X,
  M,
  g,
  W1,
  W2,
  gm = NULL,
  covars = NULL,
  min_f = 10,
  effect_scale = c("loghr", "rmst"),
  tau = NULL
)

Arguments

Details

When gm = NULL (PGC2), stage 2 uses pure NC identification. When gm is supplied (PGC2Gm), the mediator instrument helps isolate conf_MY before bridging W2.

Value

Named list (same fields as fit_unadj_mediation_surv). Returns all-NA if the first-stage partial F for G1 is below min_f.

Examples

set.seed(1)
dat <- generate_toy_data(n = 500, outcome_type = "survival",
mo_confounding = 0.8, rho_G2 = 0.3,
lambda_XM = c(1, 0), lambda_MY = c(0, 1), seed = 1)
fit_pgc_mediation2_surv(dat$surv_time, dat$surv_event, dat$X, dat$M,
dat$G[, 1], dat$W1, dat$W2, gm = dat$Gm)

Description

The original ICONIC PGC implementation, which summarises the negative-control panel as a scalar (rowMeans(W)) before bridging. This version is numerically stable and works in small samples, but the scalar bridge produces $\hat W$ proportional to the G-residualised exposure by construction, making the estimator algebraically equivalent to IV/2SLS when the instrument is valid. The proximal completeness condition (dim(W_valid) >= k) is NOT binding for this estimator.

Usage

fit_pgc_scalar(y, X, g, w, covars = NULL)

Arguments

Details

Use fit_pgc (matrix bridge) when the completeness condition is of interest. Use this function as a stable fallback when n is small relative to the number of valid controls, or as an IV-equivalent benchmark.

Value

Named list: beta, se, pvalue.

Examples

set.seed(1)
dat <- generate_toy_data(n = 200, seed = 1)
fit_pgc_scalar(dat$Y[, 1], dat$X, dat$G[, 1], rowMeans(dat$W))

Description

The original ICONIC PGC mediation implementation, which summarises the negative-control panel as a scalar (rowMeans(W)) before bridging. Like fit_pgc_scalar, this version is algebraically equivalent to the IV2SLS mediation estimator when the instrument is valid. Use fit_pgc_mediation (matrix bridge) when the completeness condition is of interest.

Usage

fit_pgc_scalar_mediation(y, X, M, g, w, covars = NULL)

Arguments

Value

Named list: NDE, NDE_se, NDE_p, NIE, NIE_se, NIE_p.

Examples

set.seed(1)
dat <- generate_toy_data(n = 200, mo_confounding = 0.8, seed = 1)
fit_pgc_scalar_mediation(dat$Y[, 1], dat$X, dat$M, dat$G[, 1], rowMeans(dat$W))

Description

Three-step bridge-function estimator with a survival outcome stage:

1. Residualise X on G -> X_resid (OLS).

2. Bridge X_resid on the FULL W matrix -> W_hat (OLS).

3. Regress the survival outcome on X + W_hat via Cox (log-HR) or RMST pseudo-observation OLS.

Usage

fit_pgc_surv(
  time,
  event,
  X,
  g,
  W,
  covars = NULL,
  effect_scale = c("loghr", "rmst"),
  tau = NULL
)

Arguments

Value

Named list: beta, se, pvalue.

Examples

set.seed(1)
dat <- generate_toy_data(n = 200, outcome_type = "survival", seed = 1)
fit_pgc_surv(dat$surv_time, dat$surv_event, dat$X, dat$G[, 1], dat$W)

Description

Regresses the binary outcome on the exposure X (plus covariates) with no instrument or negative-control adjustment. Bias reference.

Usage

fit_unadj_bin(y, X, covars = NULL, effect_scale = c("logor", "riskdiff"))

Arguments

Details

When effect_scale = "logor" (default), fits a logistic regression and returns the conditional log-odds ratio for X. When effect_scale = "riskdiff", fits a linear probability model (OLS on the 0/1 outcome), returning a risk difference — a collapsible, linear alternative to the non-collapsible odds ratio.

Value

Named list: beta, se, pvalue.

Examples

set.seed(1)
dat <- generate_toy_data(n = 200, outcome_type = "binary", seed = 1)
fit_unadj_bin(dat$y_bin, dat$X)
fit_unadj_bin(dat$y_bin, dat$X, effect_scale = "riskdiff")

Description

Stage 1: M ~ X (estimate alpha_M). Stage 2: Y ~ X + M (estimate NDE = beta_X, beta_M). NIE = alpha_M * beta_M.

Usage

fit_unadj_mediation(y, X, M, covars = NULL)

Arguments

Details

Does not adjust for unmeasured confounding. Provided as a bias reference floor, analogous to UNADJ in the total-effect setting.

Value

Named list: NDE, NDE_se, NDE_p, NIE, NIE_se, NIE_p.

Examples

set.seed(1)
dat <- generate_toy_data(n = 200, mo_confounding = 0.8, seed = 1)
fit_unadj_mediation(dat$Y[, 1], dat$X, dat$M)

Description

Stage 1 (OLS): M ~ X -> alpha_M. Stage 2 (logistic / LPM): binary outcome ~ X + M -> NDE (coef on X), beta_M (coef on M). NIE = alpha_M * beta_M. No confounding adjustment; bias reference.

Usage

fit_unadj_mediation_bin(
  y,
  X,
  M,
  covars = NULL,
  effect_scale = c("logor", "riskdiff")
)

Arguments

Value

Named list: NDE, NDE_se, NDE_p, NIE, NIE_se, NIE_p, alpha_M, alpha_se, beta_M, beta_M_se.

Examples

set.seed(1)
dat <- generate_toy_data(n = 200, outcome_type = "binary",
mo_confounding = 0.8, seed = 1)
fit_unadj_mediation_bin(dat$y_bin, dat$X, dat$M)

Description

Stage 1 (OLS): M ~ X -> alpha_M. Stage 2 (Cox / RMST): survival outcome ~ X + M -> NDE (coef on X), beta_M (coef on M). NIE = alpha_M * beta_M. No confounding adjustment; bias reference.

Usage

fit_unadj_mediation_surv(
  time,
  event,
  X,
  M,
  covars = NULL,
  effect_scale = c("loghr", "rmst"),
  tau = NULL
)

Arguments

Value

Named list: NDE, NDE_se, NDE_p, NIE, NIE_se, NIE_p, alpha_M, alpha_se, beta_M, beta_M_se.

Examples

set.seed(1)
dat <- generate_toy_data(n = 200, outcome_type = "survival",
mo_confounding = 0.8, seed = 1)
fit_unadj_mediation_surv(dat$surv_time, dat$surv_event, dat$X, dat$M)

Description

Regresses the survival outcome on the exposure X (plus covariates) with no instrument or negative-control adjustment. Bias reference.

Usage

fit_unadj_surv(
  time,
  event,
  X,
  covars = NULL,
  effect_scale = c("loghr", "rmst"),
  tau = NULL
)

Arguments

Details

When effect_scale = "loghr" (default), fits a Cox proportional-hazards model and returns the log-hazard ratio for X. When effect_scale = "rmst", regresses leave-one-out RMST pseudo-observations (Graw et al. 2009) on X via OLS, returning an effect on the restricted-mean-survival-time (time) scale — a collapsible, linear alternative to the non-collapsible hazard ratio.

Value

Named list: beta, se, pvalue.

Examples

set.seed(1)
dat <- generate_toy_data(n = 200, outcome_type = "survival", seed = 1)
fit_unadj_surv(dat$surv_time, dat$surv_event, dat$X)
fit_unadj_surv(dat$surv_time, dat$surv_event, dat$X, effect_scale = "rmst")

Description

For each scenario in the grid, generates n_iter synthetic datasets with mediator-outcome confounding (via mo_confounding), runs every mediation estimator, and summarises NDE/NIE bias / RMSE / Type I error. This is the mediation analogue of gan_sensitivity().

Usage

gan_mediation_sensitivity(
  trained_gan = NULL,
  conf_grid = c(0.2, 0.5, 0.8),
  coverage_grid = c(0.3, 0.7, 1),
  k_grid = 1,
  mo_confounding = 0.8,
  phi = 0,
  rho_G1 = 0,
  rho_G2 = 0,
  rho_pop = 0,
  lambda_XM = NULL,
  lambda_MY = NULL,
  omega_1 = NULL,
  omega_2 = NULL,
  nc_model = "proxy",
  n_iter = 50,
  n_samples = 500,
  n_features = 20,
  beta_X = 0.1,
  alpha_M = 0.5,
  beta_M = 0.3,
  base_seed = 750,
  n_cores = 1,
  outcome_type = c("continuous", "survival", "binary"),
  effect_scale = c("loghr", "rmst", "logor", "riskdiff"),
  surv_h0 = 0.1,
  surv_event_frac = 0.6,
  surv_censor_rate = NULL,
  bin_prev = 0.5
)

Arguments

Details

When phi > 0, a mediator-specific genetic instrument (Gm) is generated and the 2-stage MR estimator (IV2SLS2) is included in the results, enabling point identification of NDE/NIE under M-O confounding.

Value

A list with summary (one row per scenario x method, with conf_strength, coverage, k, mo_confounding, phi, true_NDE, true_NIE and NDE/NIE bias/RMSE/Type I columns) and grid.

Examples

sens <- gan_mediation_sensitivity(NULL, conf_grid = 0.8,
  coverage_grid = 0.7, mo_confounding = 0.8,
  n_iter = 2, n_samples = 100, n_features = 5)
head(sens$summary)

Description

For each cell in the grid (pleiotropy strength x confounding strength), generates n_iter synthetic datasets with run_single_iteration(), runs every estimator, and summarises bias / RMSE / power. Two arms are run per cell:

Usage

gan_pleiotropy_sensitivity(
  trained_gan = NULL,
  pleio_grid = c(0, 0.05, 0.1),
  conf_grid = c(0.2, 0.5, 0.8),
  tau = 0.25,
  nc_model = "proxy",
  n_iter = 50,
  n_samples = 500,
  n_features = 10,
  coverage = 0.7,
  k = 1,
  base_seed = 900,
  n_cores = 1
)

Arguments

Details

• alternative (effect_size = tau): the true total effect is tau, so power is the empirical power to detect it.

• null (effect_size = 0): the true total effect is 0, so power is the empirical Type I error rate.

The pleio parameter adds a direct G -> Y path of the requested strength, violating the exclusion restriction. IV/2SLS is consistent only when pleio = 0; any pleio > 0 introduces bias that does not shrink with sample size.

Value

A list with summary (one row per cell x arm x method, with pleio, conf_strength, arm, true_total, and the columns from summarise_results()) and grid.

Examples

sens <- gan_pleiotropy_sensitivity(NULL,
  pleio_grid = c(0, 0.10), conf_grid = 0.8,
  n_iter = 2, n_samples = 100, n_features = 5)
head(sens$summary)

Description

For each scenario in the grid, generates n_iter synthetic datasets with run_single_iteration() (using the trained generator for texture), runs every estimator, and summarises bias / RMSE / power. This is the multi-confounder, negative-control-aware generalisation of the package's parameter sweeps.

Usage

gan_sensitivity(
  trained_gan = NULL,
  conf_grid = c(0.2, 0.5, 0.8),
  coverage_grid = c(0.3, 0.7, 1),
  k_grid = 1,
  nc_model = "proxy",
  n_iter = 50,
  n_samples = 500,
  n_features = 20,
  beta_X = 0.1,
  alpha_M = 0.5,
  beta_M = 0.3,
  effect_size = NULL,
  base_seed = 700,
  n_cores = 1,
  outcome_type = c("continuous", "survival", "binary"),
  effect_scale = c("loghr", "rmst", "logor", "riskdiff"),
  surv_h0 = 0.1,
  surv_event_frac = 0.6,
  surv_censor_rate = NULL,
  bin_prev = 0.5
)

Arguments

Value

A list with summary (one row per scenario x method, with conf_strength, coverage, k, true_total and the columns from summarise_results) and grid (the scenario grid).

Examples

sens <- gan_sensitivity(NULL, conf_grid = 0.8, coverage_grid = 0.7,
  n_iter = 2, n_samples = 100, n_features = 5)
head(sens$summary)

Description

When mo_confounding > 0, the unmeasured confounder U1 also affects the mediator M, creating mediator-outcome confounding — the key extension for the mediation simulation. When mo_confounding = 0 (default), the original DGP is preserved exactly.

Usage

generate_toy_data(
  n = 500,
  n_features = 20,
  n_mediators = 1,
  beta_X = 0.1,
  alpha_M = 0.5,
  beta_M = 0.3,
  conf_str = 0.8,
  w_signal = 0.7,
  mo_confounding = 0,
  pleio = 0,
  phi = 0,
  gamma_G = 0.6,
  rho_G1 = 0,
  rho_G2 = 0,
  rho_pop = 0,
  lambda_XM = NULL,
  lambda_MY = NULL,
  omega_1 = NULL,
  omega_2 = NULL,
  feat_cor = 0,
  separate_U = NULL,
  u_strength = NULL,
  w_coverage_profile = NULL,
  outcome_type = c("continuous", "survival", "binary"),
  surv_h0 = 0.1,
  surv_event_frac = 0.6,
  surv_censor_rate = NULL,
  bin_prev = 0.5,
  seed = NULL,
  beta_Z = NULL
)