Inferência Causal

Matching

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

Matching

Considere \(n_0 >> n_1\).
Para cada unidade \(i\) (\(i = 1, \cdots, n_1\)) no grupo de tratamento, buscamos a unidade \(m(i)\) no grupo de controle tal que \(X_i = X_{m(i)}\).
No par que deu macth teremos \(e(X_i) = e(X_{m(i)})\)

Pode-se mostrar que

\[P(Z_i = 1, Z_{m(i)} = 0 | Z_i + Z_{m(i)} = 1, X_i, X_{m(i)}) = 1/2.\]

Condicionado nas covariáveis e no fato de sabermos que em cada par temos um tratamento e um controle, temos um MPE!. Ou seja, podemos analisar match perfeitos em estudos observacionais como se fossem um MPE!

Matching

E se encontrarmos \(M_i\) unidades no grupo de controle que fazem match perfeito com a unidade \(i\) no grupo de tratamento?

Quando os \(M_i\)s variam, o procedimento é chamado variable-ratio matching e pode ser analisado utilizando uma modificação do MPE, a qual apresentamos brevemente a seguir.

Matching

EMPE

Suponha que temos \(n\) matches (\(i = 1, \cdots, n\)).
Para cada \(i\) temos \(1 + M_i\) unidades (1 referente ao tratamento e \(M_i\) referente ao controle).
Assim, o número total de unidades experimentais é \(N = n + \displaystyle \sum_{i + 1}^n M_i.\)
O index \(ij\) representa a unidade \(j\) no conjunto \(i\) (\(i = 1, \cdots, n\), \(j = 1, \cdots, M_i + 1\)) com indicadora de tratamento \(Z_{ij}\) e resultados potenciais \(Y_{ij}(1)\) e \(Y_{ij}(0).\)
Dentro de cada conjunto \(i\), o pesquisador seleciona uma unidade aleatoriamente para atribuir o tratamento.

Matching

EMPE

\[Y_{ij} = Z_{ij}Y_{ij}(1) + (1 - Z_{ij})Y_{ij}(0)\]

O efeito causal médio para o grupo \(i\) é dado por

\[\tau_i = (M_i + 1)^{-1} \displaystyle \sum_{j = 1}^{M_i + 1} Y_{ij}(1) - Y_{ij}(0)\]

Cujo estimador não viesado é dado por

\[\hat{\tau}_i = \displaystyle \sum_{j = 1}^{M_i + 1} Z_{ij}Y_{ij} - M_i^{-1}\sum_{j = 1}^{M_i + 1} (1 - Z_{ij}) Y_{ij}.\]

Matching

EMPE

O efeito causal médio (dentro dos “estratos”) é dado por

\[\tau = n^{-1} \displaystyle \sum_{i = 1}^n \tau_i,\]

Com estimador não viesado dado por

\[\hat{\tau} = n^{-1} \displaystyle \sum_{i = 1}^n \hat{\tau}_i.\]

Observação:

\(\tau\) é o efeito causal médio dentro dos tratamentos, que não é a mesma coisa do que o efeito causal médio das \(N\) unidades.

Matching

EMPE

Efeito causal médio

\[\tau* = N^{-1} \displaystyle \sum_{i = 1}^n \sum_{j = 1}^{M_i + 1}(Y_{ij}(1) - Y_{ij}(0)) = \sum_{i = 1}^n \dfrac{1 + M_i}{N} \tau_i,\]

Em geral, podemos resumir ambos os casos em:

\[\tau_{\omega} = \displaystyle \sum_{i = 1}^n \omega_i \tau_i,\]

para \(\tau\) temos que \(\omega_i = 1/n\) e para \(\tau*\) teremos que \(\omega_i = (1 + M_i) / N.\)

Matching

EMPE

Com estimador não viesado e variância (conservadora):

