Introduction to Multivariable Functional Mendelian Randomization

Overview

This vignette introduces Multivariable Functional Mendelian Randomization (MV-FMR), a method to estimate time-varying causal effects of multiple correlated longitudinal exposures on health outcomes.

Key Features

Joint estimation of an arbitrary number of time-varying exposures
Automatic component selection via cross-validation
Support for continuous and binary outcomes
Bootstrap inference for confidence intervals
Instrument strength diagnostics

When to Use MV-FMR

Use joint multivariable estimation (mvfmr()) when:

You have multiple exposures that may be correlated
Exposures may share horizontal pleiotropy (common genetic variants)
There may be mediation between exposures
You want to estimate independent causal effects of each exposure

Installation

# Install from CRAN
install.packages("mvfmr")

library(mvfmr)
library(fdapace)
library(ggplot2)

Example: Joint Estimation of Two Exposures

Step 1: Simulate Data

We generate genetic instruments and two longitudinal exposures:

set.seed(473920)

# Generate exposure data
sim_data <- getX_multi_exposure(
  N = 300,              # Sample size
  J = 25,               # Number of genetic instruments
  nSparse = 10,          # Observations per subject
  n_exposures = 2        # Number of exposures (m)
)

# Check dimensions
cat("Sample size:", nrow(sim_data$details$G), "\n")
#> Sample size: 300
cat("Number of instruments:", ncol(sim_data$details$G), "\n")
#> Number of instruments: 25

Step 2: Generate Outcome

We create an outcome with linear effect from exposure 1 and exposure 2:

outcome_data <- getY_multi_exposure(
  sim_data,
  XYmodels = c("2", "2"),     # Exposure 1/2: linear beta(t) = 0.02*t
  X_effects = c(TRUE, TRUE),
  outcome_type = "continuous"
)

cat("Outcome summary:\n")
#> Outcome summary:
summary(outcome_data$Y)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#> -236.50  -71.40  -26.03  -23.83   28.53  182.65

Step 3: Functional Principal Component Analysis

We apply FPCA to each exposure to extract principal components:

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

cat("FPCA completed:\n")
#> FPCA completed:
for (k in seq_along(fpca_results)) {
  cat("  Exposure", k, ":", fpca_results[[k]]$selectK, "components selected\n")
}
#>   Exposure 1 : 3 components selected
#>   Exposure 2 : 3 components selected

Step 4: Joint Multivariable Estimation

Now we perform joint estimation using mvfmr(). All exposure-related arguments (fpca_results, max_nPC, true_effects, X_true, …) are lists or vectors of length m, one entry per exposure:

result_joint <- mvfmr(
  G = sim_data$details$G,
  fpca_results = fpca_results,
  Y = outcome_data$Y,
  outcome_type = "continuous",
  method = "gmm",
  max_nPC = c(4, 4),
  improvement_threshold = 0.001,
  bootstrap = FALSE,
  n_cores = 1,
  true_effects = c("2", "2"),
  X_true = sim_data$details$X_list,
  verbose = FALSE
)

# View results
print(result_joint)
#> 
#> Functional Multivariable MR Result
#> ==================================
#> Exposures: 2 
#> Sample size: 300 
#> Outcome: continuous 
#> Method: gmm 
#> Components selected: nPC1 = 2, nPC2 = 2 
#> 
#> Performance Metrics:
#>   Exposure 1 - MISE: 0.002458 , Coverage: 1 
#>   Exposure 2 - MISE: 0.003881 , Coverage: 0.745

Step 5: Visualize Time-Varying Effects

The estimated time-varying causal effects can be visualized using the built-in plot method:

# Plot every exposure's effect
plot(result_joint)

The solid colored lines show the estimated time-varying causal effects, with shaded bands representing 95% confidence intervals. The dashed red lines (when present) indicate the true time-varying effects used in the simulation.

Step 6: Extract Coefficients

# Estimated beta coefficients for basis functions
coef(result_joint)
#> [1] -4.035692 -1.647212 -3.634343 -1.451868

# Time-varying effects at each time point (one entry per exposure)
head(result_joint$effects[[1]])
#> [1] 0.06182589 0.07321675 0.08677954 0.10220947 0.11906258 0.13675702
head(result_joint$effects[[2]])
#> [1] 0.02302308 0.02248962 0.02740122 0.03681267 0.04977623 0.06537926

Step 7: Performance Metrics

Since we simulated data with known truth, we can evaluate performance:

