Inferência Causal

Estimador DR

Prof. Carlos Trucíos
ctrucios@unicamp.br

Instituto de Matemática, Estatística e Computação Científica (IMECC),
Universidade Estadual de Campinas (UNICAMP).

Introdução

Introdução

Sob ignorabilidade (\(Z \perp\!\!\!\perp \{ Y(1), Y(0)\} | X\)) e sobreposição (\(0 < e(X) < 1\)), temos mostrado duas fórmulas de identificação do ACE (\(\tau = \mathbb{E}[Y(1) - Y(0)]\)).

Outcome regression

\[\tau = \mathbb{E}[\mu_1(X)]- \mathbb{E}[\mu_0(X)],\]em que \(\mu_i(X) = \mathbb{E}[Y(i) | X] = \mathbb{E}[Y | Z = i, X]\)

IPW

\[\tau = \mathbb{E}\Big[\dfrac{ZY}{e(X)} \Big] - \mathbb{E}\Big[\dfrac{(1- Z)Y}{1 - e(X)} \Big],\]em que \(e(X) = P(Z = 1 | X)\) é o propensity score.

Introdução

  • Outcome regression, precisa ajustar um modelo para o resultado (\(Y\)) dado o tratamento (\(Z\)) e as covariáveis (\(X\)). É consistente se o modelo utilizado está corretamente especificado.
  • IPW precisa ajustar um modelo para o tratamento (\(Z\)) dadas as covariáveis (\(X\)). É consistente se o modelo do propensity score for corretamente especificado.
  • A combinação de ambas as abordagens motiva um novo estimador que também é consistente (conhecido como estimador augmented IPW ou Doubly Robust).

Estimador DR

Estimador DR

  • Postulemos um modelo para as médias condicionais do resultado, \(\mu_1(X, \beta_1)\) e \(\mu_0(X, \beta_0)\).
  • Por exemplo, se a média condicional for linear ou logística, os parâmetros são apenas os coeficientes da regressão.
  • Se o modelo for corretamente especificado, então \(\mu_i(X, \beta_i) = \mu_i(X)\).
  • Postulemos um modelo para o propensity score, \(e(X, \alpha)\).
  • Se o modelo for linear ou logiístico, \(\alpha\) é o coeficiente da regressão e se o modelo for corretamente especificado, \(e(X, \alpha) = e(X)\).

Estimador DR

Definimos, \[\tilde{\mu}_1^{DR} = \mathbb{E} \Big [\dfrac{Z \{Y - \mu_1(X, \beta_1)\}}{e(X, \alpha)} + \mu_1(X, \beta_1) \Big] = \mathbb{E} \Big[\dfrac{ZY}{e(X, \alpha)} - \dfrac{Z - e(X, \alpha)}{e(X, \alpha)}\mu_1(X, \beta_1) \Big],\]

\[\tilde{\mu}_0^{DR} = \mathbb{E} \Big [\dfrac{(1 - Z) \{Y - \mu_0(X, \beta_0)\}}{e(X, \alpha)} + \mu_0(X, \beta_0) \Big] = \mathbb{E} \Big[\dfrac{(1 - Z)Y}{e(X, \alpha)} - \dfrac{e(X, \alpha) - Z}{e(X, \alpha)}\mu_0(X, \beta_0) \Big],\]

  • Note que os termos do meio, aumentam \(\mu_i(X, \beta_i)\) (os estimadores obtidos por outcome regression) pelos residuais ponderados pelo inverse propsensity score.
  • Os termos finais, aumentam o estimador IPW pelos resultados imputados.
Por esse motivo, os estimadores DR são também conhecidos como IPW aumentados.

Estimador DR

Teorema

Assumindo ignorabildiade (\(Z \perp\!\!\!\perp \{ Y(1), Y(0)\} | X\)) e sobreposição (\(0 < e(X) < 1\)):

  1. Se \(e(X, \alpha) = e(X)\) ou \(\mu_1(X, \beta_1) = \mu_1(X)\), então \(\tilde{\mu}_1^{DR} = \mathbb{E}[Y(1)]\).
  2. Se \(e(X, \alpha) = e(X)\) ou \(\mu_0(X, \beta_0) = \mu_1(X)\), então \(\tilde{\mu}_0^{DR} = \mathbb{E}[Y(0)]\).
  3. Se \(e(X, \alpha) = e(X)\) ou \(\{ \mu_0(X, \beta_0), \mu_1(X, \beta_1) = \mu_1(X) \}\), então \(\tilde{\mu}_1^{DR} - \tilde{\mu}_0^{DR} = \tau\)