\[\hat{\tau}_{\omega} = \displaystyle \sum_{i = 1}^n \omega_i \hat{\tau}_i \quad e \quad \hat{V}_{\omega} = \displaystyle \sum_{i = 1}^n c_i(\hat{\tau}_i - \hat{\tau}_{\omega})^2,\] em que \(c_i = \dfrac{\omega_i^2 / (1 - 2\omega_i)}{1 + \sum_{i = 1}^n \omega_i^2 / (1 - 2\omega_i)}\)

Infelizmente não sempre temos match perfeito (\(X_i = X_{m(i)}\)), trazendo consequências negativas ao utilizar MPE ou EMPE em estudos observacionais no contexto de FRT (Guo e Rothenhausler (2023)).

\(X_i \approx X_{m(i)}\)

Mesmo se o grupo de controle for grande, pode ser o caso de não termos match perfeito mas sim \(X_i \approx X_{m(i)}\).
Seja \[m(i) = \arg \min_{k:Z_k = 0} d(X_i, X_k),\] em que \(d(X_i, X_k)\) é uma medida de distância entre \(X_i\) e \(X_k\) (tipicamente distância Euclideana ou de Mahalanobis)

Observação: Daqui em diante, considerarmos matching com reposição. Ou seja, uma mesma unidade pode ser utilizada mais do que uma vez.

Estimação pontual e correção do vies

Para cada \(i\) consideremos \[\hat{Y}_i(1) = Y_i \quad e \quad \hat{Y}_i(0) = M^{-1} \sum_{k \in J_i} Y_k,\] em que \(J_i\) é o conjunto de unidades no grupo de controle que fazen match com a unidade de tratamento \(i\).

Por exemplo, podemos calcular \(d(X_i, X_k)\) \(\forall k\) no grupo de controle e definir \(J_i\) como os índices de \(k\) com os menores \(M\) valores de \(d(X_i, X_k)\).

Estimação pontual e correção do vies

O estimador matching é dado por \[\hat{\tau}^m = n^{-1} \displaystyle \sum_{i = 1}^n (\hat{Y}_i(1) - \hat{Y}_i(0)),\]

Abadei e Imbens mostram que este estimador é viesado e eles propoem o seguinte estimador para o vies

\[\hat{B} = n^{-1} \displaystyle \sum_{i = 1}^n \hat{B}_i,\] em que \(\hat{B}_i = (2Z_i - 1)M^{-1} \displaystyle \sum_{k \in J_i}[\hat{\mu}_{1 -Z_i}(X_i) - \hat{\mu}_{1 -Z_i}(X_k)]\) com \(\hat{\mu_1}(X_i)\) e \(\hat{\mu_0}(X_i)\) sendo os valores ajustados do outcome regression.

Estimação pontual e correção do vies

Assim, o estimador já corrigido pelo vies é dado por

\[\hat{\tau}^{mbc} = \hat{\tau}^m - \hat{B} \equiv n^{-1} \displaystyle \sum_{i = 1}^n\hat{\psi}_i,\] em que \(\hat{\psi}_i = \hat{\mu}_{1}(X_i) - \hat{\mu}_0(X_i) + (2Z_i - 1)(1 + K_i/M)[Y_i - \hat{\mu}_{Z_i}(X_i)]\) com \(K_i\) o número de vezes que a unidade \(i\) é utilizada como match.

Com estimador de variâncias dado por

\[\hat{V}^{mbc} = \dfrac{1}{n^2} \displaystyle \sum_{i = 1}^n (\hat{\psi_i} - \hat{\tau}^{mbc})^2\]

Conexão com o estimador DR