cat("Performance Metrics:\n")
#> Performance Metrics:
for (k in seq_along(result_joint$effects)) {
  cat("\nExposure", k, ":\n")
  cat("  MISE:", round(result_joint$performance$MISE[[k]], 6), "\n")
  cat("  Coverage:", round(result_joint$performance$Coverage[[k]], 3), "\n")
}
#> 
#> Exposure 1 :
#>   MISE: 0.002458 
#>   Coverage: 1 
#> 
#> Exposure 2 :
#>   MISE: 0.003881 
#>   Coverage: 0.745

MISE (Mean Integrated Squared Error): measures accuracy of effect estimation
Coverage: proportion of time points where true effect falls within 95% CI

Comparison: Joint vs. Separate Estimation

We can compare joint (multivariable) with separate (univariable) estimation. For mvfmr_separate(), instruments are passed as G_list, a list of length m (here the same shared instrument matrix is reused for both exposures):

result_separate <- mvfmr_separate(
  G_list = list(sim_data$details$G, sim_data$details$G),
  fpca_results = fpca_results,
  Y = outcome_data$Y,
  outcome_type = "continuous",
  method = "gmm",
  max_nPC = c(4, 4),
  n_cores = 1,
  true_effects = c("2", "2"),
  verbose = FALSE
)
#> [1] "Processing X1"
#> [1] "Processing X2"

print(result_separate)
#> 
#> Separate Univariable MR Results
#> ================================
#> Exposures: 2 
#> Separate instruments: TRUE 
#> Outcome: continuous 
#> Method: gmm 
#> 
#> Exposure 1 :
#>   Components: 2 
#>   MSE: 0.178885 
#>   Coverage: 0 
#> 
#> Exposure 2 :
#>   Components: 2 
#>   MSE: 0.252422 
#>   Coverage: 0.157

Performance Comparison

comparison <- data.frame(
  Method = rep(c("Joint (MV-FMR)", "Separate (U-FMR)"), each = 2),
  Exposure = rep(c("X1", "X2"), times = 2),
  MISE = c(
    result_joint$performance$MISE[[1]],
    result_joint$performance$MISE[[2]],
    result_separate$exposures[[1]]$performance$MISE,
    result_separate$exposures[[2]]$performance$MISE
  ),
  Coverage = c(
    result_joint$performance$Coverage[[1]],
    result_joint$performance$Coverage[[2]],
    result_separate$exposures[[1]]$performance$Coverage,
    result_separate$exposures[[2]]$performance$Coverage
  )
)

print(comparison)
#>             Method Exposure        MISE  Coverage
#> 1   Joint (MV-FMR)       X1 0.002458118 1.0000000
#> 2   Joint (MV-FMR)       X2 0.003881328 0.7450980
#> 3 Separate (U-FMR)       X1 0.178885369 0.0000000
#> 4 Separate (U-FMR)       X2 0.252422123 0.1568627

Instrument Strength Diagnostics

Check instrument strength using F-statistics. IS() is already generic in the number of exposures/components K:

# Calculate F-statistics
K_total <- sum(result_joint$nPC_used)

PC_stacked <- do.call(cbind, lapply(seq_along(fpca_results), function(k) {
  fpca_results[[k]]$xiEst[, 1:result_joint$nPC_used[k]]
}))

fstats <- IS(
  J = ncol(sim_data$details$G),
  K = K_total,
  PC = 1:K_total,
  datafull = cbind(sim_data$details$G, PC_stacked),
  Y = outcome_data$Y
)

fstats_df <- cbind(
  "Exposure" = unlist(lapply(seq_along(result_joint$nPC_used), function(k) {
    rep(paste0("X", k), result_joint$nPC_used[k])
  })),
  as.data.frame(fstats)
)

print(fstats_df[, c("Exposure", "PC", "cFF")])
#>   Exposure PC       cFF
#> 1       X1  1 0.3913644
#> 2       X1  2 1.1379411
#> 3       X2  3 0.3986037
#> 4       X2  4 0.9953631

FF: Standard F-statistic for instrument strength
cFF: Conditional F-statistic (accounts for other exposures)
Qvalue: Hansen’s J overidentification test statistic
pvalue: P-value for Q-test (high p-value = instruments are valid)

Rule of thumb: cFF > 10 indicates strong instruments.

Binary Outcomes

MV-FMR also supports binary outcomes using the control function approach:

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

# Estimate with control function
result_binary <- mvfmr(
  G = sim_data$details$G,
  fpca_results = fpca_results,
  Y = outcome_binary$Y,
  outcome_type = "binary",
  method = "cf",  # Control function for binary
  max_nPC = c(3, 3),
  n_cores = 1,
  verbose = FALSE
)

