Inferência Causal

Efeito causal médio das unidades sob tratamento

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

Nas últimas aulas, temos focado na estimação de \(\tau = \mathbb{E}[Y(1) - Y(0)]\) em um contexto de dados observacionais.
Hoje, focaremos na estimação de:
- \(\tau_T = \mathbb{E}[Y(1) - Y(0) | Z = 1]\)
- \(\tau_C = \mathbb{E}[Y(1) - Y(0) | Z = 0]\)
Se devemos estimar \(\tau\), \(\tau_T\) ou \(\tau_C\) depende apenas de interesses práticos.

Se \(\tau_T\) e \(\tau_C\) diferirem de \(\tau\), dizemos que o efeitos causais médios são heterogeneos entre os grupos de tratameno e controle

\(\tau_T\)

O efeito causal médio sob as unidades de tratamento é dado por \[\tau_T = \mathbb{E}[Y(1) - Y(0) | Z = 1] = \underbrace{\mathbb{E}[Y(1) | Z = 1]}_{\mathbb{E}[Y | Z = 1]} - \mathbb{E}[Y(0) | Z = 1]\]

O principal problema com estimar \(\tau_T\) é o fato de \(\mathbb{E}[Y(0) | Z = 1]\) ser contrafactual.

Contudo, teremos uma forma de estimar se assumirmos que \[Z \perp\!\!\!\perp Y(0) | X \quad e \quad e(X) < 1 \qquad(1)\]

\(\tau_T\)

Teorema

Sob a suposição Equation 1, temos que \[ \mathbb{E}[Y(0) | Z = 1] = \mathbb{E}[\mathbb{E}(Y | Z = 0, X) | Z = 1] = \int \mathbb{E}[Y| Z = 0, X] f(x |Z = 1)dx \quad \text{ou} = \sum_{k = 1}^K \mathbb{E}[Y| Z = 0, X = k]P(X = k | Z = 1) \]

Isto significa que podemos então estimar \(\tau_T\) por:

\[\hat{\tau}_T = \hat{\bar{Y}}(1) - \sum_{k = 1}^K \hat{\pi}_{[k] | 1} \hat{\bar{Y}}_{[k]}(0),\] em que \(\hat{\pi}_{[k] | 1}= n_{[k]1}/n1\).

\[\hat{\tau}_T = \hat{\bar{Y}}(1) - n_1^{-1}\sum_{i = 1}^N Z_i \hat{\mu}_0(X_i),\] em que \(\hat{\mu}_0(X_i)\) são os valores ajustados da regressão \(\mathbb{E}(Y | Z = 0, X)\) utilizando apenas os dados de controle.

\(\tau_T\)

Exemplo

Seja o modelo linear para todas as unidades dado por \[\mathbb{E}[Y | Z, X] = \beta_0 + \beta_z Z + X \beta_x.\]

Então, \[\begin{align} \tau_T &= \mathbb{E}[Y | Z = 1] - \mathbb{E}[Y(0) | Z = 1] \\ &= \mathbb{E}[Y | Z = 1] - \mathbb{E}[\mathbb{E}(Y | Z = 0, X) | Z = 1] \\ &= \mathbb{E}[\mathbb{E}(Y | Z = 1, X) | Z = 1] - \mathbb{E}[\mathbb{E}(Y | Z = 0, X) | Z = 1] \\ &= \mathbb{E}[\mathbb{E}(Y | Z = 1, X) - \mathbb{E}(Y | Z = 0, X) | Z = 1] \\ &= \beta_z \end{align}\]

Se ajustarmos o modelo de regressão por MQO e obtermos \(\hat{\beta}_0\), \(\hat{\beta}_z\) e \(\hat{\beta}_x\), podemos utilizar \(\hat{\beta}_z\) como estimador de \(\tau_T\).

\(\tau_T\)

Anteriormente vimos que \(\hat{\beta}_z\) era um estimador para \(\tau\).
Agomas vimos que \(\hat{\beta}_z\) era um estimador para \(\tau_T\).
É então \(\hat{\beta}_z\) um estimador para ambos?
Sim!, o modelo linear assume que o efeito causal é constante entre unidades.

\(\tau_T\)

Exemplo

Segundo o Teorema, precisamos apenas de \(\mathbb{E}[Y | Z = 0, X]\) para poder estimar \(\tau_T\). Ou seja, precisamos apenas especificar o modelo para as unidades de controle. Quando o modelo é linear da forma \[\mathbb{E}[Y | Z = 0, X] = \beta_{0|0} + X\beta_{x|0}.\]

Então, \[\begin{align} \tau_T &= \mathbb{E}[Y | Z = 1] - \mathbb{E}[Y(0)| Z = 1] \\ &= \mathbb{E}[Y | Z = 1] - \mathbb{E}[\mathbb{E}[Y | Z = 0, X]| Z = 1] \\ &= \mathbb{E}[Y | Z = 1] - \beta_{0|0} - \mathbb{E}[X| Z = 1]\beta_{x|0} \end{align}\]

