A última aula será dedicada a demostrar a utilização de modelos de Regressão com Descontinuidade (RDD). Utilizaremos como exemplo o conjunto de dados para mensurar o efeito da idade legal em acidentes de carro (Carpenter & Dobkin 2011).
Pergunta de pesquisa: Qual o impacto de ter idade para beber alcool sobre o número de mortes no transito?
-Em um esforço para lidar com os problemas sociais e de saúde pública associados ao consumo de álcool por menores de idade, um grupo de presidentes de universidades americanas fez lobby junto aos estados para que o limite mínimo legal de idade para consumir álcool (MLDA, na sigla em inglês) retorne ao limiar da era do Vietnã, que era de 18 anos. A teoria por trás desse esforço (conhecido como Iniciativa Amethyst) é que permitir o consumo legal de álcool a partir dos 18 anos desencoraja o consumo excessivo e promove uma cultura de consumo maduro de bebidas alcoólicas. Isso contrasta com a visão tradicional de que o MLDA de 21 anos, embora uma ferramenta direta e imperfeita, reduz o acesso dos jovens ao álcool, prevenindo assim alguns danos. Angrist e Pischke (2014).
Para responder a pergunta Carpenter e Dobkin (2011) buscaram dados nos Estados Unidos, em que a idade para dirigir começa aos 16 anos mas o consumo de bebida só é legalizado a partir dos 21 anos. Assim, se torna uma oportunidade para a aplicação de RDD a fim de avaliar o impacto da MLDA nas mortes em transito.
Code
# Pacotes utilizadoslibrary(tidyverse)library(rdrobust)library(readr)library(sjPlot)library(broom)library(estimatr) set.seed(42)# Download do banco de dadosmlda_df <- readr::read_csv("https://raw.githubusercontent.com/seramirezruiz/stats-ii-lab/master/Session%206/data/mlda.csv")mlda_df <- mlda_df |>select(c(outcome, agecell, treatment, forcing))sample_n(mlda_df, 5) |>tab_df(title ="Estrutura do banco")
Estrutura do banco
outcome
agecell
treatment
forcing
31.83
22.03
1
1.03
35.83
19.07
0
-1.93
36.32
21.04
1
0.04
30.58
19.81
0
-1.19
31.54
21.95
1
0.95
No banco exposto, a variável \(outcome\) é o número de mortes no trânsito por 100 mil habitantes, \(agecell\) é a faixa de idade, \(treatment\) indica que a faixa é maior do que 21 anos e \(forcing\) é a variável padronizada a fim de 21 anos ser equivalente a zero e a variação se manter entre 2 e menos 2, sendo essa a variável de atribuição.Podemos imaginar que a variável \(D\) define o recebimento ou não do tratamento, de modo que:
\(D = 1\) se \(X\geq c\)
\(D = 0\) se \(X< c\)
Em que \(X\) é a variável de atribuição e \(C\) é o ponto de corte. A RDD é baseada em uma hipóteses:
Hipótese 1: Todos os outros fatores (que não \(D\)) estão evoluindo de maneira suave em torno de \(c\).
Essa deriva em duas outras:
Hipótese 2.1: Funções de regressões condicionais contÌnuas: \(E[Y (0)| X = x]\) e E \([Y (1)| X = x]\) são contínuas em \(x\).
Hipótese 2.2: Funções de regressões condicionais contÌnuas: \(F_{Y(0)|x}(y|x)\) e \(F_{Y(1)|x}(y|x)\) são funções continuas em \(x\) para todo \(y\).
É possível supor que o comportamento dos indivíduos entre 19 e 23 anos não sofre grandes alterações, devido a proxímidade da idade, de modo que atendemos as hipóteses necessárias. O segundo passo é determinar qual desing Caso \(c\) seja claro, ou seja um ponto específico marque o corte e nada antes dele receba o tratamento e todos depois dele recebam, temos um modelo com o set up de um modelo \(sharp\). Se a linha de corte não é clara, ou seja, existam individuos que recebem e não recebem o tratamento na mesma etapa, temos um setup \(fuzzy\).
Code
library(tidyverse)library(gridExtra)# Criando a base de dados simulada com mais pontosdata_fuzzy <-tibble(d =c(rep(50,x=0),rep(50,x=1)),x =c(rnorm(40,10, 2),rnorm(20,15,2),rnorm(40,20, 2)),eps =rnorm(100,0,15),y = x*40+ d*60+ eps)# Visualização Fuzzy setg1 <-ggplot(data_fuzzy, aes(x = x, y =as.factor(d))) +geom_point() +geom_vline(xintercept =15, color ="red", linetype ="dashed") +labs(x ="X", y ="D", title ="Dados em um fuzzy set")+theme_bw()g1.1<-ggplot(data_fuzzy, aes(x = x, y = y, color =as.factor(d))) +geom_point() +geom_vline(xintercept =15, linetype ="dotted") +theme_bw()+labs(color ="D")# Visualização Sharp setg2 <- data_fuzzy |>filter(x >15& d==1|x <15& d==0) |>ggplot(aes(x = x, y =as.factor(d))) +geom_point() +geom_vline(xintercept =15, color ="red", linetype ="dashed") +labs(x ="X", y ="D", title ="Dados em um sharp set")+theme_bw()g2.1<- data_fuzzy |>filter(x >15& d==1|x <15& d==0) |>ggplot(aes(x = x, y = y, color =as.factor(d))) +geom_point() +geom_vline(xintercept =15, linetype ="dotted") +labs(color ="D")+theme_bw()grid.arrange(g2,g1,g2.1,g1.1, ncol =2)
O primeiro passo é avaliar se aplicaremos uma RDD do tipo \(sharp\) ou do tipo \(fuzzy\). No caso, o corte na lei é claro, legalizando o consumo de alcool no momento em que o indivíduo completou 21 anos. Podemos vizualizar a propabilidade de se encaixar no tratamento a partir do seguinte gráfico:
Code
ggplot(mlda_df, aes(x = agecell, y = treatment,color =factor(treatment))) +geom_point(size =1.5, alpha =0.5,position =position_jitter(width =0, height =0.25, seed =1234)) +labs(x ="Idade", y ="Probabilidade de ser tratado") +scale_color_discrete(name =" ", labels =c("Não possui idade legal para beber", "Possui idade legal para beber")) +geom_vline(xintercept =21, linetype ="dotted") +theme_bw()+theme(legend.position ="bottom")
Assim, o uso de uma RDD \(sharp\) é o indicado. A partir disso, analisaremos então o comportamento dos dados a fim de entender qual modelo melhor capturará a variação dos dados.
Code
ggplot(mlda_df, aes(x = agecell, y = outcome)) +geom_point() +geom_vline(xintercept =21, linetype ="dotted") +labs(title ="Análise exploratória",x ="Variável forçada (Idade)",y ="Taxa de mortalidade por acidentes\n de transito (por 100,000)") +theme_minimal()
Assumimos inicialmente que a relação entre \(X\) e \(Y\) seja linear, podemos estimar a relação entre o efeito do tratamento \(\tau\) por meio da seguinte regressao linear.
\[Y = \alpha + D_{\tau} + \beta X + \epsilon\]
Se todos os fatores variam suavemente de acorco com \(X\), então \(\beta\) seria uma boa estimativa do valor de \(X\) para quem está acima da linha de corte e \(\alpha\) para aqueles abaixo da linha de corte, de modo que o efeito causal seria dado por \(\beta X\) - \(\alpha\). Para supormos isso, é necessário que todos os fatores que afetam \(Y\) variem suavimente com \(X\). Caso isso não se confirme, \(\tau\) será estimado com viés. Assim, primeiramente, rodaremos um modelo linear com retas comuns. Observe que a variável \(forcing\) é CENTRALIZADA em 0 (idade 21) e representa a distância em anos a partir dos 21 anos, enquanto o tratamento é apenas binário, acima/abaixo de 21 anos.
Code
linear_common_slope <-lm(outcome ~ treatment + forcing, data = mlda_df)linear_common_slope |>tab_model()
outcome
Predictors
Estimates
CI
p
(Intercept)
29.36
28.49 – 30.22
<0.001
treatment
4.53
2.99 – 6.08
<0.001
forcing
-3.15
-3.83 – -2.47
<0.001
Observations
48
R2 / R2 adjusted
0.703 / 0.689
De acordo com nossas suposições para modelos lineares com inclinação comum, consideramos que o efeito do tratamento não depende da variável de forçamento. Podemos formalizar esse modelo da seguinte forma:
Portanto, podemos afirmar que, com base em nosso \(\beta_1\), podemos esperar 4,53 mortes adicionais por acidentes de trânsito a cada 100.000 pessoas para aqueles que podem legalmente consumir bebidas alcoólicas. Também observamos que, para cada ano de aumento na idade, o número de mortes esperadas por 100.000 pessoas diminui em 3,15 (\(\beta_2\). Podemos extrair os valores preditos a fim de recriar o fit. Uma das melhores formas de vizualizar os resultados é dividindo o gráfico e duas janelas (bins) e plotar as linhas da regressão em ambos os casos:
Code
mlda_df$yhat_linear <-predict(linear_common_slope) # Criando uma variável com o modelo predito. mlda_df %>%ggplot(aes(x = forcing, y = yhat_linear, col =factor(treatment))) +geom_point(aes(x = forcing, y = outcome,col =factor(treatment))) +geom_vline(xintercept =0, linetype ="dotted") +labs(title ="Modelo linear com reta comum",x ="Variável forçada (Idade)",y ="Taxa de mortalidade por acidentes\n de transito (por 100,000)") +geom_line(data = mlda_df[mlda_df$forcing >=0,], color ="#cc0055",size =1) +geom_line(data = mlda_df[mlda_df$forcing <0,], color ="#696969", size =1) +scale_color_manual(name ="",values =c("#696969", "#cc0055"),labels =c("Controle", "Tratamento")) +theme_bw()
Agora, reestimaremos essa regressão permitindo que haja variação entre as retas da regressão. Isso pode ser feito por meio da introdução de um termo interativo no modelo.
linear_different_slope <-lm(outcome ~ treatment*forcing, data = mlda_df)linear_different_slope |>tab_model()
outcome
Predictors
Estimates
CI
p
(Intercept)
29.93
28.86 – 31.00
<0.001
treatment
4.53
3.02 – 6.05
<0.001
forcing
-2.57
-3.51 – -1.63
<0.001
treatment × forcing
-1.16
-2.49 – 0.17
0.085
Observations
48
R2 / R2 adjusted
0.722 / 0.703
Portanto, podemos afirmar que, dado \(\beta_1\), podemos esperar 4,53 mortes adicionais por acidentes de veículos a cada 100.000 pessoas para aqueles que podem legalmente consumir bebidas alcoólicas no limite de idade. Agora temos duas inclinações diferentes para mudanças na idade. Para indivíduos com menos de 21 anos, um aumento de um ano na idade resultaria, em média, em 2,57 mortes a menos por acidentes de veículos (\(\beta_2\)). Para aqueles que têm idade legal para beber, esperaríamos 3,73 mortes a menos por 100.000 pessoas para cada aumento de um ano na idade (\(\beta_2X_i+\beta_3\tilde{X}_i\) = - 2,57 + (-1,16). Além disso, o modelo se adequou um pouco melhor aos dados, com um \(R^2\) de 0,722, sendo 0,019 maior que o modelo anterior.
Code
mlda_df %>%ggplot(aes(x = forcing, y = outcome, col =factor(treatment))) +geom_point() +geom_vline(xintercept =0, linetype ="dotted") +geom_smooth(method ="lm", se = F) +labs(title ="Modelo linear com retas variando",x ="Variável forçada (Idade)",y ="Taxa de mortalidade por acidentes\n de transito (por 100,000)")+scale_color_manual(name ="",values =c("#696969", "#cc0055"),labels =c("Controle", "Tratamento")) +theme_bw()
Podemos também deixar o tratamento variar ao longo do tratamento, respeitando as premissas de modelos não lineares. Assim, os efeitos do tratamento \(D_i\) podem variar ao longo da variável \(X_i\). No caso, estimaremos um modelo com interações quadráticas. O modelo é especificado da seguinte forma:
quadratic <-lm(outcome ~ forcing +I(forcing^2) +# I tells R to interpret "as is" treatment +I(forcing * treatment) +I((forcing^2) * treatment),data = mlda_df)quadratic |>tidy() |>tab_df()
term
estimate
std.error
statistic
p.value
(Intercept)
29.81
0.82
36.51
0.00
forcing
-2.93
1.91
-1.53
0.13
I(forcing^2)
-0.19
0.94
-0.20
0.84
treatment
4.66
1.15
4.04
0.00
I(forcing * treatment)
-0.82
2.71
-0.30
0.76
I((forcing^2) * treatment)
0.20
1.33
0.15
0.88
Portanto, podemos dizer que, considerando nosso \(\tau\), podemos esperar 4,66 mortes adicionais por acidentes de veículos a cada 100.000 pessoas para aqueles que podem legalmente consumir bebidas alcoólicas no limite de idade. Também poderíamos calcular o valor esperado de \(Y\) em diferentes níveis de \(X\).
Digamos que queiramos saber qual é o valor esperado para pessoas de 23 anos. Como nossa variável de forçamento é 0 aos 21 anos de idade, podemos considerar 23 como 2. Além disso, pessoas de 23 anos estão acima da idade legal mínima para beber, portanto, para elas, o valor de \(D\) é 1.
Com base nisso, esperaríamos uma taxa de mortalidade por acidentes de veículos de 27,01 a cada 100.000 pessoas para pessoas de 23 anos. Podemos plotar os resultados do modelo.
Code
mlda_df$yhat_quadratic <-predict(quadratic) mlda_df %>%ggplot(aes(x = forcing, y = yhat_quadratic, #note predicted ycol =factor(treatment))) +geom_point(aes(x = forcing, y = outcome, col =factor(treatment))) +geom_vline(xintercept =0, linetype ="dotted") +labs(title ="Modelo quadratico",x ="Variável forçada (Idade)",y ="Taxa de mortalidade por acidentes\n de transito (por 100,000)") +geom_line(data = mlda_df[mlda_df$forcing >=0,], color ="#cc0055", size =1) +geom_line(data = mlda_df[mlda_df$forcing <0,], color ="#696969", size =1) +scale_color_manual(name ="",values =c("#696969", "#cc0055"),labels =c("Controle", "Tratamento")) +theme_minimal()
Podemos calcular o efeito médio causal local (Local Avarage Treatment Effect) para aqueles com mais de 21 anos por meio de uma estimação dos pontos por meio de um regressão polinomial local com correção de viés.
Code
llr <- rdrobust::rdrobust(mlda_df$outcome, mlda_df$forcing, c =0,kernel ="tri",bwselect ="mserd")summary(llr)
Sharp RD estimates using local polynomial regression.
Number of Obs. 48
BW type mserd
Kernel Triangular
VCE method NN
Number of Obs. 24 24
Eff. Number of Obs. 6 6
Order est. (p) 1 1
Order bias (q) 2 2
BW est. (h) 0.487 0.487
BW bias (b) 0.738 0.738
rho (h/b) 0.660 0.660
Unique Obs. 24 24
=============================================================================
Method Coef. Std. Err. z P>|z| [ 95% C.I. ]
=============================================================================
Conventional 4.901 2.059 2.380 0.017 [0.864 , 8.937]
Robust - - 1.881 0.060 [-0.198 , 9.674]
=============================================================================
Os resultados apontam para um efeito não significante do tratamento para aqueles com mais de 21 anos. Podemos plotar esses resultados da seguinte forma:
Code
rdplot(mlda_df$outcome, mlda_df$forcing, c =0,kernel ="tri",title ="Morte em acidentes no transito",x.label ="Idade a partir de 21",y.label ="Morte em acidentes no transito (per 100,000)")
2 Exemplo de um fuzzy set
Para dar apenas um exemplo de como operar um modelo com \(fuzzy\) set. Nesse segundo exemplo, testaremos o efeito de ter aulas particulares sobre o exame de conclusão de um curso. Nesse caso, alunos que no exame de entrada tiraram menos do que 70% tem direito a uma monitoria gratuita. Assim, teoricamente temos um \(c\) bem especificado (70%), mas o tratamento está distribuido em ambas as partes do corte. O banco está nessa estrutura
Code
tutoring <-read.csv("./tutoring_program_fuzzy.csv")sample_n(tutoring, 5) |>tab_df(title ="Estrutura do banco")
Estrutura do banco
id
entrance_exam
tutoring
tutoring_text
exit_exam
756
89.26
TRUE
Tutor
79.24
981
95.91
FALSE
No tutor
61.90
404
96.92
FALSE
No tutor
74.02
715
85.60
FALSE
No tutor
64.01
225
70.21
FALSE
No tutor
70.14
Contudo, ao fazer uma análise visual dos dados, vemos que não é isso que acontece:
Code
library(ggpubr)g1 <-ggplot(tutoring, aes(x = entrance_exam, y = tutoring_text, color = entrance_exam <=70)) +geom_point(size =1.5, alpha =0.5,position =position_jitter(width =0, height =0.25, seed =1234)) +geom_vline(xintercept =70) +labs(x ="Resultado no exame de entrada", y ="Participou do programa de monitoria") +guides(color ="none")+theme_bw()g2 <-ggplot(tutoring, aes(x = entrance_exam, y = exit_exam, color = tutoring)) +geom_point(size =1.5, alpha =0.5,position =position_jitter(width =0, height =0.25, seed =1234)) +geom_vline(xintercept =70) +labs(x ="Resultado no exame de entrada", y ="Resultado no exame de saída") +guides(color ="none")+theme_bw()ggarrange(g2,g1)
Nesse caso, 36 pessoas que deveriam ter a tutoria não a receberam e 116 pessoas que não eram elegiveis a monitoria que a receberam. Assim, o setup deixa de ser \(sharp\) e passa a ser \(fuzzy\). Assim, ao inves da propabilidade de receber o tratamento ser \(D = 1\) se \(X\geq c\) e \(D = 0\) se \(X< c\), temos apenas um pequeno aumento dessa probabilidade de um grupo para outro. Assim, utilizaremos um modelo análogo a uma variável instrumental para dividir o salto de \(Y\) sobre \(X\) com o salto da regressão do indicador do tratamento em \(X\). Para que isso seja atendido, assumimos uma terceira hipótese:
Hipótese 3: Monotonicidade - : Os valores de X que passam o ponto de corte não podem simultaneamente causar algumas unidades serem tratadas e outras rejeitarem o tratamento
Aceitando essa hipótese, dividimos nossos casos em três grupos: Compliers que respeitam o corte estimado por \(c\); Nevertakers que não recebem o tratamento, independente de estarem elegiveis ou não; e Always-takers que recebem o tratamento independente de serem elegíveis ou não.
Da mesma maneira que em variável instrumental, tinhamos que verificar que existia uma relção forte entre a variável endógena e o instrumento; neste caso temos que verifícar que de fato existe uma descontinuidade na expectativa de \(D\) condicional a \(X\). Assim, podemos formalizar em:
Que pode ser estimada por uma regressão linear de 2 estágios usando \(1{x_i-c\leq0}\) como instrumento para \(D_i\). Primeiramente, centralizaremos a variável \(X\), nesse caso a nota de entrada, em zero como estava no exemplo anterior. Logo em seguida, estimaremos a regressão em dois estágios da seguinte maneira:
Aqui, o parametro \(tutoringTRUE\) indica o tamanho do efeito, que é de 9,74 pontos no exame final. Caso tentemos estimar o efeito em uma OLS as estimativas estariam enviesadas pois estamos desconsiderando os diferentes grupos e supondo que todos são compliers. Podemos ver a diferença dos modelos na tabela a seguir:
Sem o uso de instrumentos, o modelo estima um efeito de quse dois pontos a mais da tutoria. Podemos também estimar modelos não paramétricos para estimar os resultados, do mesmo modo que fizemos com o modelo \(sharp\)
Fuzzy RD estimates using local polynomial regression.
Number of Obs. 1000
BW type mserd
Kernel Triangular
VCE method NN
Number of Obs. 238 762
Eff. Number of Obs. 170 347
Order est. (p) 1 1
Order bias (q) 2 2
BW est. (h) 12.985 12.985
BW bias (b) 19.733 19.733
rho (h/b) 0.658 0.658
Unique Obs. 238 762
First-stage estimates.
=============================================================================
Method Coef. Std. Err. z P>|z| [ 95% C.I. ]
=============================================================================
Conventional -0.708 0.073 -9.751 0.000 [-0.850 , -0.565]
Robust - - -8.424 0.000 [-0.907 , -0.565]
=============================================================================
Treatment effect estimates.
=============================================================================
Method Coef. Std. Err. z P>|z| [ 95% C.I. ]
=============================================================================
Conventional 9.683 1.893 5.116 0.000 [5.973 , 13.393]
Robust - - 4.258 0.000 [5.210 , 14.095]
=============================================================================
Assim, estimamos um efeito similar a partir de métodos não paramétricos, mantendo a significancia estatística dos resultados.
3 Limitações do método
Apesar de ser uma poderosa técnica para inferência causal, a RDD possui algumas limitações que desafiam sua validade interna e externa.
Validade interna: Nós assumimos que todos os fatores não observáveis variam de acordo com \(X\), contudo não podemos validar essa premissa por meio de testes estatísticos. Assim, é essencial entender como \(X\) foi determinado a fim de obter estimativas não viesadas.
Validade externa: Devido ao desenho para um efeito local, não é possível inferir os resultados para populações fora do entorno do ponto de corte.
Recomendações para a implementação de uma RDD:
Mostrar a distribuição de variáveis de confusão em função de \(X\)
Apresente os principais gr·Öcos de RD usando médias locais
Caso use uma especificação polinomial, faça o gráfico
Mostre sensibilidade dos resultados a diversas escolhas de \(h\) (bandwidth) na ordem do polinômio escolhido
Conduza um RD em paralelo para as covariadas
Mostre a sensibilidade dos resultados em relação a introdução de variáveis controle.