Pelo Teorem anterior, \[\tilde{\mu}_1^{DR} - \tilde{\mu}_0^{DR} = \tau,\] se ou o modelo de propensity score ou o modelo dos resultados estiver corretamente especificado (motivo pelo qual é chamado de Duplamente Robusto).
Demostração: (no quadro)

Estimador DR

A seguir, apresentamos as versões amostrais de \(\tilde{\mu}_0^{DR}\) e \(\tilde{\mu}_1^{DR}\).

Estimador

Baseados em \(\{X_i, Z_i, Y_i\}_{i = 1}^n\), podemos obter um estimador DR para \(\tau\) através dos seguintes passos:

  1. Obter os valores ajustados do propensity score: \(e(X_i, \hat{\alpha})\).
  2. Obter \(\mu_1(X_i, \hat{\beta}_1)\) e \(\mu_0(X_i, \hat{\beta}_0)\)
  3. Construir \[\hat{\tau}^{DR} = \hat{\mu}_1^{DR} - \hat{\mu}_0^{DR}\], em que \[\hat{\mu}_1^{DR} = \dfrac{1}{n} \displaystyle \sum_{i = 1}^n \Big [\dfrac{Z_i [Y_i - \mu_1(X_i, \hat{\beta}_1)]}{e(X_i, \hat{\beta}_1)} + \mu_1(X_i, \hat{\beta}_1)\Big] \quad e\]

\[\hat{\mu}_0^{DR} = \dfrac{1}{n} \displaystyle \sum_{i = 1}^n \Big [\dfrac{(1 - Z_i) [Y_i - \mu_0(X_i, \hat{\beta}_0)]}{1 - e(X_i, \hat{\beta}_0)} + \mu_0(X_i, \hat{\beta}_0)\Big]\]

Funk et al. (2011) sugerem aproximar \(\mathbb{V}(\hat{\tau}^{DR})\) vía Bootstrap não paramétrico.

Estimador DR

Implementação

Implementaremos os quatro métodos vistos até agora. O modelo utilizado para o propsensity score será o padrão (modelo logístico).

Code 
OS_est <- function(z, y, x, out_family = gaussian, alphas = c(0, 1)) {
  # Propensity score
  e <- glm(z ~ x, family = binomial)$fitted.values
  e <- pmax(alphas[1], pmin(alphas[2], e))
  
  # Fitted potencial outcomes
  outcome_1 <- glm(y ~ x, weights = z, family = out_family)$fitted.values
  outcome_0 <- glm(y ~ x, weights = 1 - z, family = out_family)$fitted.values
  
  # Outcome regression estimator
  tau_reg <- mean(outcome_1 - outcome_0)
  
  # IPW
  tau_ipw <- mean(z * y / e) - mean((1 - z) * y / (1 - e))
  tau_hajek <- mean(z * y / e) / mean(z / e) - mean((1 - z) * y / (1 - e)) / mean((1 - z) / (1 - e))
  
  # DR
  res_1 <- y - outcome_1
  res_2 <- y - outcome_0
  mu_1_dr <- mean(z * res_1 / e + outcome_1) 
  mu_0_dr <- mean((1 - z) * res_2 / (1 - e) + outcome_0)
  tau_dr <- mu_1_dr - mu_0_dr
  
  out <- c(tau_reg, tau_ipw, tau_hajek, tau_dr)
  return(out)
}

Augmented IPW

Implementação

Utilizaremos Bootstrao não paramétrico para estimar a variância dos estimadores

Code 
OS_ATE <- function(z, y, x, n_boot = 1000, out_family = gaussian, alphas = c(0, 1)) {
  taus <- OS_est(z, y, x, out_family, alphas)
  
  # Bootstrap não-paramétrico
  n <- length(z)
  x <- as.matrix(x)
  boot_est <- replicate(n_boot, {
    id_boot <- sample(1:n, n, replace = TRUE)
    OS_est(z[id_boot], y[id_boot], x[id_boot, ], out_family, alphas)
  })
  boot_se <- apply(boot_est, 1, sd)
  out <- rbind(taus, boot_se)
  rownames(out) <- c("Estim", "SE")
  colnames(out) <- c("Reg", "IPW", "Hajek", "DR")
  return(out)
}

Simulação

Simulação

Avaliaremos as propriedades em amostras finitas dos quatro estimadores implementados. O estudo considerará o seguintes cenârios:

  • Ambos, o modelo do propensity score e do resultados são corretos.
  • O modelo do propensity score é incorreto mas o do resultado correto.
  • O modelo do propensity score é correto mas o do resultado incorreto.
  • Ambos os modelos são incorretos

