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: 25Step 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.65Step 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 selectedStep 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.745Step 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.03681268 0.04977623 0.06537926Step 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