print(result_binary)

Extending to More Than Two Exposures

Nothing changes in the API when moving from 2 to m exposures: fpca_results, max_nPC, true_effects and X_true simply grow to length m.

set.seed(163918)#281046
sim_data3 <- getX_multi_exposure(N = 500, J = 50, nSparse = 10, n_exposures = 3)

outcome_data3 <- getY_multi_exposure(
  sim_data3,
  XYmodels = c("2", "2", "2"),
  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))
})

result_joint3 <- mvfmr(
  G = sim_data3$details$G,
  fpca_results = fpca_results3,
  Y = outcome_data3$Y,
  outcome_type = "continuous",
  method = "gmm",
  max_nPC = c(4, 4, 4),
  n_cores = 1,
  true_effects = c("2", "2", "2"),
  X_true = sim_data3$details$X_list,
  verbose = FALSE
)

print(result_joint3)
#> 
#> Functional Multivariable MR Result
#> ==================================
#> Exposures: 3 
#> Sample size: 500 
#> Outcome: continuous 
#> Method: gmm 
#> Components selected: nPC1 = 2, nPC2 = 2, nPC3 = 2 
#> 
#> Performance Metrics:
#>   Exposure 1 - MISE: 0.001228 , Coverage: 1 
#>   Exposure 2 - MISE: 0.004026 , Coverage: 0.706 
#>   Exposure 3 - MISE: 0.006513 , Coverage: 0.588

plot(result_joint3)

Next Steps

Learn More

Univariable FMR: See vignette("univariable-fmr") for single exposure analysis
Advanced examples: Check inst/examples/test_MV-FMR.R for complete scenarios
Manuscript: See paper for methodological details

Citation

If you use this package, please cite:

Fontana, N., Ieva, F., Zuccolo, L., Di Angelantonio, E., & Secchi, P. (2025). Unraveling time-varying causal effects of multiple exposures: integrating Functional Data Analysis with Multivariable Mendelian Randomization. arXiv preprint arXiv:2512.19064.

Session Info

sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04.1 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] ggplot2_4.0.3  fdapace_0.6.0  mvfmr_0.2.0    rmarkdown_2.32
#> 
#> loaded via a namespace (and not attached):
#>  [1] gtable_0.3.6        shape_1.4.6.1       xfun_0.61          
#>  [4] bslib_0.12.0        htmlwidgets_1.6.4   lattice_0.23-1     
#>  [7] numDeriv_2016.8-1.1 vctrs_0.7.3         tools_4.6.1        
#> [10] generics_0.1.4      parallel_4.6.1      tibble_3.3.1       
#> [13] cluster_2.1.8.3     pkgconfig_2.0.3     Matrix_1.7-6       
#> [16] data.table_1.18.6.1 checkmate_2.3.4     RColorBrewer_1.1-3 
#> [19] S7_0.2.2            lifecycle_1.0.5     compiler_4.6.1     
#> [22] farver_2.1.2        stringr_1.6.0       progress_1.2.3     
#> [25] codetools_0.2-20    htmltools_0.5.9     sys_3.4.3          
#> [28] buildtools_1.0.0    sass_0.4.10         yaml_2.3.12        
#> [31] glmnet_5.0          htmlTable_2.5.0     Formula_1.2-6      
#> [34] pracma_2.4.6        crayon_1.5.3        pillar_1.11.1      
#> [37] jquerylib_0.1.4     MASS_7.3-66         cachem_1.1.0       
#> [40] Hmisc_5.3-0         iterators_1.0.14    rpart_4.1.27       
#> [43] foreach_1.5.2       tidyselect_1.2.1    digest_0.6.39      
#> [46] stringi_1.8.9       dplyr_1.2.1         labeling_0.4.3     
#> [49] maketools_1.3.2     splines_4.6.1       fastmap_1.2.0      
#> [52] grid_4.6.1          colorspace_2.1-3    cli_3.6.6          
#> [55] magrittr_2.0.5      base64enc_0.1-6     survival_3.8-12    
#> [58] withr_3.0.3         foreign_0.8-91      prettyunits_1.2.0  
#> [61] scales_1.4.0        backports_1.5.1     otel_0.2.0         
#> [64] nnet_7.3-21         gridExtra_2.3.1     hms_1.1.4          
#> [67] evaluate_1.0.5      knitr_1.52          doParallel_1.0.17  
#> [70] rlang_1.3.0         Rcpp_1.1.2          glue_1.8.1         
#> [73] pROC_1.19.1         rstudioapi_0.19.0   jsonlite_2.0.0     
#> [76] R6_2.6.1