Se ajustarmos o modelo de regressão por MQO, temos que: \[\hat{\tau}_T = \hat{\bar{Y}}(1) - \hat{\beta}_{0|0} - \hat{\bar{X}}(1) \hat{\beta}_{x|0}\]

IPS e DR para \(\tau_T\)

Teorema

Sob Equation 1, temos que \[\mathbb{E}[Y(0) | Z = 1] = \mathbb{E} \Big [ \dfrac{e(X)}{e} \dfrac{1-Z}{1 - e(X)}Y \Big] \quad e\]

\[\tau_T = \mathbb{E}[Y | Z = 1] - \mathbb{E} \Big [ \dfrac{e(X)}{e} \dfrac{1-Z}{1 - e(X)}Y \Big],\] em que \(e = P(Z = 1)\) é a probabilidade marginal do tratamento.

Demostração: (no quadro)

IPS e DR para \(\tau_T\)

Então, temos dois estimadores IPW para \(\tau_T\):

\[\hat{\tau}_T^{HT} = \hat{\bar{Y}}(1) - n_1^{-1} \displaystyle \sum_{i = 1}^n \hat{o}(X_i)(1 - Z_i)Y_i \quad e \quad \hat{\tau}_T^{hajek} = \hat{\bar{Y}}(1) - \dfrac{\displaystyle \sum_{i = 1}^n \hat{o}(X_i)(1 - Z_i)Y_i}{\displaystyle \sum_{i = 1}^n \hat{o}(X_i)(1 - Z_i)},\]

em que \(\hat{o}(X_i) = \hat{e}(X_i) / (1 - \hat{e}(X_i))\).

IPS e DR para \(\tau_T\)

Ademais, podemos também obter um estimador DR para \(\mathbb{E}[Y(0) | Z = 1]\).

Definamos \[\tilde{\mu}_{0T} = \mathbb{E}[o(X, \alpha)(1 - Z) \{Y - \mu_0(X, \beta_0)\} + Z \mu_0 (X_1, \beta_0) ]/e,\] em que \(o(X \alpha) = e(X, \alpha) / [1 - e(X, \alpha)].\)

Teorema

Sob Equation 1, se \(e(X, \alpha) = e(X)\) ou \(\mu_0(X, \beta_0) = \mu_0(X)\), então \[\tilde{\mu}_{0T} = \mathbb{E}[Y(0) | Z = 1].\]

IPS e DR para \(\tau_T\)

Estimador DR para \(\tau_T\)

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

Obter os valores ajustados do propensity score (\(e(X_i, \hat{\alpha})\)) e do odds (\(o(X_i, \hat{\alpha}) = e(X_i, \hat{\alpha}) / [1 - e(X_i, \hat{\alpha})]\)).
Obter os valores ajustados do resultado sob controle (\(\mu_0(X_i, \hat{\beta}_0)\)).
Construir o estimador DR para \(\tau_T\) como \(\hat{\tau}_T^{DR} = \hat{\bar{Y}}(1) - \hat{\mu}_{0T}^{DR},\) em que \[\hat{\mu}_{0T}^{DR} = \dfrac{1}{n} \sum_{i = 1}^n [o(X_i, \hat{\alpha}) (1 - Z_i) \{Y_i - \mu_0(X_i, \hat{\beta}_0) \} + Z_i\mu_0(X_i, \hat{\beta}_0)]\]

IPS e DR para \(\tau_T\)

Observação

Note que, podemos re-escrever \(\hat{\tau}_T^{DR}\) como \[\hat{\tau}_T^{reg} - \dfrac{1}{n_1} \sum_{i = 1}^n o(X_i, \hat{\alpha}) (1 - Z_i) \{Y_i - \mu_0(X_i, \hat{\beta}_0)\} \quad \text{ou}\]

\[\hat{\tau}_T^{HT} - \dfrac{1}{n_1}\sum_{i = 1}^n \{o(X_i, \hat{\alpha})(1 - Z_i) + Z_i\} \mu_0(X_i, \hat{\beta}_0).\]

Para estimar a variância de \(\hat{\tau}_T^{DR}\), assim como no caso para \(\hat{\tau}^{DR}\), podemos utilizar Bootstrap.

Exemplo

Implementação

A seguir, vamos implementar 5 estimadores para \(\tau_T\):

2 outcome regression,
2 IPW,
DR.

Implementação

Implementação

Code

