Univariable Functional Mendelian Randomization

Example: Single Exposure Analysis

Step 1: Simulate Data

set.seed(12345)

# Generate exposure data with a single exposure (m = 1)
sim_data <- getX_multi_exposure(
  N = 300,
  J = 25,
  nSparse = 10,
  n_exposures = 1
)

cat("Data simulated:\n")
#> Data simulated:
cat("  Sample size:", nrow(sim_data$details$G), "\n")
#>   Sample size: 300
cat("  Instruments:", ncol(sim_data$details$G), "\n")
#>   Instruments: 25

Step 2: Generate Outcome

We create an outcome where the exposure has a linear causal effect:

outcome_data <- getY_multi_exposure(
  sim_data,
  XYmodels = "2",     # Linear effect: beta(t) = 0.02*t
  X_effects = TRUE,
  outcome_type = "continuous"
)

cat("Outcome summary:\n")
#> Outcome summary:
summary(outcome_data$Y)
#>     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
#> -106.837  -30.302   -3.945   -3.029   24.274  111.770

Step 3: FPCA for the Single Exposure

fpca1 <- FPCA(
  sim_data$exposures[[1]]$Ly_sim,
  sim_data$exposures[[1]]$Lt_sim,
  list(dataType = 'Sparse', error = TRUE, verbose = FALSE)
)

cat("FPCA completed:\n")
#> FPCA completed:
cat("  Components selected:", fpca1$selectK, "\n")
#>   Components selected: 2
cat("  Variance explained:",
    round(sum(fpca1$lambda[1:fpca1$selectK]) / sum(fpca1$lambda) * 100, 1), "%\n")
#>   Variance explained: 100 %

Step 4: Univariable Estimation

Estimate the causal effect using mvfmr_separate() with a G_list of length 1:

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

print(result)
#> 
#> Separate Univariable MR Results
#> ================================
#> Exposures: 1 
#> Separate instruments: TRUE 
#> Outcome: continuous 
#> Method: gmm 
#> 
#> Exposure 1 :
#>   Components: 2 
#>   MSE: 0.002314 
#>   Coverage: 0.745

Step 5: Visualize Time-Varying Effect

plot(result)

The solid colored line shows the estimated time-varying causal effect, with a shaded band representing the 95% confidence interval. The dashed red line (when present) indicates the true time-varying effect used in the simulation.

Step 6: Extract Results

# Coefficients for basis functions
cat("Estimated coefficients:\n")
#> Estimated coefficients:
print(round(coef(result, exposure = 1), 4))
#> [1] -4.0138  1.1379

# Time-varying effect curve
cat("\nFirst 10 time points of beta(t):\n")
#> 
#> First 10 time points of beta(t):
head(result$exposures[[1]]$effect, 10)
#>  [1] 0.1332247 0.1355505 0.1394325 0.1448319 0.1517528 0.1602287 0.1702996
#>  [8] 0.1819805 0.1952264 0.2098995

Step 7: Performance Metrics

cat("Performance:\n")
#> Performance:
cat("  MISE:", round(result$exposures[[1]]$performance$MISE, 6), "\n")
#>   MISE: 0.002314
cat("  Coverage:", round(result$exposures[[1]]$performance$Coverage, 3), "\n")
#>   Coverage: 0.745
cat("  Components used:", result$exposures[[1]]$nPC_used, "\n")
#>   Components used: 2

Binary Outcomes

U-FMR also works with binary outcomes:

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

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

print(result_binary)
cat("Cases:", sum(outcome_binary$Y == 1), "\n")
cat("Controls:", sum(outcome_binary$Y == 0), "\n")

Advanced Topics

Available Effect Models

The package includes 10 pre-defined time-varying effect shapes:

"0": No effect (null)
"1": Constant (beta = 0.1)
"2": Linear increasing
"3": Linear decreasing
"4": Early life effect
"5": Late life effect
"6": Early decreasing
"7": Late increasing
"8": Quadratic (U-shape)
"9": Cubic

Bootstrap Inference

