Section 5.3: Multivariate Endogenous Variables

library(divR)
library(ggplot2)

Introduction

In Section 5.3 of the paper, we demonstrate the capability of DIV to handle multivariate endogenous variables. In many real-world applications, several treatment variables may be affected by the same unobserved confounding factors. DIV’s generator-based architecture naturally extends to such cases.

Data Generating Process

We consider a setting with two endogenous variables, \(X_1\) and \(X_2\), and a continuous instrument \(Z\): - \(X_1 = \text{log}(2Z + H + \epsilon_{X1} + 7)\) - \(X_2 = Z \cdot (2H - 0.5\epsilon_{X2})\) - \(Y = X_1 + 2X_2 + 2H + \epsilon_Y\)

Here, both \(X_1\) and \(X_2\) are functions of the same hidden confounder \(H\), which also affects the outcome \(Y\).

g_log <- function(Z, H, eps_X) log(2 * Z + H + eps_X + 7)
g_mult <- function(Z, H, eps_X) Z * (2 * H - 0.5 * eps_X)
f_lin <- function(X1, X2, H, eps_Y) X1 + 2*X2 + 2*H + eps_Y

set.seed(1998)
n <- 2000
H <- rnorm(n)
Z <- runif(n, 0, 3)
eps_X1 <- rnorm(n); eps_X2 <- rnorm(n); eps_Y <- rnorm(n)

X1 <- g_log(Z, H, eps_X1)
X2 <- g_mult(Z, H, eps_X2)
Y <- f_lin(X1, X2, H, eps_Y)

X <- cbind(X1, X2)

Fit divR Model

We fit the model using the matrix of endogenous variables \(X\).

model <- divR(Z = matrix(Z, ncol=1),
              X = X,
              Y = matrix(Y, ncol=1),
              num_epochs = 2000,
              silent = TRUE)

Partial Interventional Mean

To visualize the results, we examine the Partial Interventional Mean. We fix one variable (\(X_2\)) at its median value and vary the other (\(X_1\)) across its range.

\[E[Y | do(X_1 = x_1, X_2 = \text{median}(X_2))]\]

x2_fixed <- median(X2)
x1_grid <- seq(min(X1), max(X1), length.out = 100)
X_test <- cbind(x1_grid, rep(x2_fixed, 100))

Y_causal_pred <- predict(model, Xtest = X_test, type = "mean", nsample = 500)

# True partial mean via simulation
get_true_mean <- function(x1, x2, n_rep = 5000) {
  H_rep <- rnorm(n_rep); eps_Y_rep <- rnorm(n_rep)
  Y_rep <- f_lin(rep(x1, n_rep), rep(x2, n_rep), H_rep, eps_Y_rep)
  mean(Y_rep)
}
true_means <- sapply(x1_grid, function(x1) get_true_mean(x1, x2_fixed))

plot_df <- data.frame(
  X1 = x1_grid,
  Predicted = as.vector(Y_causal_pred),
  True = true_means
)

ggplot(plot_df, aes(x = X1)) +
  geom_line(aes(y = Predicted, color = "DIV Estimate"), linewidth = 1) +
  geom_line(aes(y = True, color = "True Partial Mean"), linetype = "dashed", linewidth = 1) +
  scale_color_manual(values = c("DIV Estimate" = "darkgoldenrod1", "True Partial Mean" = "black")) +
  labs(title = "Section 5.3: Partial Effect of X1",
       subtitle = paste("X2 fixed at its median value (", round(x2_fixed, 2), ")"),
       x = "X1", y = "E[Y | do(X1, X2)]",
       color = "Method") +
  theme_minimal()

Conclusion

As discussed in Section 5.3, DIV successfully recovers the partial interventional mean function even when the model must account for the joint distribution of multiple endogenous variables. This capability makes DIV a powerful tool for analyzing complex systems with multiple treatments.