Example: Joint Estimation of Two Exposures
Step 1: Simulate Data
We generate genetic instruments and two longitudinal exposures:
set.seed(12345)
# Generate exposure data
sim_data <- getX_multi_exposure(
N = 300, # Sample size
J = 25, # Number of genetic instruments
nSparse = 10 # Observations per subject
)
# 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 X1 and quadratic effect from X2:
outcome_data <- getY_multi_exposure(
sim_data,
X1Ymodel = "2", # Linear: β₁(t) = 0.02 * t
X2Ymodel = "5", # Late-age effect: β₂(t) = 0.1 * (t > 30)
X1_effect = TRUE,
X2_effect = TRUE,
outcome_type = "continuous"
)
cat("Outcome summary:\n")
#> Outcome summary:
summary(outcome_data$Y)
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> -111.070 -31.849 -5.190 -3.012 26.394 112.621Step 3: Functional Principal Component Analysis
We apply FPCA to both exposures to extract principal components:
# FPCA for exposure 1
fpca1 <- FPCA(
sim_data$X1$Ly_sim,
sim_data$X1$Lt_sim,
list(dataType = 'Sparse', error = TRUE, verbose = FALSE)
)
# FPCA for exposure 2
fpca2 <- FPCA(
sim_data$X2$Ly_sim,
sim_data$X2$Lt_sim,
list(dataType = 'Sparse', error = TRUE, verbose = FALSE)
)
cat("FPCA completed:\n")
#> FPCA completed:
cat(" Exposure 1:", fpca1$selectK, "components selected\n")
#> Exposure 1: 3 components selected
cat(" Exposure 2:", fpca2$selectK, "components selected\n")
#> Exposure 2: 4 components selectedStep 4: Joint Multivariable Estimation
Now we perform joint estimation using mvfmr():
result_joint <- mvfmr(
G = sim_data$details$G,
fpca_results = list(fpca1, fpca2),
Y = outcome_data$Y,
outcome_type = "continuous",
method = "gmm",
max_nPC1 = 4,
max_nPC2 = 4,
improvement_threshold = 0.001,
bootstrap = FALSE,
n_cores = 1,
true_effects = list(model1 = "2", model2 = "5"),
X_true = list(X1_true = sim_data$details$X1, X2_true = sim_data$details$X2),
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.001899 , Coverage: 0.824
#> Exposure 2 - MISE: 0.001591 , Coverage: 0.922Step 5: Visualize Time-Varying Effects
The estimated time-varying causal effects can be visualized using the built-in plot method:
# Plot both effects
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] -3.7775154 1.5166093 -0.3552057 0.4624341
# Time-varying effects at each time point
head(result_joint$effects$effect1)
#> [1] 0.02665083 0.03882021 0.05180349 0.06570553 0.08061147 0.09654441
head(result_joint$effects$effect2)
#> [1] -0.04923079 -0.05503448 -0.06005135 -0.06404133 -0.06679739 -0.06817082Step 7: Performance Metrics
Since we simulated data with known truth, we can evaluate performance:
cat("Performance Metrics:\n")
#> Performance Metrics:
cat("\nExposure 1 (Linear effect):\n")
#>
#> Exposure 1 (Linear effect):
cat(" MISE:", round(result_joint$performance$MISE1, 6), "\n")
#> MISE: 0.001899
cat(" Coverage:", round(result_joint$performance$Coverage1, 3), "\n")
#> Coverage: 0.824
cat("\nExposure 2 (Quadratic effect):\n")
#>
#> Exposure 2 (Quadratic effect):
cat(" MISE:", round(result_joint$performance$MISE2, 6), "\n")
#> MISE: 0.001591
cat(" Coverage:", round(result_joint$performance$Coverage2, 3), "\n")
#> Coverage: 0.922- MISE (Mean Integrated Squared Error): measures accuracy of effect estimation
- Coverage: proportion of time points where true effect falls within 95% CI