mvfmr: Multivariable Functional Mendelian Randomization

@misc{fontana2025unravelingtimevaryingcausaleffects,
      title={Unraveling time-varying causal effects of multiple exposures: integrating Functional Data Analysis with Multivariable Mendelian Randomization}, 
      author={Nicole Fontana and Francesca Ieva and Luisa Zuccolo and Emanuele Di Angelantonio and Piercesare Secchi},
      year={2025},
      eprint={2512.19064},
      archivePrefix={arXiv},
      primaryClass={stat.AP},
      url={https://arxiv.org/abs/2512.19064}, 
}

install.packages("mvfmr")

install.packages("devtools") # Install devtools if not already installed
devtools::install_github("NicoleFontana/mvfmr")

# Run the manuscript simulations
demo("tests_manuscript", package = "mvfmr")

# Run the MV-FMR tutorial
source(system.file("examples", "test_MV-FMR.R", package = "mvfmr"))

# Or just locate it to open/copy/adapt it:
system.file("examples", "test_MV-FMR.R", package = "mvfmr")

# Run the U-FMR tutorial
source(system.file("examples", "test_U-FMR.R", package = "mvfmr"))

library(mvfmr)
library(fdapace)

# Step 1: Simulate exposure data (m = 2 exposures)
set.seed(12345)
sim_data <- getX_multi_exposure(
  N = 1000,           # Sample size
  J = 50,             # Number of genetic instruments
  nSparse = 10,        # Sparse observations per subject
  n_exposures = 2      # Number of exposures (m)
)

# Step 2: Generate outcome
outcome_data <- getY_multi_exposure(
  sim_data,
  XYmodels = c("2", "8"),     # Linear effect for exposure 1, quadratic for exposure 2
  X_effects = c(TRUE, TRUE),
  outcome_type = "continuous"
)

# Step 3: Functional PCA for each exposure
fpca_results <- lapply(sim_data$exposures, function(exp_k) {
  FPCA(
    exp_k$Ly_sim,
    exp_k$Lt_sim,
    list(dataType = 'Sparse', error = TRUE, verbose = FALSE)
  )
})

# Step 4: Joint estimation with MV-FMR
result <- mvfmr(
  G = sim_data$details$G,
  fpca_results = fpca_results,
  Y = outcome_data$Y,
  outcome_type = "continuous",
  method = "gmm",
  max_nPC = c(10, 10),
  bootstrap = TRUE,
  n_bootstrap = 100
)

# View results
print(result)
summary(result)
plot(result)  # Visualize time-varying effects for every exposure

# Extract coefficients and effects
coef(result)
result$effects[[1]]  # Time-varying effect for exposure 1
result$effects[[2]]  # Time-varying effect for exposure 2

# Separate estimation (U-FMR for each exposure independently), reusing
# sim_data / outcome_data / fpca_results simulated in Example 1
result_separate <- mvfmr_separate(
  G_list = list(sim_data$details$G, sim_data$details$G),
  fpca_results = fpca_results,
  Y = outcome_data$Y
)

# Compare performance: joint (`result`, from Example 1) vs. separate
result$performance
result_separate$exposures[[1]]$performance
result_separate$exposures[[2]]$performance

library(mvfmr)
library(fdapace)

# Step 1: Simulate a single exposure
set.seed(12345)
sim_data <- getX_multi_exposure(
  N = 1000,
  J = 50,
  nSparse = 10,
  n_exposures = 1
)

# Step 2: Generate outcome
outcome_data <- getY_multi_exposure(
  sim_data,
  XYmodels = "2",     # Linear effect
  X_effects = TRUE,
  outcome_type = "continuous"
)

# Step 3: FPCA for the (only) exposure
fpca1 <- FPCA(
  sim_data$exposures[[1]]$Ly_sim,
  sim_data$exposures[[1]]$Lt_sim,
  list(dataType = 'Sparse', error = TRUE, verbose = FALSE)
)

# Step 4: Univariable estimation
result <- mvfmr_separate(
  G_list = list(sim_data$details$G),  # A list of length 1
  fpca_results = list(fpca1),
  Y = outcome_data$Y,
  outcome_type = "continuous",
  method = "gmm",
  max_nPC = 10
)

# View results
print(result)
coef(result, exposure = 1)
result$exposures[[1]]$effect  # Time-varying effect

set.seed(2026)
sim_data3 <- getX_multi_exposure(N = 1000, J = 50, nSparse = 10, n_exposures = 3)

outcome_data3 <- getY_multi_exposure(
  sim_data3,
  XYmodels = c("2", "5", "8"),
  outcome_type = "continuous"
)

fpca_results3 <- lapply(sim_data3$exposures, function(exp_k) {
  FPCA(exp_k$Ly_sim, exp_k$Lt_sim, list(dataType = 'Sparse', error = TRUE, verbose = FALSE))
})

# Joint estimation across all 3 exposures
result_joint3 <- mvfmr(
  G = sim_data3$details$G,
  fpca_results = fpca_results3,
  Y = outcome_data3$Y,
  outcome_type = "continuous",
  method = "gmm",
  max_nPC = c(10, 10, 10),
  true_effects = c("2", "5", "8"),
  X_true = sim_data3$details$X_list
)

print(result_joint3)
plot(result_joint3)  # One panel per exposure

# Separate estimation across all 3 exposures
result_separate3 <- mvfmr_separate(
  G_list = list(sim_data3$details$G, sim_data3$details$G, sim_data3$details$G),
  fpca_results = fpca_results3,
  Y = outcome_data3$Y,
  max_nPC = c(10, 10, 10),
  true_effects = c("2", "5", "8")
)

# Access any exposure by index (1..m)
result_joint3$effects[[3]]
coef(result_separate3, exposure = 3)

library(mvfmr)
library(fdapace)

# Step 1: Simulate exposure data (individual-level)
set.seed(12345)
sim_data <- getX_multi_exposure(
  N = 5000,           # Exposure sample size
  J = 30,             # Number of genetic instruments (SNPs)
  nSparse = 10,
  n_exposures = 2
)

# Perform FPCA on longitudinal exposures
fpca_results <- lapply(sim_data$exposures, function(exp_k) {
  FPCA(exp_k$Ly_sim, exp_k$Lt_sim, list(dataType = 'Sparse', error = TRUE, verbose = FALSE))
})

# Step 2: Get outcome GWAS summary statistics (from a separate study)
# Simulate obtaining summary statistics from a separate GWAS
# (this mimics what you'd get from a published GWAS)
by_outcome <- rnorm(30, mean = 0.02, sd = 0.01)  # SNP-outcome associations
sy_outcome <- runif(30, 0.005, 0.015)            # Standard errors
ny_outcome <- 100000                             # GWAS sample size

# Step 3: Two-sample MV-FMR estimation
result_twosample <- fmvmr_twosample(
  G_exposure = sim_data$details$G,   # Genotypes from the exposure sample
  fpca_results = fpca_results,       # FPCA from the exposures
  by_outcome = by_outcome,           # GWAS betas (from the outcome study)
  sy_outcome = sy_outcome,           # GWAS standard errors
  ny_outcome = ny_outcome,           # GWAS sample size
  max_nPC = c(3, 3),
  verbose = TRUE
)

# Step 4: View results
print(result_twosample)

# Extract time-varying effects
result_twosample$effects[[1]]
result_twosample$effects[[2]]

getX_multi_exposure(
  N = 1000,                  # Sample size
  J = 50,                    # Number of genetic instruments
  nSparse = 10,               # Observations per subject
  n_exposures = 2,            # Number of exposures (m)
  shared_effect = TRUE,       # Whether all exposures share the same time-varying confounding
  separate_G = FALSE,         # Whether to use separate instruments per exposure
  shared_G_proportion = 0.15  # Proportion of shared instruments (0-1, if separate_G = TRUE)
)

getX_multi_exposure_mediation(
  N = 1000,                  # Sample size
  J = 50,                    # Number of genetic instruments
  nSparse = 10,               # Observations per subject
  n_exposures = 2,            # Number of exposures (m)
  mediation_strength = NULL,  # m x m matrix: entry [j, k] (j < k) is the strength
                              # with which exposure j mediates its effect onto
                              # exposure k. Default: NULL = no mediation.
  mediation_type = "linear"   # "linear", "nonlinear", "time_varying" (scalar or
                              # m x m matrix mirroring mediation_strength)
)

getY_multi_exposure(
  RES,                         # Output from getX_multi_exposure() or getX_multi_exposure_mediation()
  XYmodels = NULL,             # Length-m vector of effect models, one per exposure (see below); default '1' for all
  X_effects = NULL,            # Length-m logical vector: include each exposure's effect?; default TRUE for all
  outcome_type = "continuous"  # "continuous" or "binary"
)

mvfmr(
  G,                                     # Genetic instrument matrix (N x J)
  fpca_results,                          # List of length m of FPCA objects, one per exposure
  Y,                                     # Outcome vector
  outcome_type = "continuous",           # "continuous" or "binary"
  method = "gmm",                        # "gmm", "cf" (control function), or "cf-lasso"
  nPC = NA,                              # Fixed number of components per exposure (length 1 or m; NA = select automatically)
  max_nPC = NA,                          # Maximum number of components per exposure (length 1 or m)
  improvement_threshold = 0.001,         # Minimum CV improvement required to add a component
  bootstrap = FALSE,                     # Whether to compute bootstrap confidence intervals
  n_bootstrap = 100,                     # Number of bootstrap replicates
  n_cores = parallel::detectCores() - 1, # Number of CPU cores for parallel computations
  true_effects = NULL,                   # Length-m vector of true effect model codes (simulation only)
  X_true = NULL,                         # Length-m list of true X curves (simulation only)
  verbose = FALSE                        # Print progress and diagnostic messages
)

mvfmr_separate(
  G_list,                # List of length m of genetic instrument matrices, one per exposure
                         # (use a list of length 1 to analyze a single exposure)
  fpca_results,          # List of length m of FPCA objects, same length as G_list
  Y,                     # Outcome vector
  outcome_type = "continuous",
  method = "gmm",
  nPC = NA,
  max_nPC = NA,
  improvement_threshold = 0.001,
  bootstrap = FALSE,
  n_bootstrap = 100,
  n_cores = parallel::detectCores() - 1,
  true_effects = NULL,
  X_true = NULL,
  verbose = FALSE
)

fmvmr_twosample(
  G_exposure,            # Genetic instrument matrix from the exposure sample (N x J)
  fpca_results,          # List of length m of FPCA objects
  by_outcome,            # Vector of SNP-outcome betas from the outcome GWAS, length J
  sy_outcome,            # Vector of standard errors for SNP-outcome effects, length J
  ny_outcome,            # Sample size of the outcome GWAS
  max_nPC = NA,          # Maximum number of components per exposure (length 1 or m)
  true_effects = NULL,   # Length-m vector of true effect model codes (simulation only)
  verbose = TRUE
)

fmvmr_separate_twosample(
  G_list,                # List of length m of genetic instrument matrices
  fpca_results,          # List of length m of FPCA objects
  by_outcome_list,       # List of length m of SNP-outcome beta vectors
  sy_outcome_list,       # List of length m of SNP-outcome standard error vectors
  ny_outcome,            # Outcome GWAS sample size
  max_nPC = NA,
  true_effects = NULL,
  verbose = TRUE
)

IS(
  J,                     # Number of genetic instruments
  K,                     # Number of exposures/components
  PC,                    # Vector of indices indicating which columns in datafull correspond to the principal components
  datafull,              # Data frame containing instruments (first J columns) and principal components (subsequent columns) [G, X]
  Y                      # Optional outcome vector; if provided, Q-statistic for overidentification is calculated
)

result <- mvfmr(...)
names(result)

result <- mvfmr_separate(...)
names(result)

# Generate binary outcome
outcome_binary <- getY_multi_exposure(
  sim_data,
  XYmodels = c("2", "8"),
  outcome_type = "binary"
)

# Estimate with control function
result <- mvfmr(
  G = sim_data$details$G,
  fpca_results = list(fpca1, fpca2),
  Y = outcome_binary$Y,
  outcome_type = "binary",
  method = "cf"
)

result <- mvfmr(
  G = G,
  fpca_results = list(fpca1, fpca2),
  Y = Y,
  max_nPC = c(10, 10),          # Search up to 10 components per exposure
  improvement_threshold = 0.01  # Stop if improvement < 1%
)

# View selected components
result$nPC_used

result <- mvfmr(
  G = G,
  fpca_results = list(fpca1, fpca2),
  Y = Y,
  bootstrap = TRUE,
  n_bootstrap = 200  # Number of bootstrap replicates
)

# Bootstrap confidence intervals available in:
result$confidence_intervals

result <- mvfmr(
  G = G,
  fpca_results = list(fpca1, fpca2),
  Y = Y,
  n_cores = 4  # Use 4 cores for cross-validation
)

# Generate data where exposure 1 mediates onto exposure 2
mediation_strength <- matrix(0, 2, 2)
mediation_strength[1, 2] <- 0.5

sim_mediation <- getX_multi_exposure_mediation(
  N = 1000,
  J = 50,
  n_exposures = 2,
  mediation_strength = mediation_strength,
  mediation_type = "linear"
)

outcome <- getY_multi_exposure(
  sim_mediation,
  XYmodels = c("2", "1"),  # Direct effect of exposure 1; effect of exposure 2 (mediator)
  outcome_type = "continuous"
)

fpca_results <- lapply(sim_mediation$exposures, function(exp_k) {
  FPCA(exp_k$Ly_sim, exp_k$Lt_sim, list(dataType = 'Sparse', error = TRUE, verbose = FALSE))
})

# Estimate with MV-FMR to capture mediation
result <- mvfmr(
  G = sim_mediation$details$G,
  fpca_results = fpca_results,
  Y = outcome$Y
)

# After FPCA
K_total <- sum(sapply(fpca_results, function(f) f$selectK))

PC_stacked <- do.call(cbind, lapply(fpca_results, function(f) f$xiEst[, 1:f$selectK]))

fstats <- IS(
  J = ncol(G),
  K = K_total,
  PC = 1:K_total,
  datafull = cbind(G, PC_stacked)
)

# View conditional F-statistics (cFF)
print(fstats)