att_est <- function(z, y, x, out_family = gaussian, Utruncps = 1) {
  # Outcome regression
  ace_reg_1 <- lm(y ~ z + x)$coef[2]
  outcome_0 <- glm(y ~ x, weights = (1 - z), family = out_family)$fitted.values
  ace_reg_2 <- mean(y[z == 1]) - mean(outcome_0[z == 1])
  
  # IPW
  n <- length(z)
  n_1 <- sum(z)
  ps <- glm(z ~ x, family = binomial)$fitted.values
  ps <- pmin(ps, Utruncps)
  odds <- ps / (1 - ps)
  ace_ipw_1 <- mean(y[z == 1]) - mean(odds * (1 - z) * y) * n / n_1
  ace_ipw_2 <- mean(y[z == 1]) - mean(odds * (1 - z) * y) / mean(odds * (1 - z))
  
  # DR
  res <- y - outcome_0
  ace_dr <- ace_reg_2 - mean(odds * (1 - z) * res) * n / n_1
  
  out <- c(ace_reg_1, ace_reg_2, ace_ipw_1, ace_ipw_2, ace_dr)
  return(out)
}

Implementação

Implementação

Code

os_att <- function(z, y, x, n_boot = 1000, out_family = gaussian, Utruncps = 1) {
  point_estimate <- att_est(z, y, x, out_family, Utruncps)
  
  # NP Bootstrap
  n <- length(z)
  boot_est <- replicate(n_boot, {
    id_boot <- sample(1:n, n, replace = TRUE)
    x <- as.matrix(x)
    att_est(z[id_boot], y[id_boot], x[id_boot, ], out_family, Utruncps)
  })
  boot_sd <- apply(boot_est, 1, sd)
  out <- rbind(point_estimate, boot_sd)
  colnames(out) <- c("Reg1", "Reg2", "HT", "Hajek", "DR")
  rownames(out) <- c("Estim", "SE")
  return(out)
}

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. Em particula, os pesquisadores querem estimar \(\tau_T\) 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)])
resultados_trunc <- os_att(z, y, x, Utruncps = 0.9)
resultados_without_trunc <- os_att(z, y, x)
resultados_trunc

            Reg1       Reg2         HT      Hajek         DR
Estim 0.06124785 -0.3507181 -0.5968341 -0.1922400 -0.2295022
SE    0.22583596  0.2574767  0.5846391  0.3114981  0.2777448

Code

resultados_without_trunc

            Reg1       Reg2         HT      Hajek         DR
Estim 0.06124785 -0.3507181 -1.9920599 -0.3507046 -0.1871062
SE    0.21902572  0.2489611  0.7191328  0.3314690  0.2770278

Como esperado, o estimador HT é sensível ao truncamento. Ademais, Reg1 é bastante diferentes das outras estimativas (Por quê?), já Reg2 é mais parecido com os outros estimadores. Repare que os resultados são diferentes dos resultados obtidos na aula anterir (quando estimamos \(\tau\)) sugerindo a existencia de heterogeneidade entre \(\tau_T\) e \(\tau\).

Outros estimadores

Li et al. (2018) propõem uma classe geral de estimadores para efeitos causais. Considere o seguinte caso geral:

\[\tau^h = \dfrac{\mathbb{E}[h(X) \tau(X)]}{\mathbb{E}[h(X)]},\] em que \(\mathbb{E}(h(X)) \neq 0\).

Assumindo ignorabilidade,

\[\tau^h = \dfrac{\mathbb{E}[h(X) \{ \mu_1(X) - \mu_0(X) \}]}{\mathbb{E}[h(X)]},\]

O que motiva o seguinte estimador,

\[\hat{\tau}^h = \displaystyle \sum_{i = 1}^n h(X_i) [\hat{\mu}_1(X_i) - \hat{\mu}_0(X_i)] \Big / \displaystyle \sum_{i = 1}^n h(X_i).\]

Outros estimadores

Pode-se mostrar que, sob as suposições de ignorabilidade e sobreposição,

\[\tau^h = \mathbb{E} \Big[\dfrac{ZYh(X)}{e(X)} - \dfrac{(1-Z)Yh(X)}{1 - e(X)} \Big] \Big / \mathbb{E}[h(X)],\] o que motiva o estimador IPW para \(\tau^ h\). Note que as unidades de tratamento são ponderadas por \(h(X)/e(X)\) e as unidades de controle por \(h(X)/(1-e(X))\). Li et al. (2018) fornecem um resumo para cada um dos casos particulares de interesse:

	\(h(X)\)	Pesos
\(\tau\)	1	\(1/e(X)\) e \(1/[1 - e(X)]\)
\(\tau_T\)	\(e(X)\)	1 e \(e(X)/[1 - e(X)]\)
\(\tau_C\)	\(1 - e(X)\)	\([1 - e(X)]/e(X)\) e 1
\(\tau_O\)	\(e(X) [1 - e(X)]\)	\(1 - e(X)\) e \(e(X)\)

Referências

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