library(wooldridge)
library(bacondecomp)
library(did)
library(haven)
library(ggplot2)
library(tidyverse)
library(modelsummary)
library(fixest)
library(lfe)Aula de laboratório - Diferenças em diferenças
Diferenças em diferenças (DID)
Modelos DID são atualmente os mais utilizados em estudos observacionais para investigar os efeitos de alguma política em uma variável de interesse. Este tipo de design faz uso das variações tanto entre grupo de indivíduos quanto no tempo para obter o efeito médio do tratamento sobre os tratados - ATT.
\[ \begin{align} DID &= \underbrace{\mathbb{E}[Y_{ist}|s = \text{NJ}, t = \text{Nov}] - \mathbb{E}[Y_{ist}|s = \text{NJ}, t = \text{Fev}]}_{\Delta_t^{GT}} - \Big(\underbrace{\mathbb{E}[Y_{ist}|s = \text{PA},t = \text{Nov}] - \mathbb{E}[Y_{ist}| s = \text{PA},t = \text{Fev}]}_{\Delta_t^{GC}}\Big) \\ \end{align} \]
onde estamos fazendo uso do exemplo Card e Krueger (1994). Pensilvânia (PA) é o estado de controle e Nova Jersey (NJ) o grupo tratado.
Em um formato de regressão, o estimador DID é obtido regredindo a variável de interesse em uma dummy para o grupo tratamento, uma outra dummy para o período após a intervenção e um termo de interação entre os dois, que denota o efetivo tratamento.
\[Y_{st} = \alpha + \beta TRAT_s + \gamma POS_t + \delta(TRAT_s \times POS_t) + \varepsilon_{st}\]
Nesta especificação \(\hat\delta = \text{DID}\) e a vantagem de fazer a regressão é que podemos fazer inferência sobre o parâmetro estimado diretamente.
Vamos fazer um exemplo encontrado no livro-texto do Wooldridge, exemplo 13.3 os efeitos de um incinerador de lixo no preço das moradias.
data("kielmc")
# DID manual
tratados <- kielmc[kielmc$nearinc == 1, ]
controles <- kielmc[kielmc$nearinc == 0, ]
delta_t <- mean(tratados[tratados$year == 1981, "rprice"]) -
mean(tratados[tratados$year == 1978, "rprice"])
delta_c <- mean(controles[controles$year == 1981, "rprice"]) -
mean(controles[controles$year == 1978, "rprice"])
did_manual <- delta_t - delta_cEstime o efeito com um DID em forma de regressão
Código
did_reg <- feols(price ~ nearinc * y81,
data = kielmc,
vcov = "HC1"
)
etable(did_reg) did_reg
Dependent Var.: price
Constant 82,517.2*** (1,882.4)
nearinc -18,824.4** (5,996.5)
y81 49,385.2*** (4,275.9)
nearinc x y81 -21,131.8* (10,070.7)
_______________ _____________________
S.E. type Heteroskedastic.-rob.
Observations 321
R2 0.35617
Adj. R2 0.35008
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Digamos que a hipótese de tendências paralelas seja válida somente após controlar pela idade do imóvel e idade ao quadrado, distância para uma rodovia, área do terreno e da moradia, número de quartos e banheiros. Gostaríamos também que o modelo fosse especificado em logarítimos, para ter uma ideia de variação percentual do preço da moradia em relação a presença do incinerador.
Calcule o efeito no cenário acima descrito
Código
did_ctr <- feols(log(rprice) ~ nearinc * y81 + sw0(age + agesq + lintst + log(land) + log(area) + rooms + baths),
data = kielmc,
vcov = "HC1"
)
etable(did_ctr) did_ctr.1 did_ctr.2
Dependent Var.: log(rprice) log(rprice)
Constant 11.29*** (0.0251) 7.652*** (0.5262)
nearinc -0.3399*** (0.0623) 0.0322 (0.0633)
y81 0.1931*** (0.0405) 0.1621*** (0.0272)
nearinc x y81 -0.0626 (0.0947) -0.1315* (0.0599)
age -0.0084*** (0.0016)
agesq 3.76e-5*** (1.05e-5)
lintst -0.0614. (0.0350)
log(land) 0.0998** (0.0328)
log(area) 0.3508*** (0.0621)
rooms 0.0473** (0.0171)
baths 0.0943*** (0.0281)
_______________ ___________________ ____________________
S.E. type Heteroskedast.-rob. Heteroskedasti.-rob.
Observations 321 321
R2 0.24597 0.73256
Adj. R2 0.23884 0.72394
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Replicando Cheng and Hoekstra (2013), quase 😄
Cheng e Hoekstra (2013) avaliaram o impacto que uma reforma de armas teve sobre a violência e para ilustrar vários princípios e práticas em relação ao timing diferencial. Vamos discutir esses princípios no contexto deste artigo.
Considerando que antes a autodefesa letal era legal apenas dentro de casa, uma nova lei, “Stand Your Ground”, estendeu esse direito a outros locais públicos. Entre 2000 e 2010, vinte e um estados expandiram explicitamente o estatuto da doutrina do castelo, estendendo os locais fora de casa onde a força letal poderia ser usada legalmente.
Cheng e Hoekstra (2013) escolheram um projeto de diferença em diferenças para seu projeto, onde a lei da doutrina do castelo era o tratamento e o tempo era diferenciado entre os estados. Sua equação de estimativa foi
\[Y_{it}=\alpha+\delta D_{it}+\gamma X_{it}+\sigma_i+\tau_t+\varepsilon_{it}\]
Cheng e Hoekstra (2013) não encontraram evidências de que as leis da doutrina do castelo dissuadissem ofensas violentas, mas descobriram que isso aumentava os homicídios. Um aumento líquido de 8% nas taxas de homicídios se traduz em cerca de seiscentos homicídios adicionais por ano nos vinte e um estados que o adotaram.
read_data <- function(df) {
full_path <- paste("https://raw.github.com/scunning1975/mixtape/master/",
df,
sep = ""
)
df <- read_dta(full_path)
return(df)
}
castle <- read_data("castle.dta")
#--- global variables
crime1 <- c(
"jhcitizen_c", "jhpolice_c",
"murder", "homicide",
"robbery", "assault", "burglary",
"larceny", "motor", "robbery_gun_r"
)
# demographics
demo <- c(
"emo", "blackm_15_24", "whitem_15_24",
"blackm_25_44", "whitem_25_44"
)
# variables dropped to prevent colinearity
dropped_vars <- c(
"r20004", "r20014",
"r20024", "r20034",
"r20044", "r20054",
"r20064", "r20074",
"r20084", "r20094",
"r20101", "r20102", "r20103",
"r20104", "trend_9", "trend_46",
"trend_49", "trend_50", "trend_51"
)
# Trend columns to keep
lintrend <- castle %>%
select(starts_with("trend")) %>%
colnames() %>%
# remove due to colinearity
subset(., !. %in% dropped_vars)
# Region columns to keep
region <- castle %>%
select(starts_with("r20")) %>%
colnames() %>%
# remove due to colinearity
subset(., !. %in% dropped_vars)
exocrime <- c("l_lacerny", "l_motor")
spending <- c("l_exp_subsidy", "l_exp_pubwelfare")
xvar <- c(
"blackm_15_24", "whitem_15_24", "blackm_25_44", "whitem_25_44",
"l_exp_subsidy", "l_exp_pubwelfare",
"l_police", "unemployrt", "poverty",
"l_income", "l_prisoner", "l_lagprisoner"
)
law <- c("cdl")Por exemplo, vamos conferir a evolução dos homicídios na Florida (sid = 10) contra a Califórnia (sid = 5)
fl_ca <- castle |>
filter(sid %in% c(5, 10)) |>
mutate(
sid = as.factor(sid),
year_fct = as.factor(year)
)
ggplot(fl_ca, aes(year_fct, l_homicide, group = sid, color = sid)) +
geom_line() +
scale_color_discrete(labels = c("California", "Florida")) +
geom_vline(xintercept = "2005", linetype = 2) +
labs(
x = "Ano",
y = "Log homicídios por 100K hab.",
title = "Log homicídios FL vs CA"
) +
guides(color = guide_legend(title = "Estado")) +
theme_classic()
Vamos conferir também quais estados passaram este tipo de lei e quando
cdl_sid <- castle |>
group_by(sid, state) |>
summarise(sum_post = sum(post)) |>
filter(sum_post > 0) |>
pull(sid)`summarise()` has grouped output by 'sid'. You can override using the `.groups`
argument.cdl_states <- castle |>
select(sid, state, year, post) |>
filter(sid %in% cdl_sid, post == 0) |>
group_by(sid) |>
slice(n()) |>
arrange(year)
kableExtra::kbl(cdl_states[, c("state", "year")],
col.names = c("Estado", "Ano CDL")
) |>
kableExtra::kable_classic(full_width = FALSE)| Estado | Ano CDL |
|---|---|
| Florida | 2005 |
| Alabama | 2006 |
| Alaska | 2006 |
| Arizona | 2006 |
| Georgia | 2006 |
| Indiana | 2006 |
| Kansas | 2006 |
| Kentucky | 2006 |
| Louisiana | 2006 |
| Michigan | 2006 |
| Mississippi | 2006 |
| Oklahoma | 2006 |
| South Carolina | 2006 |
| South Dakota | 2006 |
| Missouri | 2007 |
| North Dakota | 2007 |
| Tennessee | 2007 |
| Texas | 2007 |
| Ohio | 2008 |
| West Virginia | 2008 |
| Montana | 2009 |
Este tipo de design é o que chama-se de adoção escalonada ao tratamento (staggered adoption). Existem diversas coortes que aderem ao tratamento em momentos distintos
dd_formula <- as.formula(
paste(
"l_homicide ~ ",
paste(
paste(xvar, collapse = " + "),
paste(region, collapse = " + "),
paste(lintrend, collapse = " + "),
paste("post", collapse = " + "),
sep = " + "
),
"| year + sid | 0 | sid"
)
)
# Fixed effect regression using post as treatment variable
dd_reg <- felm(dd_formula, weights = castle$popwt, data = castle)
msummary(dd_reg,
output = "kableExtra",
fmt = 4,
coef_map = "post",
gof_map = c("nobs", "adj.r.squared", "vcov.type")
) |>
kableExtra::kable_classic()| (1) | |
|---|---|
| post | 0.0769 |
| (0.0336) | |
| Num.Obs. | 550 |
| R2 Adj. | 0.951 |
| Std.Errors | by: sid |
Agora, vamos ir além do estudo deles e implementar um estudo de evento. Para fazer isso, usamos uma variável “time_til”, que é o número de anos até ou depois que o estado recebeu o tratamento. Usando essa variável, criamos os leads (que serão os anos anteriores ao tratamento) e os lags (os anos pós-tratamento). Felizmente o dataset providenciado já possui estas dummies.
event_study_formula <- as.formula(
paste(
"l_homicide ~ + ",
paste(
paste(region, collapse = " + "),
paste(paste("lead", 1:9, sep = ""), collapse = " + "),
paste(paste("lag", 1:5, sep = ""), collapse = " + "),
sep = " + "
),
"| year + state | 0 | sid"
),
)
event_study_reg <- felm(event_study_formula, weights = castle$popwt, data = castle)
msummary(event_study_reg,
output = "kableExtra",
fmt = 4,
coef_omit = "^(?!lead|lag)", # lookahed regex. starts with lead or lag
gof_map = c("nobs", "adj.r.squared", "vcov.type")
) |>
kableExtra::kable_classic()| (1) | |
|---|---|
| lead1 | -0.0255 |
| (0.0327) | |
| lead2 | 0.0191 |
| (0.0308) | |
| lead3 | 0.0122 |
| (0.0346) | |
| lead4 | -0.0041 |
| (0.0461) | |
| lead5 | 0.0050 |
| (0.0466) | |
| lead6 | 0.0093 |
| (0.0590) | |
| lead7 | -0.1372 |
| (0.0854) | |
| lead8 | -0.3036 |
| (0.0808) | |
| lead9 | -0.2607 |
| (0.0446) | |
| lag1 | 0.0777 |
| (0.0279) | |
| lag2 | 0.0824 |
| (0.0449) | |
| lag3 | 0.1050 |
| (0.0524) | |
| lag4 | 0.0787 |
| (0.0583) | |
| lag5 | 0.1724 |
| (0.0554) | |
| Num.Obs. | 550 |
| R2 Adj. | 0.932 |
| Std.Errors | by: sid |
É comum plotar esses estudos de eventos. Utiliza-se gráficos tipo point-range, onde a estimativa pontual e o intervalo de confiança são mostrados. Vamos fazer isso agora.
# order of the coefficients for the plot
plot_order <- c(
"lead9", "lead8", "lead7",
"lead6", "lead5", "lead4", "lead3",
"lead2", "lead1", "lag1",
"lag2", "lag3", "lag4", "lag5"
)
# grab the clustered standard errors
# and average coefficient estimates
# from the regression, label them accordingly
# add a zero'th lag for plotting purposes
leadslags_plot <- tibble(
sd = c(event_study_reg$cse[plot_order], 0),
mean = c(coef(event_study_reg)[plot_order], 0),
label = c(-9:-1, 1:5, 0)
)
# This version has a point-range at each
# estimated lead or lag
# comes down to stylistic preference at the
# end of the day!
leadslags_plot %>%
ggplot(aes(
x = label, y = mean,
ymin = mean - 1.96 * sd,
ymax = mean + 1.96 * sd
)) +
geom_hline(
yintercept = coef(dd_reg)["post"],
color = "red"
) +
geom_pointrange() +
theme_classic() +
xlab("Anos antes e depois da expansão da lei doutrina do castelo") +
ylab("Log homicídios por 100K hab.") +
geom_hline(
yintercept = 0,
linetype = "dashed"
) +
geom_vline(
xintercept = 0,
linetype = "dashed"
)
# This version includes
# an interval that traces the confidence intervals
# of your coefficients
leadslags_plot %>%
ggplot(aes(
x = label, y = mean,
ymin = mean - 1.96 * sd,
ymax = mean + 1.96 * sd
)) +
# this creates a red horizontal line
geom_hline(
yintercept = coef(dd_reg)["post"],
color = "red"
) +
geom_line() +
geom_point() +
geom_ribbon(alpha = 0.2) +
theme_minimal() +
# Important to have informative axes labels!
xlab("Anos antes e depois da expansão da lei doutrina do castelo") +
ylab("Log homicídios por 100K hab.") +
geom_hline(yintercept = 0) +
geom_vline(xintercept = 0)
Uma literatura recente tem demonstrado que esta especificação pode ser viesda para o efeito do tratamento sobre os tratados, ou seja, TWFE com adocação escalonada pode não ser equivalente ao estimador de DID. Vamos utilizar a mesma especificação, porém utilizando o estimador de Callaway & Sant’anna encontrado no pacote did.
# Precisamos adaptar o data.frame para a funcao att_gt
castle_cs <- replace_na(castle, list(effyear = 0)) |>
select(sid, year, effyear, popwt, l_homicide, r20001)
did_gt <- att_gt(
yname = "l_homicide",
tname = "year",
idname = "sid",
gname = "effyear",
weightsname = "popwt",
# xformla = ~r20001,
control_group = "notyettreated",
data = castle_cs
)
did_es <- aggte(did_gt, type = "dynamic")
ggdid(did_es) +
scale_x_continuous(n.breaks = 10) +
geom_hline(
yintercept = coef(dd_reg)["post"],
color = "red"
) +
labs(
x = "Anos antes e depois da expansão da lei doutrina do castelo",
y = "Log homicídios por 100K hab."
)
Vejam que as estimativas mudaram ligeiramente apenas. Muito embora a especificação TWFE deva ser viesada em teoria, este viés aparentemente é pequeno. Podemos averiguar esta situação através da decomposição de Goodman-Bacon. Se nenhum dos pesos atribuídos as estimativas for negativo, então o viés do TWFE tende a ser pequeno (pesos negativos exacerbam o viés além de prejudicar uma interpretação causal do estimador TWFE).
bacon_d <- bacon(l_homicide ~ post,
data = castle,
id_var = "sid",
time_var = "year"
) type weight avg_est
1 Earlier vs Later Treated 0.07708 -0.02858
2 Later vs Earlier Treated 0.02411 0.04563
3 Treated vs Untreated 0.89881 0.07844