O estimador DR é dado por \[\hat{\tau}^{DR} = \hat{\tau}^{Reg} + \dfrac{1}{n}\displaystyle \sum_{i = 1}^n \Big[ \dfrac{Z_i \hat{R_i}}{\hat{e}(X_i)} - \dfrac{(1 - Z_i) \hat{R_i}}{1 - \hat{e}(X_i)} \Big ],\] em que \(\hat{\tau}^{Reg} = n^{-1} \displaystyle \sum_{i = 1}^n [\hat{\mu}_1(X_i) - \hat{\mu}_0(X_i)]\) e \(\hat{R}_{i} = \left\{ \begin{array}{ll} Y_i - \hat{\mu}_1(X_i) & \text{se } Z_i = 1,\\ Y_i - \hat{\mu}_0(X_i) & \text{se } Z_i = 0. \end{array}\right.\)

Pode-se mostrar que \[\hat{\tau}^{mbc} = \hat{\tau}^{Reg} + \dfrac{1}{n}\displaystyle \sum_{i = 1}^n \Big[ \Big( 1 + \dfrac{K_i}{M} \Big ) Z_i \hat{R_i} - \Big( 1 + \dfrac{M}{K_i} \Big )(1 - Z_i) \hat{R_i}\Big ]\]

Matching para \(\tau_T\)

Suponha agora que estamos interessados em \(\tau_T = \mathbb{E}(Y | Z = 1) - \mathbb{E}(Y(0) | Z = 0)\)
O estimador matching já corrigido pelo vies, é dado por \[\hat{\tau}_T^{mbc} = \hat{\tau}_T^m - \hat{B}_T,\] em que
- \(\hat{\tau}_T^m = n_1^{-1} \displaystyle \sum_{i = 1}^n Z_i[Y_i - \hat{Y}_i(0)]\),
- \(\hat{B}_T = n_1^{-1} \displaystyle \sum_{i = 1}^n Z_i\hat{B}_{T, i}\) e
- \(\hat{B}_{T,i} = M^{-1}\displaystyle \sum_{k \in J_i} [\hat{\mu}_0(X_i) - \hat{\mu}_0 (X_k)]\)

Matching para \(\tau_T\)

Ademais, se fizermos \(\hat{\psi}_{T,i} = Z_i [Y_i - \hat{\mu}_0(X_i)] - (1 - Z_i)K_i/M[Y_i - \hat{\mu}_0(X_i)]\), podemos re-escrever o estimador como \[\hat{\tau}_T^{mbc} = n_1^{-1}\displaystyle \sum_{i = 1}^n \hat{\psi}_{T,i},\] o que motiva o estimador de variância dado por \[\hat{V}_T^{mbc} = \dfrac{1}{n_1^2} \displaystyle \sum_{i = 1}^n (\hat{\psi}_{T,i} - \hat{\tau}_T^{mbc})^2.\]

É possível mostrar que \(\hat{\tau}_T^{mdc} = \hat{\tau}_T^{Reg}-n_1^{-1} \displaystyle \sum_{i = 1}^n \dfrac{K_i}{M}(1 - Z_i)\hat{R}_i.\)

Exemplos

Estudo experimental vs. observacional

Lalonde (1986) está interessado no efeito causal de um programa de treinamento sob o salario. Ele compara os resultados utilizando dados experimentais (dataset lalonde do pacote Matching) vs. dados observacionais (cps1re74.csv).

Estudo Experimental

Code

library(Matching)
library(car)
data(lalonde)
y <- lalonde$re78
z <- lalonde$treat
modelo_ee <- lm(y ~ z)
v_ehw_ee <- hccm(modelo_ee)
Estimate <- coef(modelo_ee)
Std_Error <- sqrt(diag(v_ehw_ee))
t <- Estimate / Std_Error
p_valor <- pnorm(abs(t), lower.tail = FALSE)
tabela_ee <- round(cbind(Estimate,  Std_Error, t, p_valor), 3)
tabela_ee

            Estimate Std_Error      t p_valor
(Intercept) 4554.802   340.749 13.367   0.000
z           1794.343   672.682  2.667   0.004

Estudo Observacional

Code

data <- read.csv("datasets/cps1re74.csv", sep = " ")
y <- data$re78
z <- data$treat
modelo_eo <- lm(y ~ z)
v_ehw_eo <- hccm(modelo_eo)
Estimate <- coef(modelo_eo)
Std_Error <- sqrt(diag(v_ehw_eo))
t <- Estimate / Std_Error
p_valor <- pnorm(abs(t), lower.tail = FALSE)
tabela_eo <- round(cbind(Estimate,  Std_Error, t, p_valor), 3)
tabela_eo

             Estimate Std_Error       t p_valor
(Intercept) 14855.639    76.371 194.519       0
z           -8506.495   584.999 -14.541       0

Reanalisaremos esses dados (mas agora também considerando as covariáveis)

Exemplos

Lalonde: dados experimentais

Code

library(Matching)
library(car)
data(lalonde)
y <- lalonde$re78
z <- lalonde$treat
x <- as.matrix(lalonde[, c("age", "educ", "black", "hisp", "married", "nodegr", "re74", "re75")])

# Experimentos randomizados
neyman_ols <- lm(y ~ z)
fisher_ols <- lm(y ~ z + x)  # ANCOVA
x_c <- scale(x)
lin_ols <- lm(y ~ z*x_c)

out <- c(neyman_ols$coef[2], fisher_ols$coef[2], lin_ols$coef[2],
         sqrt(hccm(neyman_ols, type = "hc2")[2, 2]),
         sqrt(hccm(fisher_ols, type = "hc2")[2, 2]),
         sqrt(hccm(lin_ols, type = "hc2")[2, 2]))

out <- matrix(out, nrow = 3, ncol = 2)
colnames(out) <- c("Est", "SE")
rownames(out) <- c("Neyman", "Fisher", "Lin")
out

            Est       SE
Neyman 1794.343 670.9967
Fisher 1676.343 677.0493
Lin    1621.584 694.7217

Exemplos

Lalonde: dados experimentais

Analisaremos os dados experimentais como se fossem observacionais e para isto, utilizaremos Matching

Code

matching_model <- Match(Y = y, Tr = z, X = x, BiasAdjust = TRUE)
summary(matching_model)


Estimate...  2119.7 
AI SE......  876.42 
T-stat.....  2.4185 
p.val......  0.015583 

Original number of observations..............  445 
Original number of treated obs...............  185 
Matched number of observations...............  185 
Matched number of observations  (unweighted).  268

Ambos, o valor estimado e seu desvio padrão aumentaram. Contudo, qualitativamente as conclusões são as mesmas.

Exemplos

Lalonde: dados observacionais

Code

data <- read.csv("datasets/cps1re74.csv", sep = " ")
y <- data$re78
z <- data$treat
x <- as.matrix(data[, c("age", "educ", "black", "hispan", "married", "nodegree", "re74", "re75")])

# Analisamos os dados como se fossem dados experimentais
neyman_ols <- lm(y ~ z)
fisher_ols <- lm(y ~ z + x)  # ANCOVA
x_c <- scale(x)
lin_ols <- lm(y ~ z*x_c)

out <- c(neyman_ols$coef[2], fisher_ols$coef[2], lin_ols$coef[2],
         sqrt(hccm(neyman_ols, type = "hc2")[2, 2]),
         sqrt(hccm(fisher_ols, type = "hc2")[2, 2]),
         sqrt(hccm(lin_ols, type = "hc2")[2, 2]))

out <- matrix(out, nrow = 3, ncol = 2)
colnames(out) <- c("Est", "SE")
rownames(out) <- c("Neyman", "Fisher", "Lin")
out

             Est        SE
Neyman -8506.495  583.4426
Fisher   704.697  618.2429
Lin    -3689.852 2987.1832

Exemplos

Lalonde: dados observacionais

Agora utilizaremos um método próprio para dados observacionais.

Code

matching_model <- Match(Y = y, Tr = z, X = x, BiasAdjust = TRUE)
summary(matching_model)


Estimate...  2182.4 
AI SE......  869.39 
T-stat.....  2.5103 
p.val......  0.012064 

Original number of observations..............  16177 
Original number of treated obs...............  185 
Matched number of observations...............  185 
Matched number of observations  (unweighted).  227

As conclusões obtidas ao analisarmos os dados observacionais com método apropriados são as mesmas obtids quando analisamos dados experimentais.

Observação

Quando \(X\) tem dimensão grande, matching sobre da maldição da dimensionalidade.
Uma alternativa é realizar o matching baseado no propensity score.
desenvolvem toda a teoria que justifica esta alternativa.

Referências

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