Simulação

Simulação

Code 
mc <- 500
n <- 500
bias <- matrix(NA, ncol = 4, nrow = mc)
estim_se <- matrix(NA, ncol = 4, nrow = mc)
for (i in 1:mc) {
  x <- matrix(rnorm(n * 2), n, 2)
  x1 <- cbind(1, x)
  beta_z <- c(0, 1, 1)
  e <- 1 / (1 + exp(-as.vector(x1 %*% beta_z)))
  z <- rbinom(n, 1, e)
  beta_y1 <- c(1, 2, 1)
  beta_y0 <- c(1, 2, 1)
  y1 <- rnorm(n, x1 %*% beta_y1)
  y0 <- rnorm(n, x1 %*% beta_y0)
  y <- z * y1 + (1 - z) * y0 
  ce <- OS_ATE(z, y, x)
  bias[i,] <-  ce[1, ] - 0
  estim_se[i, ] <- ce[2, ]
}
# Bias, True SD, Estim SD
resultados <- rbind(apply(bias, 2, mean), 
                    apply(bias, 2, sd), 
                    apply(estim_se, 2, mean))
colnames(resultados) <- c("Reg", "IPW", "Hajek", "DR")
rownames(resultados) <- c("Bias", "True Var", "Estim Var")
resultados
                  Reg         IPW      Hajek          DR
Bias      0.005645807 0.004843971 0.01945095 0.005782348
True Var  0.107839132 0.313688570 0.28074983 0.123059224
Estim Var 0.104973769 0.256120961 0.22858887 0.120900388

Simulação

Simulação

Code 
mc <- 500
n <- 500
bias <- matrix(NA, ncol = 4, nrow = mc)
estim_se <- matrix(NA, ncol = 4, nrow = mc)
for (i in 1:mc) {
  x <- matrix(rnorm(n*2), n, 2)
  x1 <- cbind(1, x, exp(x))
  beta_z <- c(-1, 0, 0, 1, -1)
  e <- 1/(1 + exp(-as.vector(x1 %*% beta_z)))
  z <- rbinom(n, 1, e)
  beta_y1 <- c(1, 2, 1, 0, 0)
  beta_y0 <- c(1, 2, 1, 0, 0)
  y1 <- rnorm(n, x1 %*% beta_y1)
  y0 <- rnorm(n, x1 %*% beta_y0)
  y <- z*y1 + (1 - z)*y0
  ce <- OS_ATE(z, y, x)
  bias[i,] <-  ce[1, ] - 0
  estim_se[i, ] <- ce[2, ]
}
# Bias, True SD, Estim SD
resultados <- rbind(apply(bias, 2, mean), 
                    apply(bias, 2, sd), 
                    apply(estim_se, 2, mean))
colnames(resultados) <- c("Reg", "IPW", "Hajek", "DR")
rownames(resultados) <- c("Bias", "True Var", "Estim Var")
resultados
                  Reg        IPW      Hajek          DR
Bias      0.003172167 -0.7754680 -0.6661254 -0.01207538
True Var  0.111727438  1.9938424  0.5684628  0.25101684
Estim Var 0.118663850  0.6675723  0.3870083  0.19131786

Simulação

Simulação

Code 
mc <- 500
n <- 500
bias <- matrix(NA, ncol = 4, nrow = mc)
estim_se <- matrix(NA, ncol = 4, nrow = mc)
for (i in 1:mc) {
  x <- matrix(rnorm(n*2), n, 2)
  x1 <- cbind(1, x, exp(x))
  beta_z <- c(0, 1, 1, 0, 0)
  e <- 1 / (1 + exp(-as.vector(x1 %*% beta_z)))
  z <- rbinom(n, 1, e)
  beta_y1 <- c(1, 0, 0, 0.2, -0.1)
  beta_y0 <- c(1, 0, 0, -0.2, 0.1)
  y1 <- rnorm(n, x1 %*% beta_y1)
  y0 <- rnorm(n, x1 %*% beta_y0)
  y <- z*y1 + (1 - z)*y0
  ce <- OS_ATE(z, y, x)
  bias[i,] <-  ce[1, ]  - 0.2*exp(1/2)
  estim_se[i, ] <- ce[2, ]
}
# Bias, True SD, Estim SD
resultados <- rbind(apply(bias, 2, mean), 
                    apply(bias, 2, sd), 
                    apply(estim_se, 2, mean))
