Aula 8 - Inferência Causal 2
IESP-UERJ
Aleatorização
Como vimos na aula passada o problema da inferência causal é a impossibilidade de obeservarmos todos reultados potenciais de um tratamento. Por resultados potenciais queremos dizer situações contafatuais onde observamos o valor de uma variável resposta \(Y\) de uma observação \(i\) que foi exposta a um ‘tratamento’ \(T\) e o valor desta variável caso esta mesma observação não tivesse recebido esse tratamento. O efeito causal é a diferença entre os resultados potenciais. Como o efeito causal pode variar entre pessoas, estamos interessados no efeito causal médio (em inglês, Average Treatment Effect ou ATE).
Mostrei que, ao estimar o efeito causal médio, estamos sujeitos a um viés que ocorre toda vez que os resultados potenciais (\(Y_{0i}\) e \(Y_{1i}\)) dependerem do tratamento \(T\). Essa dependência ocorre quando algum fator faz com que certos tipos de pessoas tenham maior ou menor probabilidade de receber o tratamento ou o controle, daí a alcunha "viés de seleção’. Se pudermos controlar estes fatores o viés de seleção desaparece. Ilustrei isso com algumas simulações. Fora do mundo simulado é difícil controlar todas as variáveis de confusão. Quando omitimos essas variáveis em noso modelo geramos o viés de seleção, ou seja temos o que se chama de Viés de Variável Omitida. As equações abaixo ajudam a ilustrar o porque disso.
Vamos supor que temos um efeito direto de BF em AP e a correta especificação de nosso modelo seja
\[ AP_{i} = \beta_{0} + \beta_{1} BF_{i} + \beta_{2} x_{i} + \epsilon_{i} \]
Se não observamos \(x_{i}\) estaremos ajustando o seguinte modelo
\[ AP_{i} = \beta_{0}^{*} + \beta_{1}^{*} BF_{i} + \epsilon_{i}^{*} \]
Se BF tem relação com \(x_{i}\) teríamos algo como
\[ x_{i} = \gamma_{0} + \gamma_{1} BF_{i} + \nu_{i} \]
O que substituindo acima nos dá
\[ AP_{i} = \beta_{0} + \beta_{2} \gamma_{0} + (\beta_{1} + \beta_{2}\gamma_{1}) BF_{i} + \epsilon_{i} + \nu_{i} \]
Então quando estimamos o coeficiente sem \(x_{i}\), ou seja, com essa variável omitida temos
\[\beta_{1}^{*} = \beta_{1} + \beta_{2}^{*}\gamma_{1}\]
Como fazemos para nos livrar desse viés se não observamos a variável omitida? Se atribuimos o tratamento de maneira aleatória tornamos esse tratamento independente do resultado potencial, obtendo
\[ \begin{align} \text{E}[Y_{i}\mid T_{i} = 1] - \text{E}[Y_{i}\mid T_{i} = 0] &= \text{E}[Y_{1i}\mid T_{i} = 1] - \text{E}[Y_{0i}\mid T_{i} = 0]\\ &= \text{E}[Y_{1i}\mid T_{i} = 1] + \text{E}[Y_{0i}\mid T_{i} = 1] \end{align} \]
Como \(Y_{i} \perp T_{i} \mid X\) ao controlar pelas variáveis de confusão \(X\) os resultados potencias de quem recebeu o tratamento se distribui da mesma maneira aos de quem não recebeu. Por isso as duas linhas da equação são equivalentes. Isso acontece por que ao atribuir o tratamento de forma aleatória tendemos a anular diferenças sistemáticas que poderiam causar algum viés.
Se atribuírmos o tratamento de forma aleatória podemos identificá-lo usando nossa estratégia de modelagem:
\[ Y_{i} \sim f(y_{i}|\theta_{i},\alpha) \\ \theta_{i} = g(\alpha + \rho T) \]
Onde \(\alpha = Y_{0i}\) e \(\rho\) é o efeito do tatamento \((Y_{1i} - Y_{0i})\).
Como podemos aleatorizar o tratamento? Podemos fazer um experimento onde certas pessoas recebem o tratamento e outras não recebem por sorteio. Um exemplo clássico é o experimento STAR onde procurou-se medir o efeito do tamanho da sala de aula na aprendizagem. Esse experimento foi feito no Tennessee (EUA) e consistiu em atribuir 11.600 alunos do infantil a classes de diferentes tamanhos por meio de sorteio. Esses alunos foram acompanhados até o terceiro ano e comparou-se o desemepnho desses alunos em um teste (ver Angrist e Pischke, p. 17). A pesquisa constatou que alunos em classes menores obtiveram um rendimento 5% maior no teste.
Outro exemplo é o experimento feito por King et al (2007) para ver os efeitos do programa de seguridade social no México. Esse programa de universalização de saúde envolvia a construção de hospitais e o acesso à medicamentos e outras intervenções médicas. O país foi dividido em mais de 7 mil clusters de saúde e King e sua equipe escolheram 148 destes clusters para participar do experimento sorteando 74 para receber o tratamento e 74 para controle. Esses 148 clusters envolveram mais de 500 mil pessoas se tornando o maior experimento aleatorizado em saúde pública.
A aleatorização garante que características individuais que poderiam causar viés estejam balanceadas nos grupos, o que pode ser constatado comparando-se as médias dessas características. Além da aleatorização em intervenções como a do STAR e do SPS podemos nos valer de uma amostra aleatória da população e atribuir parte desta amostra ao grupo de tratamento e outra parte ao grupo de controle. Isso é o que fazemos quando usamos experimentos de survey, por exemplo.
Ainda outra forma de aleatorização se dá quando temos os chamados experimentos naturais. Zimmerman (2003) se aproveitou do fato de que em determinada faculdade os alunos residentes no campus eram atribuídos de forma aleatória, com relação ao desempenho, a quartos compartilhados, para medir o efeito dos pares (companheiros de quarto) neste desempenho, constatando que alunos no meio da distribuição do teste desempenhavam pior quando tinham companheiors de quarto com baixo desempenho.
Outro exemplo é o do mexicano Progressa (programa de trasferência de renda condicional). As comunidades que recebiam o programa eram escolhidas pelo nível de pobreza. Depois das comunidades consideradas muito pobres e pobres as comunidades de renda média foram alocadas por quintis de renda, dos piores aos melhores. Isso permitiu o uso de RDD, desenho experimental que vermeos mais à frente.
Em todos esses exemplos, além do efeito do tratamento também foram incluídas variáveis de pré-tratamento. Essas variáveis não são variáveis de confusão, pois, como a atribuição do tratamento é feita de modo aleatório, não há necessidade de controlar por elas. Essas variáveis de pré tratamento, como sexo e cor, por exemplo, assim como sua interação com o tratmento, ajudam a melhorar a eficiência dos estimadores, desde que não sejam afetadas pelo tratamento.
Dados observacionais
Quando não podemos recorrer a experimentos podemos nos valer de alguma forma de ajuste estatístico para nos aproximar da aleatorização. Esses ajustes procuram garantir a premissa da independência condicional dos resultados potenciais com relação ao tratamento que permite a identificação do efeito causal. Essa premissa, também chamada de ignorabilidade ou seleção em observáveis, garante que, se controlamos por todas as variáveis que possam causar viés, podemos tratar nossas observações como se fossem resultado de um experimento.
Para garantir a ignorabilidade nosso objetivo é comparar unidades o mais parecidas o possível. Não é necessário que a probabilidade dos dois tratamentos seja a mesma, mas a probabilidade de tratamento, condicional às variáveis de confusão, Pr(\(T = 1 \mid X\)) deve ser a mesma para todos as observações. Se garantimos a ignorabilidade podemos expressar o efeito causal como
\[ \delta_{TOT} = \text{E}(\delta_X \mid T_{i} = 1)\]
onde
\[ \delta_{X} \equiv \text{E}(y_{i}\mid X{i} = x, T_{i} = 1) - \text{E}(y_{i}\mid X{i} = x, T_{i} = 0)\]
Isso é, se assumirmos a ignorabilidade, a diferença entre tratados e não tratados, condicional aos possíveis valores que variáveis de controle assumem quando \(T = 1\), nos dá o efeito causal do tratamento nos tratados (\(\delta_{TOT}\)). Se \(X\) assume valores fora daqueles abrangidos pelo tratamento temos que extrapolar de nossos dados e estimamos o efeito médio do tratamento (\(\delta_{ATE}\)). Usando o exemplo do BF podemos associar o \(\delta_{TOT}\) ao efeito do tratamento comparando pessoas elegíveis para participar no programa que receberam ou não o BF. \(\delta_{ATE}\) nos dá o efeito entre os que receberam o tratamento (BF) e todos os demais.
O balancemento da distribuição dos valores dos controles é feito por matching, isto é, compara-se observações de controle que possuem distribuições de covariantes muito parecidas com as das observações que recebem o tratamento. Quando isso não ocorre, isto é, quando não há um balanceamento no tratamento ou quando não há sobreposição nas covariantes dos dois grupos exige-se muito do modelo forçando ele a extrapolar para estas regiões. Quando usamos regressão linear, por exemplo, e as observações não são balanceadas pode surgir um viés em regiões fora da região do TOT e o efeito causal pode ser subestimado ou superestimado.
Na prática o matching é feito aplicando-se algum algoritmo que calacula a ‘distância’ entre observações onde as dimensões são as covariantes. Entre os métodos mais utilizados estão o Mahalanobis Distance Matching (MDM), que escolhe o par conforme a distância seja menor que um certo valor \(\sqrt{(X_{i} - X_{j})S^{-1}(X_{i} - X_{j})} < \delta\), onde \(S\) é a matriz de covariânca de \(X\), e o Coarsened Exact Matching (CEM) onde se compara subespaços de X conforme se agrupa as categorias ou valores de uma ou mais variáveis (juntando quem tem ensino médio e superior completo, por exemplo) e o matching se dá quando \(C_{\delta}X_{i} = C_{\delta}X_{j}\), sendo \(\delta\) aqui o grau de granularidade. Esses algoritmos são implementados em pacotes como o MatchIt que também usa como método o matching exato ou o KNN entre outros.
library(MatchIt)
data("lalonde")
m.out2 <- matchit(treat ~ age + educ + race + nodegree +
married + re74 + re75, data = lalonde,
distance = "mahalanobis", replace = TRUE,
exact = ~ married + race)
m.out2## A matchit object
## - method: 1:1 nearest neighbor matching with replacement
## - distance: Mahalanobis
## - number of obs.: 614 (original), 263 (matched)
## - target estimand: ATT
## - covariates: age, educ, race, nodegree, married, re74, re75summary(m.out2)##
## Call:
## matchit(formula = treat ~ age + educ + race + nodegree + married +
## re74 + re75, data = lalonde, distance = "mahalanobis", exact = ~married +
## race, replace = TRUE)
##
## Summary of Balance for All Data:
## Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
## age 25.8162 28.0303 -0.3094 0.4400 0.0813
## educ 10.3459 10.2354 0.0550 0.4959 0.0347
## raceblack 0.8432 0.2028 1.7615 . 0.6404
## racehispan 0.0595 0.1422 -0.3498 . 0.0827
## racewhite 0.0973 0.6550 -1.8819 . 0.5577
## nodegree 0.7081 0.5967 0.2450 . 0.1114
## married 0.1892 0.5128 -0.8263 . 0.3236
## re74 2095.5737 5619.2365 -0.7211 0.5181 0.2248
## re75 1532.0553 2466.4844 -0.2903 0.9563 0.1342
## eCDF Max
## age 0.1577
## educ 0.1114
## raceblack 0.6404
## racehispan 0.0827
## racewhite 0.5577
## nodegree 0.1114
## married 0.3236
## re74 0.4470
## re75 0.2876
##
##
## Summary of Balance for Matched Data:
## Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
## age 25.8162 24.8973 0.1284 0.8265 0.0330
## educ 10.3459 10.5784 -0.1156 0.9046 0.0179
## raceblack 0.8432 0.8432 0.0000 . 0.0000
## racehispan 0.0595 0.0595 -0.0000 . 0.0000
## racewhite 0.0973 0.0973 0.0000 . 0.0000
## nodegree 0.7081 0.5892 0.2616 . 0.1189
## married 0.1892 0.1892 0.0000 . 0.0000
## re74 2095.5737 1780.4016 0.0645 1.8106 0.0429
## re75 1532.0553 1161.2419 0.1152 2.0101 0.0258
## eCDF Max Std. Pair Dist.
## age 0.1676 0.3400
## educ 0.1189 0.3038
## raceblack 0.0000 0.0000
## racehispan 0.0000 0.0000
## racewhite 0.0000 0.0000
## nodegree 0.1189 0.3329
## married 0.0000 0.0000
## re74 0.2000 0.2483
## re75 0.0811 0.2086
##
## Percent Balance Improvement:
## Std. Mean Diff. Var. Ratio eCDF Mean eCDF Max
## age 58.5 76.8 59.5 -6.2
## educ -110.3 85.7 48.4 -6.8
## raceblack 100.0 . 100.0 100.0
## racehispan 100.0 . 100.0 100.0
## racewhite 100.0 . 100.0 100.0
## nodegree -6.8 . -6.8 -6.8
## married 100.0 . 100.0 100.0
## re74 91.1 9.7 80.9 55.3
## re75 60.3 -1462.3 80.8 71.8
##
## Sample Sizes:
## Control Treated
## All 429. 185
## Matched (ESS) 34.96 185
## Matched 78. 185
## Unmatched 351. 0
## Discarded 0. 0plot(summary(m.out2))
m.data1 <- match.data(m.out2)
head(m.data1)## treat age educ race married nodegree re74 re75 re78 weights
## NSW1 1 37 11 black 1 1 0 0 9930.0460 1
## NSW2 1 22 9 hispan 0 1 0 0 3595.8940 1
## NSW3 1 30 12 black 0 0 0 0 24909.4500 1
## NSW4 1 27 11 black 0 1 0 0 7506.1460 1
## NSW5 1 33 8 black 0 1 0 0 289.7899 1
## NSW6 1 22 9 black 0 1 0 0 4056.4940 1fit1 <- lm(re78 ~ treat + age + educ + race + married + nodegree + re74 + re75, data = m.data1, weights = weights)
display(fit1)## lm(formula = re78 ~ treat + age + educ + race + married + nodegree +
## re74 + re75, data = m.data1, weights = weights)
## coef.est coef.se
## (Intercept) -3191.25 4446.48
## treat 697.42 995.40
## age 44.55 64.63
## educ 625.45 306.99
## racehispan 1547.77 1934.02
## racewhite 1538.40 1538.74
## married 390.75 1236.24
## nodegree 808.93 1379.16
## re74 0.03 0.13
## re75 0.18 0.21
## ---
## n = 263, k = 10
## residual sd = 7283.91, R-Squared = 0.04Controle em variáveis pós-tratamento e efeitos de mediação
Não se deve controlar por variáveis que são afetadas pelo tratamento. Se estamos interessado no efeito da escolaridade na renda, para pegar um exemplo de Angrit e Pischke, e se a escolaridade também afeta a ocupação futura da pessoa, não devemos controlar por esta variável, pois isso faz com que estejamos comparando coisas diferentes. Ao controlar pela ocupação estamos comparando tratados e não tratados com a mesma ocupação. No entanto, como a ocupação depende do tratamento, os não tratados tem na base ocupação diferente dos tratados. O efeito causal está, portanto, desbalanceado.
Outro problema ocorre quando controlamos por efeitos intermediários, ou por mediações. Suponha, por exemplo, que o efeito de BF sobre a aprovação presidencial passa pelo seu efeito na ideologia, medida pela identificação partidária. Um pesquisador pode querer saber o quanto o efeito do BF na verdade se deve à ideologia. Ele pode achar que se o efeito do BF cai em 50% ao controlar por ideologia pode significar que metade do efeito do BF em AP se deve a esta variável. Isso não é necessariamente verdade, pois, como ocorre ao controlar por variáveis pós-tratamento, a ideologia pré-tratamento nos grupos controle e tratamento seria, necessariamente, diferente.
Propensity Score
O teorema do propensity score nos diz que se os resultados potenciais são independentes do tratamento condicional ao vetor de covariantes \(X\) eles também serão independentes a um escalar resultante de uma função destas covariantes, o propensity score \(P(X) \equiv \text{E}(T_{i} = 1 \mid X) = P(T_{i} = 1 \mid X)\).
Esse teorema nos diz que, assim como com o viés de variável omitida, se controlarmos por covariantes que afetam a probabilidade de tratamento, nos livramos do viés de seleção. O teorema diz mais, ele diz que só precisamos controlar por esta probabilidade, isto é, fazemos o matching de unidades tratadas com unidades com mesmo propensity score.
A questão é como modelar e estimar \(P(X)\), como parameterizar as covariantes, ou o quão flexível o modelo deve ser. Uma possibilidade, apontada por Angrist e Piscke, é utilizar um logit com polinômios das covariantes. Outra possibilidade é utilizar métodos não-paramétricos.
Um dos problemas com o propensity score é que, assintóticamente, ele não é tão eficiente quanto uma regressão em todas as covariantes. Por outro lado ele pode apresentar ganhos de eficiência em amostras finitas, ainda mais se aproximado por métodos não paramétricos. O propensity score, ainda segundo Angrist e Pischke, nos força a dar maior atenção à modelagem. No final do cap. 3 esses autores comparam o propensity score com a regressão nas covariantes e mostram a vantagem desta última.
Quando se usa alguma técnica de matching, em geral, procura-se pares para as unidades tratadas formando “duplas”. Se temos mais observações que não receberam tratamento do que unidades tratadas ao estabelecer as duplas descartamos observações que não deram matching e reduzimos nossos dados apenas à faixa de valores abrangidas pelo tratamento procurando garantir sobreposição, além do balanceamento.
#
# Exemplo de Propensity Score Gelman e Hill (2007) pp. 207 - 209
#
cc2 <- read.csv("cc2.csv")
library ("arm")
## Calculando o propensity score matches
# Primeiro prevemos quem recebe o tratamento com base nos preditores que consideramos necessários
# para garantir a ignorabilidade.
# A ideia é ir testando modelos para ver qual o melhor match.
# Primeiro modelo
ps.fit.1 <- glm (treat ~ as.factor(educ) + as.factor(ethnic) + b.marr +
work.dur + prenatal + momage + sex + first + preterm + age +
dayskidh + bw + unemp.rt, data=cc2, family=binomial(link="logit"))
# (...)
# modelo final
ps.fit.2 <- glm (treat ~ bwg + as.factor(educ) + bwg:as.factor(educ) +
as.factor(ethnic) + b.marr + as.factor(ethnic):b.marr + work.dur +
prenatal + preterm + age + momage + sex + first, data=cc2,
family=binomial(link="logit"))
# Desse modelo tiramos os valores preditos que serão usados para montar o matching
pscores <- predict(ps.fit.2, type = "link")
# A função metching encontra a unidade controle com o score mais próximo da de tratamento e forma pares.
matches <- matching (z=cc2$treat, score=pscores)
matched <- cc2[matches$match.ind,]
# Podemos verificar o balanceamento pós matching
# Só visual
summary(matched[matched$treat == 1,c("educ", "ethnic","age","momage", "sex","bwg")])## educ ethnic age momage
## Min. :1.000 Min. :1.00 Min. :52.35 Min. :13.00
## 1st Qu.:1.000 1st Qu.:2.00 1st Qu.:55.25 1st Qu.:20.00
## Median :2.000 Median :2.00 Median :56.64 Median :24.00
## Mean :1.966 Mean :2.31 Mean :56.67 Mean :24.44
## 3rd Qu.:3.000 3rd Qu.:3.00 3rd Qu.:58.22 3rd Qu.:28.00
## Max. :4.000 Max. :3.00 Max. :60.86 Max. :41.00
## sex bwg
## Min. :0.0000 Min. :0.0000
## 1st Qu.:0.0000 1st Qu.:0.0000
## Median :1.0000 Median :0.0000
## Mean :0.5069 Mean :0.4897
## 3rd Qu.:1.0000 3rd Qu.:1.0000
## Max. :1.0000 Max. :1.0000summary(matched[matched$treat == 0,c("educ", "ethnic","age","momage", "sex","bwg")])## educ ethnic age momage
## Min. :1.00 Min. :1.000 Min. : 5.00 Min. :17.00
## 1st Qu.:1.00 1st Qu.:2.000 1st Qu.: 30.25 1st Qu.:21.00
## Median :2.00 Median :2.000 Median : 59.00 Median :24.00
## Mean :1.81 Mean :2.359 Mean : 56.09 Mean :23.69
## 3rd Qu.:2.00 3rd Qu.:3.000 3rd Qu.: 84.00 3rd Qu.:26.00
## Max. :4.00 Max. :3.000 Max. :104.00 Max. :31.00
## sex bwg
## Min. :0.0000 Min. :0.0000
## 1st Qu.:0.0000 1st Qu.:1.0000
## Median :1.0000 Median :1.0000
## Mean :0.5069 Mean :0.7897
## 3rd Qu.:1.0000 3rd Qu.:1.0000
## Max. :1.0000 Max. :1.0000covs <- c('bw', 'preterm', 'dayskidh', 'sex', 'first', 'age', 'black', 'hispanic', 'white', 'b.marr', 'lths', 'hs', 'ltcoll', 'college', 'work.dur', 'prenatal', 'momage')
cov_names <- c('birth weight', 'weeks preterm', 'days in hospital', 'male', 'first born', 'age', 'black', 'hispanic', 'white', 'unmarried at birth', 'less than high school', 'high school graduate', 'some college', 'college graduate', 'worked during pregnancy', 'had no prenatal care', 'age at birth')
bal_nr <- balance(rawdata=cc2[,covs], treat=cc2$treat, matched=matches$cnts, estimand='ATT')## Balance diagnostics assume that the estimand is the ATTbal_nr## Balance Statistics for Unmatched Data
## --
## treat control unstd.diff abs.std.diff ratio
## bw 2008.65 3335.27 -1326.62 4.68 NA
## preterm 6.07 1.18 4.89 2.55 NA
## dayskidh 14.69 4.17 10.52 0.93 NA
## sex 0.51 0.50 0.01 0.02 NA
## first 0.48 0.42 0.06 0.12 NA
## age 56.67 56.35 0.32 0.16 NA
## black 0.50 0.28 0.22 0.44 NA
## hispanic 0.09 0.21 -0.12 0.41 NA
## white 0.40 0.50 -0.10 0.20 NA
## b.marr 0.43 0.69 -0.26 0.52 NA
## lths 0.43 0.30 0.13 0.27 NA
## hs 0.28 0.42 -0.14 0.31 NA
## ltcoll 0.17 0.19 -0.03 0.07 NA
## college 0.12 0.08 0.04 0.11 NA
## work.dur 0.59 0.62 -0.03 0.06 NA
## prenatal 0.96 0.99 -0.03 0.15 NA
## momage 24.44 23.75 0.69 0.12 NA
## --
##
## Balance Statistics for Matched Data
## --
## treat control unstd.diff abs.std.diff ratio
## bw 2008.65 2602.50 -593.85 2.10 2.31
## preterm 6.07 5.43 0.64 0.33 1.19
## dayskidh 14.69 10.54 4.14 0.37 1.28
## sex 0.51 0.51 0.00 0.00 --
## first 0.48 0.46 0.02 0.04 --
## age 56.67 56.09 0.59 0.30 14.66
## black 0.50 0.43 0.07 0.14 --
## hispanic 0.09 0.10 -0.01 0.04 --
## white 0.40 0.46 -0.06 0.12 --
## b.marr 0.43 0.48 -0.05 0.10 --
## lths 0.43 0.48 -0.04 0.08 --
## hs 0.28 0.32 -0.03 0.08 --
## ltcoll 0.17 0.13 0.04 0.10 --
## college 0.12 0.08 0.04 0.12 --
## work.dur 0.59 0.58 0.01 0.01 --
## prenatal 0.96 0.98 -0.02 0.12 --
## momage 24.44 23.69 0.75 0.13 0.55
## --# Reproduzindo o gráfico 10.3 (fonte: de https://avehtari.github.io/ROS-Examples/Childcare/childcare.html)
pts <- bal_nr$diff.means.raw[,4]
pts2 <- bal_nr$diff.means.matched[,4]
K <- length(pts)
idx <- 1:K
main <- 'Absolute Standardized Difference in Means'
#mar <- c(8, 6, 6, 7)
#par(mar=mar)
#maxchar <- max(sapply(cov_names, nchar))
#min.mar <- par('mar')
#mar[2] <- max(min.mar[2], trunc(mar[2] + maxchar/10)) + mar[2] + 0.5
#par(mar=mar)
pts <- rev(pts)
pts2 <- rev(pts2)
longcovnames <- rev(cov_names)
plot(c(pts,pts2), c(idx,idx),
bty='n', xlab='', ylab='',
xaxt='n', yaxt='n', type='n',
main=main, cex.main=1.2)
abline(v=0, lty=2)
points(pts, idx, cex=1)
points(pts2, idx, pch=19, cex=1)
axis(3)
axis(2, at=1:K, labels=longcovnames[1:K],
las=2, hadj=1, lty=0, cex.axis=1.2)
# Tendo criado e checado o matching, ajustamos uma regressão nesses preditores e numa dummy oara o tratamento
reg.ps <- lm (ppvtr.36 ~ treat + hispanic + black + b.marr + lths + hs +
ltcoll + work.dur + prenatal + momage + sex + first + preterm + age +
dayskidh + bw, data=matched)
display(reg.ps) # temos um efeito de 10.2 com erro padrão de 1.6## lm(formula = ppvtr.36 ~ treat + hispanic + black + b.marr + lths +
## hs + ltcoll + work.dur + prenatal + momage + sex + first +
## preterm + age + dayskidh + bw, data = matched)
## coef.est coef.se
## (Intercept) 104.46 9.57
## treat 10.28 1.59
## hispanic -16.53 2.50
## black -15.80 1.61
## b.marr 1.41 1.60
## lths -10.84 2.84
## hs -7.13 2.70
## ltcoll -4.88 2.89
## work.dur 5.46 1.49
## prenatal -3.74 3.86
## momage -0.15 0.18
## sex 0.18 1.35
## first 4.58 1.51
## preterm 0.31 0.40
## age -0.01 0.04
## dayskidh -0.22 0.06
## bw 0.00 0.00
## ---
## n = 580, k = 17
## residual sd = 16.14, R-Squared = 0.37Em um artigo de 2019 King e Nielsen sugerem que o uso do propensity score deve ser evitado, pois ele pode acabar gerando desbalanceamento, viés e maior dependência do modelo, isto é, tudo o que propõe evitar. Isso se dá por que o propensity score procura aproximar um experimento completamente aleatório. Em experimentos deste tipo alguma falta de balanceamento pode ocorrer pois a atribuição do tratamento é aleatória com relação a \(X\). Para evitar esse desbalanceamento utiliza-se o a aleatorização por blocos, isto é, o tratamento é aleatório dentro de alguma categoria (homens e mulheres, por exemplo). A aleatorização por blocos garante o máximo balancemaneto nas variáveis que compuseram os blocos. As técnicas de matching que vimos anteriormente se aproximam deste tipo de experimento, por isso, essas técnicas vão ser mais balaneceadas que o propensity score. O problema se agrava quando, ao se aproximar do experimento completamente aleatório, o porpensity score vai descartar observações não utilizadas no matching de forma também aleatória e isso, conforme King e Nielsen, gera mais desbalanceamento, enviesamento e dependência de modelo.
Regressão Descontínua
Apesar de usarmos o matching para garantir a comparação de unidades semelhantes, garantindo balanceamento e sobreposição, há uma maneira de utilizar a falta de sobreposição para estimar um efeito causal. Se uma covariante utilizada para selecionar quem recebe o tratamento for a responsável por uma clara falha na sobreposição podemos usar uma análise da descontinuidade em uma regressão nesta variável para estabelecer um efeito causal.
Pense por exemplo na adoção da urna eletrônica e seu efeito sobre o número de votos inválidos (Hidalgo, 2014). Nas eleições de 1998 apenas municípios com mais de 40 mil eleitores receberam a urna eletrônica. Assim o tratamento foi aplicado para as unidades onde \(x \ge 40 mil\) . Se assumimos que a relação entre o número de votos inválidos e a população é contínua na região de corte a descontinuidade mede o efeito do tratamento.
Como, por desenho, não temos sobreposição no ponto de corte, somos obrigados a extrapolar nosso modelo para além dos dados. O que geralmente se faz para mitigar esta extrapolação é limitar a análise a observações próximas ao ponto de corte supondo que estas unidades serão bem balanceadas com relação a possíveis variáveis de confusão.
A escolha do corte e da forma funcional tem importantes consequências como podemos ver no exemplo a seguir.
# Exemplo fictício. Tratamento para valores de x acima de 10.
x <- rnorm(200,10,1)
z <- rnorm(200,15,1)
y <- 50 + 0.5*x + 10*I(x > 10) + 2*z + rnorm(200)
plot(y~x, pch = 16, cex = 0.7, col = "darkgrey")
lines(lowess(y[x<10] ~ x[x<10]))
lines(lowess(y[x>10]~x[x>10]))
head(lowess(y[x>10]~x[x>10])$y,1) - tail(lowess(y[x<10]~x[x<10])$y,1)## [1] 10.69119mean(head(lowess(y[x>10]~x[x>10])$y,30)) - mean(tail(lowess(y[x<10]~x[x<10])$y,30))## [1] 10.72691lm <- lm(y ~ x + z)
display(lm)## lm(formula = y ~ x + z)
## coef.est coef.se
## (Intercept) 13.44 3.93
## x 4.33 0.22
## z 2.20 0.23
## ---
## n = 200, k = 3
## residual sd = 3.20, R-Squared = 0.73lm <- lm(y[x>9.5 & x<10.5]~ x[x>9.5 & x<10.5] + z[x>9.5 & x<10.5])
display(lm)## lm(formula = y[x > 9.5 & x < 10.5] ~ x[x > 9.5 & x < 10.5] +
## z[x > 9.5 & x < 10.5])
## coef.est coef.se
## (Intercept) -92.31 13.70
## x[x > 9.5 & x < 10.5] 15.10 1.33
## z[x > 9.5 & x < 10.5] 2.07 0.40
## ---
## n = 71, k = 3
## residual sd = 3.11, R-Squared = 0.73lm <- lm(y ~ x + I(x>10) + z)
display(lm)## lm(formula = y ~ x + I(x > 10) + z)
## coef.est coef.se
## (Intercept) 47.43 1.34
## x 0.49 0.10
## I(x > 10)TRUE 10.15 0.22
## z 2.17 0.06
## ---
## n = 200, k = 4
## residual sd = 0.92, R-Squared = 0.98lm1 <- lm(y[x>10]~x[x>10] + z[x>10])
lm2 <- lm(y[x<10]~x[x<10] + z[x<10])
curve(coef(lm1)[1] + coef(lm1)[2]*x + coef(lm1)[3]*mean(z), from = 10, to = 13, add = T, col = 2)
curve(coef(lm2)[1] + coef(lm2)[2]*x + coef(lm2)[3]*mean(z), from = 7, to = 10, add = T, col = 2)
(coef(lm1)[1] + coef(lm1)[2]*10 + coef(lm1)[3]*mean(z)) - (coef(lm2)[1] + coef(lm2)[2]*10 + coef(lm2)[3]*mean(z))## (Intercept)
## 10.134dt1 <- data.frame(rep(1,10), seq(10,10.5,length.out = 10), rep(mean(z),10))
dt2 <- data.frame(rep(1,10),seq(9.5,10,length.out = 10), rep(mean(z),10))
mean(predict(lm1,dt1)) - mean(predict(lm2,dt2))## [1] 11.2522y <- 50 + 0.5*x + 0.02*x^2 -0.008*x^3 + 10*I(x > 10) + 2*z + rnorm(200)
plot(y~x, pch = 16, cex = 0.7, col = "darkgrey")
lines(lowess(y[x<10] ~ x[x<10]))
lines(lowess(y[x>10]~x[x>10]))
head(lowess(y[x>10]~x[x>10])$y,1) - tail(lowess(y[x<10]~x[x<10])$y,1)## [1] 9.913642mean(head(lowess(y[x>10]~x[x>10])$y,30)) - mean(tail(lowess(y[x<10]~x[x<10])$y,30))## [1] 9.503712lm <- lm(y ~ x + I(x>10) + z) # modelo reduzido
display(lm)## lm(formula = y ~ x + I(x > 10) + z)
## coef.est coef.se
## (Intercept) 64.86 1.40
## x -1.68 0.11
## I(x > 10)TRUE 10.22 0.23
## z 2.03 0.07
## ---
## n = 200, k = 4
## residual sd = 0.96, R-Squared = 0.96lm <- lm(y ~ poly(x,3) + I(x>10) + z) # modelo correto
display(lm)## lm(formula = y ~ poly(x, 3) + I(x > 10) + z)
## coef.est coef.se
## (Intercept) 47.87 0.99
## poly(x, 3)1 -22.74 1.73
## poly(x, 3)2 -3.30 0.93
## poly(x, 3)3 -2.62 1.10
## I(x > 10)TRUE 9.82 0.26
## z 2.05 0.07
## ---
## n = 200, k = 6
## residual sd = 0.92, R-Squared = 0.96lm <- lm(y ~ poly(x,3) + I(x>10) + I(x>10)*poly(x,3) + z) # modelo saturado
display(lm)## lm(formula = y ~ poly(x, 3) + I(x > 10) + I(x > 10) * poly(x,
## 3) + z)
## coef.est coef.se
## (Intercept) 47.04 1.70
## poly(x, 3)1 -38.26 28.04
## poly(x, 3)2 -17.08 20.99
## poly(x, 3)3 -10.03 8.35
## I(x > 10)TRUE 10.88 1.69
## z 2.06 0.07
## poly(x, 3)1:I(x > 10)TRUE 8.87 34.09
## poly(x, 3)2:I(x > 10)TRUE 17.73 26.22
## poly(x, 3)3:I(x > 10)TRUE 7.46 10.66
## ---
## n = 200, k = 9
## residual sd = 0.92, R-Squared = 0.96lm1 <- lm(y[x>10]~x[x>10] + z[x>10])
lm2 <- lm(y[x<10]~x[x<10] + z[x<10])
curve(coef(lm1)[1] + coef(lm1)[2]*x + coef(lm1)[3]*mean(z), from = 10, to = 13, add = T, col = 2)
curve(coef(lm2)[1] + coef(lm2)[2]*x + coef(lm2)[3]*mean(z), from = 7, to = 10, add = T, col = 2)
(coef(lm1)[1] + coef(lm1)[2]*10 + coef(lm1)[3]*mean(z)) - (coef(lm2)[1] + coef(lm2)[2]*10 + coef(lm2)[3]*mean(z))## (Intercept)
## 10.17807dt1 <- data.frame(rep(1,10), seq(10,10.5,length.out = 10), rep(mean(z),10))
dt2 <- data.frame(rep(1,10),seq(9.5,10,length.out = 10), rep(mean(z),10))
mean(predict(lm1,dt1)) - mean(predict(lm2,dt2))## [1] 7.646167Variáveis instrumentais
Até agora nossa estatégia para assegurar a premissa de independência condicional dos resultados potenciais com relação ao tratamento foi controlar por variáveis de confusão, \((y_{0i},y_{1i}) \perp T_{i} \mid X_{i}\). Com a aleatorização ou experimentos controlados essa premissa era diretamente assegurada, pois a natureza aleatória do tratamento garante a independência. Com dados observacionais isso só ocorre com a inclusão das variáveis omitidas garantindo o balanceamento e a sobreposição destas variáveis nos grupos tratado e controle.
O que acontece quando não podemos controlar por variáveis omitidas? Vamos supor, usando o exemplo em Angrist e Pischke, que estamos interessados no efeito da escolaridade na renda. Variáveis que podem influir tanto na escolaridade quanto na renda são o ambiente familiar, as abilidades inatas, a motivação pessoal etc. Essas variáveis não são facilmente mensuradas, mas sabemos que elas podem afetar substantivamente nossa conclusão sobre o efeito da escolaridade na renda. Como podemos lidar com isso?
Uma estratégia é substituirmos a escolaridade por uma variável fortemente correlacionada com ela, mas que seja independente da abilidade, motivação e ambiente familiar. Ao mesmo tempo, essa variável só deve impactar a renda por meio do seu impacto na escolaridade. Chamamos uma variável com estas características de variável instrumental.
Usando a notação de Angrist e Pischke a relação entre a variável instrumental, o tratamento e a variável de interesse se dá em duas etapas. Primeiro temos a relação entre o tratamento (\(T\)) e a variávels instrumental (\(z\))1, o primeiro estágio e, então a regressão da variável de interesse (\(Y\)) no instrumento:
\[ \begin{align} T_{i} &= X'\pi_{10} + \pi_{11}z_{i} + \xi_{1i} \\ Y_{i} &= \alpha'X_{i} + \rho[X'\pi_{10} + \pi_{11}z_{i} + \xi_{1i}] + \eta_{i} \\ &= X_{i}[\alpha + \rho\pi_{10}] + \rho\pi_{11}z_{i} + [\rho\xi_{1i} + \eta_{i}] \\ &= X'\pi_{20} + \pi_{21}z_{i} + \xi_{2i} \end{align} \]
Na última linha da equação \(\pi_{20} = \alpha + \rho\pi_{10}\), \(\eta_{2i} = \rho\eta_{1i} + \eta_{i}\) e \(\pi_{21} = \rho\pi_{11}\), Desta última igualdade resulta que o efeito causal do tratamento é dado por \(\rho = \frac{\pi_{21}}{\pi_{11}}\).
Se rearranjarmos a última linha da equação temos
\[Y_{i} = \alpha'X_{i} + \rho[X'\pi_{10} + \pi_{11}z_{i}] + \xi_{2i}\] O que mostra que o efeito causal de \(T\) em \(Y\) pode ser obtido fazendo a regressão de \(Y_{i}\) em \(X_{i}\) e nos valores ajustados da regressão de \(T_{i}\) em \(X_{i}\) \(z_{i}\). Na prática obtemos os valores ajustados
\[ \hat{T_{i}} = X' \pi_{10} + \pi_{11} z_{i} \] e substituimos na equação
\[Y_{i} = \alpha'X_{i} + \rho \hat{T_{i}} + [\eta_{i} + \rho(T_{i} - \hat{T_{i}})]\] Essa equação é conhecida como 2SLS (two-stage least square) pois envolve a estimativa dos valores preditos \(\hat{t}\) num primeiro estágio e, então, o uso desses valores para obter \(Y_{i}\) num segundo estágio. Em geral, devido ao cálculo dos erros padrão e de outros ajustes (ver Gelman e Hill pg. 224), a estimativa é feita utilizando-se algum agoritmo específico como o ivreg() de Fox, Kleiber e Zeileis (2020)
library(ivreg)## Registered S3 methods overwritten by 'car':
## method from
## influence.merMod lme4
## cooks.distance.influence.merMod lme4
## dfbeta.influence.merMod lme4
## dfbetas.influence.merMod lme4data("SchoolingReturns", package = "ivreg")
summary(SchoolingReturns[, 1:8])## wage education experience ethnicity smsa
## Min. : 100.0 Min. : 1.00 Min. : 0.000 other:2307 no : 864
## 1st Qu.: 394.2 1st Qu.:12.00 1st Qu.: 6.000 afam : 703 yes:2146
## Median : 537.5 Median :13.00 Median : 8.000
## Mean : 577.3 Mean :13.26 Mean : 8.856
## 3rd Qu.: 708.8 3rd Qu.:16.00 3rd Qu.:11.000
## Max. :2404.0 Max. :18.00 Max. :23.000
## south age nearcollege
## no :1795 Min. :24.00 no : 957
## yes:1215 1st Qu.:25.00 yes:2053
## Median :28.00
## Mean :28.12
## 3rd Qu.:31.00
## Max. :34.00m_ols <- lm(log(wage) ~ education + poly(experience, 2) + ethnicity + smsa + south,
data = SchoolingReturns)
summary(m_ols)##
## Call:
## lm(formula = log(wage) ~ education + poly(experience, 2) + ethnicity +
## smsa + south, data = SchoolingReturns)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.59297 -0.22315 0.01893 0.24223 1.33190
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.259820 0.048871 107.626 < 2e-16 ***
## education 0.074009 0.003505 21.113 < 2e-16 ***
## poly(experience, 2)1 8.931699 0.494804 18.051 < 2e-16 ***
## poly(experience, 2)2 -2.642043 0.374739 -7.050 2.21e-12 ***
## ethnicityafam -0.189632 0.017627 -10.758 < 2e-16 ***
## smsayes 0.161423 0.015573 10.365 < 2e-16 ***
## southyes -0.124862 0.015118 -8.259 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3742 on 3003 degrees of freedom
## Multiple R-squared: 0.2905, Adjusted R-squared: 0.2891
## F-statistic: 204.9 on 6 and 3003 DF, p-value: < 2.2e-16# Como discutido acima escolaridade (education) é endógena, assim como experiencia (que é função linear de escolaridade)
# Para lidar com isso usamos o instrumento "proximidade" geográfica a uma escola para escolaridade e "idade" para experiencia
m_iv <- ivreg(log(wage) ~ education + poly(experience, 2) + ethnicity + smsa + south |
nearcollege + poly(age, 2) + ethnicity + smsa + south,
data = SchoolingReturns)
summary(m_iv)##
## Call:
## ivreg(formula = log(wage) ~ education + poly(experience, 2) +
## ethnicity + smsa + south | nearcollege + poly(age, 2) + ethnicity +
## smsa + south, data = SchoolingReturns)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.82400 -0.25248 0.02286 0.26349 1.31561
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.48522 0.67538 6.641 3.68e-11 ***
## education 0.13295 0.05138 2.588 0.009712 **
## poly(experience, 2)1 9.14172 0.56350 16.223 < 2e-16 ***
## poly(experience, 2)2 -0.93810 1.58024 -0.594 0.552797
## ethnicityafam -0.10314 0.07737 -1.333 0.182624
## smsayes 0.10798 0.04974 2.171 0.030010 *
## southyes -0.09818 0.02876 -3.413 0.000651 ***
##
## Diagnostic tests:
## df1 df2 statistic p-value
## Weak instruments (education) 3 3003 8.008 2.58e-05 ***
## Weak instruments (poly(experience, 2)1) 3 3003 1612.707 < 2e-16 ***
## Weak instruments (poly(experience, 2)2) 3 3003 174.166 < 2e-16 ***
## Wu-Hausman 2 3001 0.841 0.432
## Sargan 0 NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4032 on 3003 degrees of freedom
## Multiple R-Squared: 0.1764, Adjusted R-squared: 0.1747
## Wald test: 148.1 on 6 and 3003 DF, p-value: < 2.2e-16O coeficiente causal encontrado acima, \(\rho\), era constante para todas as observações. Assumir que o efeito é o mesmo para todos é algo um tanto irrealista, algumas unidades podem responder melhor a um tratamento do que outras, por exemplo. Uma maneira de lidar com isso é permitir que o efeito varie com os indivíduos, isto é, agregamos um subscrito à \(\rho_{i}\). Uma consequência disso é que se o efeito é específico a uma determinada população ele pode ter validade interna, mas não ter muita validade externa, isto é, podemos estar seguro de estimar o efeito correto do tratamento em uma amostra, mas podemos fracassar ao utilizar esse efeito para predizer novoso valores. Já enfrentamos isso quando falamos de validação, sobre o trade-off entre viés e variância, MSE treino e MSE teste.
Quando temos efeitos heterogêneos devemos assumir algumas premissas sobre nosso instrumento:
- 1: Independência \((Y_{i}(T_{0i},0), Y_{i}(T_{1i},1), T_{1i}, T_{0i}) \perp z_{i}\);
- 2: Exclusão de Restrição \(Y_i(d,0) = Y_i(d,1) \equiv Y_{di} \text{ para } d = 0,1\);
- 3: Monotonicidade $ T_{1i} - T_{0i} 0 i$;
- 4: Existência do primeiro estágio \(E(T_{1i} - T_{0i}) \neq 0\)
A premissa da independência exige ignorabilidade não só do tratamento, como também do instrumento, o que é um tanto razoável. A premissa da existência do primeiro estágio nos diz que o instrumento deve ter algum efeito sobre o tratamento. Utilizando o exemplo de Gelman e Hill, do efeito de um programa de TV sobre o desempenho escolar de crianças, o tratamento - exposição ao programa - deve ser afetado pelo instrumento - encorajamento a assistir ao programa. Naquele caso 90% dos encorajados assistiam regularmente ao programa contra 55% dos não encorajados. O instrumento, portanto, só afetou 35% dos encorajados, os outrs 55% assistiriam de qualquer maneira.
O efeito do tratamento só se aplica aos que eram afetados pelo instrumento. Isso divide as observações em três grupos:
- Always takers
- Compliers
- Never takers
Os compliers (aderentes) são aqueles com base nos quais podemos identificar o efeito causal. Os always takers e never takers não alteram o resultado potencial seja qual for o tratamento ou instrumento. A restrição de exclusão exige que o tratamento ou instrumento não altere o resultado desses não aderentes (no exemplo de Gelman e Hill o encorajamento não deve alterar o desempenho dos não aderentes). Além disso o efeito do instrumento só pode ter uma direção (monotonicidade). O encorajamento para ver o programa não pode fazer com que o indivíduo sucetível (aderente) deixe de ver o programa.
Com essas premissas justificamos o teorema do LATE (Local Average Treatment Effect). Supondo que o efeito causal é uma dummy (estimador de Wald) podemos expressá-lo como a razão da diferença nos efeitos potenciais esperados com relação ao tratamento e ao instumento
\[ \frac{E[Y_{i} \mid z_{i} = 1] - E[Y_{i} \mid z_{i} = 0]}{E[T_{i} \mid z_{i} = 1] - E[T_{i} \mid z_{i} = 0]} = E[Y_{1i} - Y_{0i} \mid T_{1i} > T_{0i}] \\ = E[\rho_i \mid \pi_{1i} > 0 ] \]
O que mostra que \(\rho_{i}\) é o efeito causal quando o efeito de \(Z_{i}\) em \(T_{i}\) existe e é maior que zero para todos. Mais uma vez, isso quer dizer que estamos estimando o efeito do tratamento apenas nos que foram afetados pelo instrumento, ou o efeito médio do tratamento nos tratados (ATT). Utilizando o estimador de Wald, podemos estimar a proporção de aderentes por meio do coeficiente do instrumento na regressão do tratamento neste instrumento. A regressão da variável de interesse no instrumento nos dá o efeito do tratamento em todos os tratados, o que se chama de intent to treat (ITT). A razão entre esse dois coeficientes é o ATT. Enquanto o ITT nos dá uma ideia do efeito do tratamento levando em conta a possibilidade de não adesão ele só pode ser generalizado se na próxima aplicação do tratamento tivermos o mesmo número de aderentes. Já o ATT nos dá o efeito do tratamento, sem se preocupar com o número de aderentes. Cada uma destas estimativas olha para o mesmo problema de ângulos diferentes.
Voltando à questão da validade interna e validade externa, o ATT nos dá a efetividade de nosso tratamento. No caso das crianças encorajadas a assistir ao programa de TV importa apenas se o desempenho destas crianças melhorou. O ATT, ou o ITT, não nos diz nada sobre o efeito de futuros encorajamentos para assistir a outros programas, não há valor preditivo nestes estimadores. Para poder prever é necessário saber por que o tratamento leva a determinado resultado, isto é, é necessário teoria (ver Angrist e Pischke, pg. 156).
No jargão da inferência causal, as variáveis \(T\) e \(Y\) são chamadas de variáveis endógenas, pois são definidas dentro do sistema de equações. As variáveis \(X_{i}\) e \(z_{i}\) são chamadas de variáveis exógenas.↩