# Get robust confidence intervals via bootstrap
result_boot <- mvfmr_separate(
  G_list = list(sim_data$details$G),
  fpca_results = list(fpca1),
  Y = outcome_data$Y,
  outcome_type = "continuous",
  bootstrap = TRUE,
  n_bootstrap = 100,
  max_nPC = 4,
  verbose = FALSE
)

# Bootstrap CIs are stored in result_boot$exposures[[1]]$...

Two-Sample Design

If you have GWAS summary statistics instead of individual-level outcome data:

# Simulate GWAS summary statistics
by_outcome <- rnorm(25, 0.02, 0.01)
sy_outcome <- runif(25, 0.005, 0.015)

result_2sample <- fmvmr_separate_twosample(
  G_list = list(sim_data$details$G),
  fpca_results = list(fpca1),
  by_outcome_list = list(by_outcome),
  sy_outcome_list = list(sy_outcome),
  ny_outcome = 50000,
  max_nPC = 3,
  verbose = FALSE
)

print(result_2sample)

Next Steps

Learn More

Multivariable FMR: See vignette("multivariable-fmr") for joint estimation
Complete examples: Check inst/examples/test_U-FMR.R for all 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.3.0 (2023-04-21)
#> Platform: aarch64-apple-darwin20 (64-bit)
#> Running under: macOS Ventura 13.6
#> 
#> Matrix products: default
#> BLAS:   /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRblas.0.dylib 
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.11.0
#> 
#> locale:
#> [1] C/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#> 
#> time zone: Europe/Rome
#> tzcode source: internal
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] ggplot2_3.5.2 fdapace_0.5.9 mvfmr_0.2.0  
#> 
#> loaded via a namespace (and not attached):
#>  [1] gtable_0.3.6        shape_1.4.6         xfun_0.41          
#>  [4] bslib_0.6.1         htmlwidgets_1.6.4   lattice_0.21-8     
#>  [7] numDeriv_2016.8-1.1 vctrs_0.6.5         tools_4.3.0        
#> [10] generics_0.1.3      parallel_4.3.0      tibble_3.2.1       
#> [13] fansi_1.0.6         highr_0.10          cluster_2.1.4      
#> [16] pkgconfig_2.0.3     Matrix_1.6-4        data.table_1.14.10 
#> [19] checkmate_2.3.1     lifecycle_1.0.4     farver_2.1.1       
#> [22] compiler_4.3.0      stringr_1.5.1       progress_1.2.3     
#> [25] munsell_0.5.0       codetools_0.2-19    htmltools_0.5.7    
#> [28] sass_0.4.8          yaml_2.3.8          glmnet_4.1-8       
#> [31] htmlTable_2.4.2     Formula_1.2-5       pracma_2.4.4       
#> [34] pillar_1.9.0        crayon_1.5.3        jquerylib_0.1.4    
#> [37] MASS_7.3-58.4       cachem_1.0.8        Hmisc_5.1-1        
#> [40] iterators_1.0.14    rpart_4.1.19        foreach_1.5.2      
#> [43] tidyselect_1.2.0    digest_0.6.33       stringi_1.8.3      
#> [46] dplyr_1.1.4         labeling_0.4.3      splines_4.3.0      
#> [49] fastmap_1.1.1       grid_4.3.0          colorspace_2.1-0   
#> [52] cli_3.6.2           magrittr_2.0.3      base64enc_0.1-3    
#> [55] survival_3.5-7      utf8_1.2.4          withr_2.5.2        
#> [58] foreign_0.8-84      prettyunits_1.2.0   scales_1.3.0       
#> [61] backports_1.5.0     rmarkdown_2.25      nnet_7.3-18        
#> [64] gridExtra_2.3       hms_1.1.3           evaluate_0.23      
#> [67] knitr_1.45          doParallel_1.0.17   rlang_1.1.6        
#> [70] Rcpp_1.1.0          glue_1.6.2          pROC_1.18.5        
#> [73] rstudioapi_0.15.0   jsonlite_1.8.8      R6_2.5.1           
#> [76] plyr_1.8.9