colnames(resultados) <- c("Reg", "IPW", "Hajek", "DR")
rownames(resultados) <- c("Bias", "True Var", "Estim Var")
resultados
                  Reg          IPW       Hajek          DR
Bias      -0.04640363 0.0008687522 0.002644223 0.001847562
True Var   0.11470270 0.1729404572 0.167957141 0.164342718
Estim Var  0.11303603 0.1511667243 0.139691558 0.140591620

Simulação

Simulação

Code 
mc <- 500
n <- 500
bias <- matrix(NA, ncol = 4, nrow = mc)
estim_se <- matrix(NA, ncol = 4, nrow = mc)
for (i in 1:mc) {
  x <- matrix(rnorm(n*2), n, 2)
  x1 <- cbind(1, x, exp(x))
  beta_z <- c(-1, 0, 0, 1, -1)
  e <- 1 / (1 + exp(-as.vector(x1 %*% beta_z)))
  z <- rbinom(n, 1, e)
  beta_y1 <- c(1, 0, 0, 0.2, -0.1)
  beta_y0 <- c(1, 0, 0, -0.2, 0.1)
  y1 <- rnorm(n, x1 %*% beta_y1)
  y0 <- rnorm(n, x1 %*% beta_y0)
  y <- z*y1 + (1 - z)*y0
  ce <- OS_ATE(z, y, x)
  bias[i,] <-  ce[1, ] - 0.2*exp(1/2)
  estim_se[i, ] <- ce[2, ]
}
# Bias, True SD, Estim SD
resultados <- rbind(apply(bias, 2, mean), 
                    apply(bias, 2, sd), 
                    apply(estim_se, 2, mean))
colnames(resultados) <- c("Reg", "IPW", "Hajek", "DR")
rownames(resultados) <- c("Bias", "True Var", "Estim Var")
resultados
                  Reg       IPW       Hajek        DR
Bias      -0.06756789 0.1077739 -0.05767143 0.1441579
True Var   0.12494688 0.2606137  0.18246527 0.2630190
Estim Var  0.12725755 0.2407535  0.16386663 0.2357870

Exemplo

Exemplo

Exemplo

Dataset disponível aqui. Chang et al. (2016) estudam se a participação no programa da merenda nas escolas (School_meal) leva a um aumento no índice de massa muscular (BMI) das crianças. 12 outras covariáveis também estão disponíveis no dataset.

Code 
dados <- read.csv("datasets/nhanes_bmi.csv")[, -1]
z <- dados$School_meal
y <- dados$BMI
x <- scale(dados[, -c(1, 2)])
causal_effetcs = OS_ATE(z, y, x, n_boot = 1000)
round(causal_effetcs, 3)
         Reg    IPW  Hajek     DR
Estim -0.017 -1.516 -0.156 -0.019
SE     0.220  0.484  0.242  0.223
         Reg    IPW  Hajek     DR
Estim -0.017 -0.713 -0.054 -0.043
SE     0.233  0.410  0.244  0.239

Comentários Finais

Comentários Finais

  • Pode-se mostrar que:
    • \(\hat{\tau}^{DR} = \hat{\tau}^{Reg} + \dfrac{1}{n} \displaystyle \sum_{i = 1}^n \dfrac{Z_i (Y_i - \mu_1(X_i, \hat{\beta}_1))}{e(X_i, \hat{\alpha})} - \dfrac{1}{n} \sum_{i = 1}^n \dfrac{(1 - Z_i)(Y_i - \mu_0(X_i, \hat{\beta}_0))}{1 - e(X_i, \hat{\alpha})}.\)
    • \(\hat{\tau}^{DR} = \hat{\tau}^{IPW} - \dfrac{1}{n} \displaystyle \sum_{i = 1}^n \dfrac{Z_i - e(X, \hat{\alpha})}{e(X_i, \hat{\alpha})} \mu_1(X_i, \hat{\beta}_1) + \dfrac{1}{n} \displaystyle \sum_{i = 1}^n \dfrac{e(X_i, \hat{\alpha})-Z_i}{1 - e(X_i, \hat{\alpha})} \mu_0(X_i, \hat{\beta}_0)\)
  • O estimador RD foi proposto por Scharfstein et al. (1999) mas recentemente ressucitado por Chernozhukov et al. (2018) com o nome de double machine learning. A ideia básica é substituir o modleo utilizado no propensity score e no modelo para o resultado, por modelo de aprendizado de máquina.

Referências

  • Peng Ding (2023). A First Course in Causal Inference. Capítulo 12.

Carlos Trucíos (IMECC/UNICAMP) | ME920/MI628 - Inferência Causal | ctruciosm.github.io