Diff_in_Diff-exemplos

Author

Igor Viveiros e João Cardoso

Published

August 10, 2023

1 Objetivo

Esta aula será focada em discutir exemplos da aplicação da DiD em dados reais. Ao todo, discutiremos três exemplos: Card e Krueger (1994), Fraser et al (2020).
Serão utilizados os seguintes pacotes:

Code

library(dplyr)
library(readr)
library(stringr)
library(ggplot2)
library(tidyr)
library(sjlabelled)
library(ggrepel)
library(scales)
library(ggpubr)
library(plm)
library(lmtest)
library(sjPlot)
library(DT)
library(magrittr)
library(foreign)
library(dataverse)
library(Synth) 
# Pacotes especificos para a última análise
#install.packages("remotes")
#remotes::install_github("kylebutts/did2s")
library(did2s)

2 Exemplo 1: Replicação de Card e Krueger (1994)

O estudo tem como objetivo identificar o efeito do aumento do salário-mínimo nos empregos. A hipótese levantada pela literatura é que um aumento no salário-mínimo leva a um aumento do desemprego. A partir de um desenho quase-experimental com uso de DiD, Card e Krueger (1994) demonstraram o contrário, que o aumento do salário-mínimo levou a uma redução do desemprego. Em 1992, o estado americano de Nova Jersey (NJ) aumentou o salário-mínimo de $4,25 para $5,05 por hora, abrindo espaço para testar a teoria. Como controle, os pesquisadores utilizaram o estado vizinho da Pennsylvania (PA), em que não houve alteração na política. Foram conduzidos surveys em ambos os estados antes e depois da aplicação da política em NJ, sendo representativo dos restaurantes de fast-food. Vale destacar que esse desenho é um quase-experimento pois as unidades não são iguais e o tratamento não foi aplicado de forma aleatória.
O primeiro passo para a análise é limpar e transformar os dados. O banco de dados foi baixado diretamente do site dos autores e transformado em um formato necessário para realizar as análises em DiD. O formato final do banco é estruturado da seguinte forma:

Code

# Criar arquivo e diretório temporário
tfile_path <- tempfile()
tdir_path <- tempdir()

# Download zip file
download.file("http://davidcard.berkeley.edu/data_sets/njmin.zip", 
              destfile = tfile_path)

# Unzip
unzip(tfile_path, exdir = tdir_path)

# Ler o codebook
codebook <- read_lines(file = paste0(tdir_path, "/codebook"))

# Criando um vetor com o nome das variáveis
variable_names <- codebook %>%
  `[`(8:59) %>% # Nome das variáveis a partir do elemento 8 (planilha)
  `[`(-c(5, 6, 13, 14, 32, 33)) %>% # Removendo elementos sem nome de variáveis
  str_sub(1, 8) %>% # Selecionando apenas as variáveis com 8 elementos
  str_squish() %>% # Removendo espaço em branco
  str_to_lower() # Transformando todos os nomes em letras minusculas

# Criando um vetor com o nome das variáveis
variable_labels <- codebook %>%
  `[`(8:59) %>% # Nome das variáveis a partir do elemento 8 (planilha)
  `[`(-c(5, 6, 13, 14, 32, 33)) %>% # Removendo elementos sem nome de variáveis
  sub(".*\\.[0-9]", "", .) %>%
  `[`(-c(5:10))  %>% # Elementos que serão combinados
  str_squish() # Removendo espaço em branco
  
# Região
variable_labels[41] <- "region of restaurant"

# Lendo os dados "cru"
data_raw <- read_table2(paste0(tdir_path, "/public.dat"),
                        col_names = FALSE)

# Adicionando o nome das variáveis
data_mod <- data_raw %>%
  select(-X47) %>% # Removendo colunas vazias
  `colnames<-`(., variable_names) %>% # Nomeando as variáveis
  mutate_all(as.numeric) %>% # Tratar todas as variáveis como númericas
  mutate(sheet = ifelse(sheet == 407 & chain == 4, 408, sheet)) # Duplicando o id 407

# Dados processados em formato wide
data_mod <- data_mod %>%
  # Label dos valores de chain 
  mutate(chain = case_when(chain == 1 ~ "bk",
                           chain == 2 ~ "kfc",
                           chain == 3 ~ "roys",
                           chain == 4 ~ "wendys")) %>%
  # Label dos valores de state
  mutate(state = case_when(state == 1 ~ "New Jersey",
                           state == 0 ~ "Pennsylvania")) %>%
  # Dummy de região
  mutate(region = case_when(southj == 1 ~ "southj",
                            centralj == 1 ~ "centralj",
                            northj == 1 ~ "northj",
                            shore == 1 ~ "shorej",
                            pa1 == 1 ~ "phillypa",
                            pa2 == 1 ~ "eastonpa")) %>%
  # Label dos valores de meals
  mutate(meals = case_when(meals == 0 ~ "none",
                           meals == 1 ~ "free meals",
                           meals == 2 ~ "reduced price meals",
                           meals == 3 ~ "both free and reduced price meals")) %>%
  # Label dos valores de meals2
  mutate(meals2 = case_when(meals2 == 0 ~ "none",
                            meals2 == 1 ~ "free meals",
                            meals2 == 2 ~ "reduced price meals",
                            meals2 == 3 ~ "both free and reduced price meals")) %>%
  # Label dos valores de status2 
  mutate(status2 = case_when(status2 == 0 ~ "refused second interview",
                             status2 == 1 ~ "answered 2nd interview",
                             status2 == 2 ~ "closed for renovations",
                             status2 == 3 ~ "closed permanently",
                             status2 == 4 ~ "closed for highway construction",
                             status2 == 5 ~ "closed due to Mall fire")) %>%
  mutate(co_owned = if_else(co_owned == 1, "yes", "no")) %>%
  mutate(bonus = if_else(bonus == 1, "yes", "no")) %>%
  mutate(special2 = if_else(special2 == 1, "yes", "no")) %>%
  mutate(type2 = if_else(type2 == 1, "phone", "personal")) %>%
  select(-southj, -centralj, -northj, -shore, -pa1, -pa2) %>% # Incluindo dummies de região
  mutate(date2 = lubridate::mdy(date2)) %>% # Convertendo a data
  rename(open2 = open2r) %>% # Inserindo nome da onda 1
  rename(firstinc2 = firstin2) %>% # Inserindo nome da onda 1
  sjlabelled::set_label(variable_labels) # Adicionando as labels guardadas

# Variáveis estruturais
structure <- data_mod %>%
  select(sheet, chain, co_owned, state, region)

# Variáveis da onda 1
wave1 <- data_mod %>%
  select(-ends_with("2"), - names(structure)) %>%
  mutate(observation = "February 1992") %>%
  bind_cols(structure) 

# Variáveis da onda 2
wave2 <- data_mod %>%
  select(ends_with("2")) %>%
  rename_all(~str_remove(., "2"))  %>%
  mutate(observation = "November 1992") %>%
  bind_cols(structure) 

# Banco de dados final
card_krueger_1994 <- bind_rows(wave1, wave2) %>%
  select(sort(names(.))) %>% # Ordenando as colunas por ordem alfabética
  sjlabelled::copy_labels(data_mod) # Restaurando a label das variáveis

card_krueger_1994_mod <- card_krueger_1994 %>%
  mutate(emptot = empft + nmgrs + 0.5 * emppt,
         pct_fte = empft / emptot * 100) |> 
  select(c(observation, state, emptot, everything()))

set.seed(12361)
sample_n(card_krueger_1994_mod, 5) |> tab_df(title = "Estrutura do banco")

Estrutura do banco
observation	state	emptot	bonus	chain	co_owned	date	empft	emppt	firstinc	hrsopen	inctime	meals	ncalls	nmgrs	nregs	nregs11	open	pctaff	pentree	pfry	psoda	region	sheet	special	status	type	wage_st	pct_fte
November 1992	New Jersey	22.50	NA	wendys	no	1992-11-17	10	17	0.13	12.50	14	reduced price meals	1	4	2	2	10.50	NA	1.05	0.98	1.05	northj	165	no	answered 2nd interview	phone	5.05	44.44
February 1992	Pennsylvania	20.50	yes	roys	yes	NA	10	15	0.17	16.00	13	reduced price meals	0	3	5	3	6.00	0	0.73	0.90	1.05	eastonpa	496	NA	NA	NA	5.00	48.78
February 1992	New Jersey	24.00	no	roys	no	NA	10	20	NA	12.50	26	none	0	4	5	3	9.00	NA	NA	NA	NA	centralj	135	NA	NA	NA	NA	41.67
February 1992	New Jersey	25.50	no	roys	yes	NA	15	15	NA	17.00	17	reduced price meals	0	3	4	4	6.00	NA	1.06	NA	1.02	northj	306	NA	NA	NA	NA	58.82
February 1992	New Jersey	6.50	no	kfc	no	NA	0	7	0.75	12.00	NA	both free and reduced price meals	1	3	3	2	11.00	40	2.45	1.06	1.06	northj	73	NA	NA	NA	4.25	0.00

Nesse caso, nossa variável de tempo é $observation$, nossa variável que especifica o estado é $state$ e nossa variável dependente é $emptot$. Para argumentar a existencia de tendencias paralelas, os autores utilizam a distribuição de salários dentro de fast foods pelos estados. Os histogramas a seguir mostram o percentual de restaurantes por valor do valor de sálario em ambos os estados estudados:

Code

hist.feb <- card_krueger_1994_mod %>%
  filter(observation == "February 1992") %>%
  ggplot(aes(wage_st, fill = state)) +
  geom_histogram(aes(y=c(..count..[..group..==1]/sum(..count..[..group..==1]),
                         ..count..[..group..==2]/sum(..count..[..group..==2]))*100),
                 alpha=0.5, position = "dodge", bins = 23) +
  labs(title = "Fevereiro 1992", x="Distribuição do salário", y = "Percentual de restaurantes", fill = "") +
  scale_fill_grey()

hist.nov <- card_krueger_1994_mod %>%
  filter(observation == "November 1992") %>%
  ggplot(aes(wage_st, fill = state)) +
  geom_histogram(aes(y=c(..count..[..group..==1]/sum(..count..[..group..==1]),
                         ..count..[..group..==2]/sum(..count..[..group..==2]))*100),
                 alpha = 0.5, position = "dodge", bins = 23) +
  labs(title = "Novembro 1992", x="Distribuição do salário", y = "Percentual de restaurantes", fill="") +
  scale_fill_grey()

ggarrange(hist.feb, hist.nov, ncol = 2, 
          common.legend = TRUE, legend = "bottom")

Antes de usar uma abordagem de DiD com modelos de regressão OLS, podemos tentar calcular a diferença das diferenças nas médias a partir de uma abordagem de First Differences (FD). Assim, o ATT (Avarage Treatment Effect) é obtido subtraindo as médias do grupo tratamento antes e depois da intervenção pelas médias do grupo controle antes e depois da intervenção.

\[ATT = (\mu NJ_{t_0} - \mu NJ_{t_1}) - (\mu PA_{t_0} - \mu PA_{t_1})\]

Code

# Calculando as médias pelo grupo
differences <- card_krueger_1994_mod %>%
  group_by(observation, state) %>%
  summarise(emptot = mean(emptot, na.rm = TRUE))

# Grupo tratado (NJ) antes do tratamento
njfeb <- differences[1,3]

# Grupo controle (PA) antes do tratamento
pafeb <- differences[2,3]

# Grupo tratado (NJ) depois do tratamento
njnov <- differences[3,3]

# Grupo controle (PA) depois do tratamento
panov <- differences[4,3]

# Calculando o ATT
ATT <- (njnov-njfeb)-(panov-pafeb)

print(paste("ATT = ", round(ATT, 3)))

[1] "ATT =  2.754"

Assim, o efeito causal apontado pelo calculo de DiD a partir de FD é de um aumento de 2,754 pontos no emprego. Plotando visualmente:

Code

# Calculando o outcome do contrafactual
nj_counterfactual <- tibble(
  observation = c("February 1992","November 1992"), 
  state = c("New Jersey (Counterfactual)","New Jersey (Counterfactual)"),
  emptot = as.numeric(c(njfeb, njfeb-(pafeb-panov)))
  ) 

# Pontos para a ocorrência da intevenção
intervention <- tibble(
    observation = c("Intervention", "Intervention", "Intervention"),
    state = c("New Jersey", "Pennsylvania", "New Jersey (Counterfactual)"),
    emptot = c(19.35, 22.3, 19.35)
  ) 

# Unindo os bancos
did_plotdata <- bind_rows(differences, 
                          nj_counterfactual, 
                          intervention)

# Plot
did_plotdata %>%
  mutate(label = if_else(observation == "November 1992", as.character(state), NA_character_)) %>%
  ggplot(aes(x=observation,y=emptot, group=state)) +
  geom_line(aes(color=state), size=1.2) +
  geom_vline(xintercept = "Intervention", linetype="dotted", 
             color = "black", size=1.1) + 
  scale_color_brewer(palette = "Accent") +
  scale_y_continuous(limits = c(17,24)) +
  ggrepel::geom_label_repel(aes(label = label),
                   nudge_x = 0.5, nudge_y = -0.5,
                   na.rm = TRUE) +
  guides(color=FALSE) +
  labs(x="", y="FTE Employment (mean)") +
  annotate(
    "text",
    x = "November 1992",
    y = 19.6,
    label = "{Difference-in-Differences}",
    angle = 90,
    size = 3
  )+
  theme_bw()

Agora, estimaremos este modelo por meio de uma regressão OLS simples: \[Y = \beta_0 + \beta_1 Time + \beta_2 Treated + \beta (Time \times Treated) +\epsilon\]

Code

# Criando as dummies de tempo e tratamento
card_krueger_1994_mod <- mutate(card_krueger_1994_mod,
                                time = ifelse(observation == "November 1992", 1, 0),
                                treated = ifelse(state == "New Jersey", 1, 0)
                                )

# Estimando o modelo OLS
did_model <- lm(emptot ~ time + treated + time:treated, data = card_krueger_1994_mod)

did_model|>  tab_model()

	#full-time employees
Predictors	Estimates	CI	p
(Intercept)	23.33	21.23 – 25.44	<0.001
time	-2.17	-5.14 – 0.81	0.154
treated	-2.89	-5.23 – -0.55	0.016
time × treated	2.75	-0.56 – 6.07	0.103
Observations	794
R² / R² adjusted	0.007 / 0.004

O modelo aponta que o ATT do aumento do salário mínimo é de um aumento de 2,75, contudo, o resultado não é estatisticamente significante. Contudo, devemos avaliar se os dados estão corretamente estimados por meio de testes, a fim de avaliar se o modelo deve incorporar efeitos fixos ou aleatórios e o uso de erros robustos devido a heterosedasticidade. Primeiramente, vamos testar se devemos usar um modelo de Efeitos Fixos por meio de um teste F:

Code

# Defindo os dados em painel
panel <- pdata.frame(card_krueger_1994_mod, "sheet")

# Construindo o modelo within (Efeitos Fixos)
EF <- plm(emptot ~ time + treated + time:treated, 
               data = panel, model = "within")

# Construindo o modelo pooled (Equivalente a OLS)
PO <- plm(emptot ~ time + treated + time:treated, 
               data = panel,
          model = "pooling")

# Teste F, lembrar que EF vem sempre primeiro
pFtest(EF,PO)


    F test for individual effects

data:  emptot ~ time + treated + time:treated
F = 3.3237, df1 = 408, df2 = 382, p-value < 2.2e-16
alternative hypothesis: significant effects

O teste aponta para a existencia de efeitos específicos, $\alpha_i$, logo, devemos utilizar um modelo de efeitos fixos ou aleatórios. Para definir qual dos dois modelos, realizaremos o teste de Hausman

Code

# Construindo o modelo random (Efeitos Aleatórios)
EA <- plm(emptot ~ time + treated + time:treated, 
               data = panel,
          model = "random")

# Teste de Hausman, lembrar que EF vem sempre primeiro
phtest(EF,EA)


    Hausman Test

data:  emptot ~ time + treated + time:treated
chisq = 0.63463, df = 2, p-value = 0.7281
alternative hypothesis: one model is inconsistent

O resultado aponta para o uso de Efeitos Fixos no lugar de Efeitos Aleatórios. Não precisamos testar efeitos de tempo, pois o painel é composto por apenas dois pontos no tempo que são capturados pelo tempo. Assim, estimamos o modelo com Efeitos Fixos:

Code

EF |> tab_model()

	emptot
Predictors	Estimates	CI	p
time [2]	-2.28	-4.32 – -0.25	0.028
time [2] × treated	2.75	0.48 – 5.02	0.018
Observations	794
R² / R² adjusted	0.015 / -1.045

Observamos que houve uma redução dos erros padrões, que agora apontam significancia estatistica para o resultado. Agora avaliaremos se existe heterocedasticidade nos resultados. Podemos fazer isso pela análise gráfica:

Code

plot <- tibble(fitted(EF),resid(EF))

ggplot(data = plot, aes(x = `fitted(EF)`, y = `resid(EF)`^2)) + 
         geom_point() +
         labs(title = "Resíduos ao quadrado") +
         geom_smooth(method = "loess", se = FALSE)+
  theme_bw()

A análise visual aponta para a existencia de heterocedasticidade nos dados. Vamos conduzir um teste de Breusch-Pagan para avaliar se isso se confirma:

Code

bptest(EF, studentize=F)


    Breusch-Pagan test

data:  EF
BP = 15.526, df = 3, p-value = 0.001418

Confirmando a análise gráfica, o teste aponta para a existência de heterocedasticidade dos dados. Assim, aplicaremos uma matrix de correção de Arellano para obter os erros padrões robustos.

Code

# Obtendo os erros padrões robustos
ARL  <- vcovHC(EF, method = "arellano")

EF.ARL <- EF; EF.ARL$vcov <- ARL

tab_model(EF, EF.ARL,
          dv.labels = c("Modelo sem erros padrões robustos","Modelo com erros padrões robustos"), 
          show.se = T,
          show.ci = F,
          digits = 4
)

	Modelo sem erros padrões robustos			Modelo com erros padrões robustos
Predictors	Estimates	std. Error	p	Estimates	std. Error	p
time [2]	-2.2833	1.0355	0.028	-2.2833	1.2449	0.067
time [2] × treated	2.7500	1.1544	0.018	2.7500	1.3342	0.040
Observations	794			794
R² / R² adjusted	0.015 / -1.045			0.015 / -1.045

Assim, estimamos o efeito do tratamento com os erros padrões corretos, indicando significancia estatística a um p-valor menor que 0,05. Nesse modelo, contudo não é possível testar a existência de tendências paralelas por termos apenas dois períodos no tempo.

3 Exemplo 2: Replicação de Fraser et al (2020)

O segundo exemplo foca em como o processo de organização social impacta a emissão de gás carbono. A teoria aponta que o aumento da organização social levaria a uma redução na emissão de gás carbono. Assim, o estudo mobilizou dados de diferentes cidades na região de Tóquio, avaliando sua capacidade de organização social e a classificando em alta (1) e baixa (0). Assim, utilizaremos a DiD para estimar o efeito ao longo do tempo. A base de dados a ser utilizada segue a seguinte estrutura:

Code

#Limpando os arquivos do exemplo anterior
rm(list=ls())

# Carregando os dados
Sys.setenv("DATAVERSE_SERVER" = "dataverse.harvard.edu")
  
# Importando o banco de dados do Havard Dataverse
dat <- get_dataframe_by_name(
    filename  = "dataset.tab",
    dataset   = "10.7910/DVN/E8PEW7")

# Limpando as cidades
dat <- dat %>%
  # Selecionando apenas as cidades na região de Tokyo
  filter(pref == "Tokyo") %>%
  # Calculando a emissão de carbono em kilotons por 1000 residentes
  mutate(emissions = emissions / pop * 1000) %>%
  # Calculando a densidade populacional
  mutate(pop_density = pop / area_inhabitable) %>%
  # Criando uma contagem anual de 1, 2, 3, 4, etc.
  mutate(counter = year - min(year)) %>%
  # Criando um tratamento se a cidade é ou foi uma cidade com alta capacidade e a contando como tratada pelo resto do periodo
  mutate(treat = if_else(fin_str_index > median(fin_str_index), 1, NA_real_)) %>%
  group_by(muni_code) %>%
  fill(treat, .direction = "up") %>%
  mutate(treat = if_else(is.na(treat), 0, treat)) %>%
  ungroup() %>%
  # Garantidno um fator para efeitos fixos
  mutate(year = factor(year),
         muni_code = factor(muni_code)) %>%
  # Selecionando apenas as variáveis de interesse:
   select(
    muni, # nome do município
    muni_code, # código único de cinco dígitos para identificação de cada cidade
    year, # ano
    counter, # número de anos passados desde 2005 (0, 1, 2, ...),
    emissions, # emissões em quilotons por 1.000 habitantes
    treat, # Alta capacidade de organização (1/0)
    fin_str_index, # capacidade governamental (numérica) e um conjunto de covariáveis que também podem influenciar as emissões, incluindo:
    pop_density, # densidade populacional
    income_per_capita, # renda per capita
    pop_age_65_plus, # % de pessoas com 65 anos ou mais
    pop_college, # % com ensino superior completo
    value_agr_mill, # produção agrícola em milhões de ienes por habitante
    value_manuf_mill, # produção industrial em milhões de ienes por habitante
    value_commerce_mill # produção comercial em milhões de ienes por habitante
)


set.seed(12361)
sample_n(dat, 5) |> tab_df(title = "Estrutura do banco")

Estrutura do banco
muni	muni_code	year	counter	emissions	treat	fin_str_index	pop_density	income_per_capita	pop_age_65_plus	pop_college	value_agr_mill	value_manuf_mill	value_commerce_mill
Tokyo-to Tama-shi	13224	2005	0	3.86	1	1.22	74.01	1758.41	15.75	22.97	0.00	0.15	1.64
Tokyo-to Chuo-ku	13102	2005	0	33.35	0	0.52	96.94	2985.35	16.26	38.18	0.21	1.03	429.50
Tokyo-to Sumida-ku	13107	2008	3	5.53	1	0.56	180.08	1785.14	21.31	16.60	0.22	1.46	7.13
Tokyo-to Koto-ku	13108	2010	5	5.29	1	0.53	115.38	1941.81	19.11	21.25	0.22	0.66	10.17
Tokyo-to Arakawa-ku	13118	2017	12	4.09	1	0.50	208.92	1796.20	23.05	15.34	0.34	0.43	2.79

Sendo $muni$ o nome do município, $municode$ o código do múnicipio, $year$ o ano da observação, $counter$ a passagem de ano a partir de 2005, que servirá como variável de exposição e indicará o aumento do efeito ao longo do tempo, $emissions$ é a variável dependente e $treat$ é a variável que indica se o grupo tem ou não uma alta organização social. O restante das variáveis são controles que a literatura aponta como importantes para a emissão de gás carbono pela população. O gráfico a seguir

Code

# Identificação de tendencias paralelas

dat |> 
  group_by(year,treat) |> 
  summarise(media = mean(emissions)) |> 
  mutate(treat = as.factor(treat)) |> 
  ggplot(aes(x = year, y = media, group = treat, color = treat))+
  geom_path()+
  labs(title = "Variação da emissão ao longo do tempo", y = "Emissão de carbono em quilotons por 1000 habitantes", x = "Ano", color = "Tratamento")+
  theme_bw()

Primeiramente, estimaremos o modelo naive (sem controles), a fim de ver como a variável dependente pode ser afetada pelo modelo. Diferente do modelo anteriror, essa especificação conta que o efeito aumenta com o passar do tempo de forma regular. A especificação do modelo pode ser dada por:

\[Y = \beta_0 + \beta_1 Counter + \beta_2 Treat + \beta (Counter \times Treat)+\epsilon\]

Code

naive <- lm(emissions ~ treat * counter, 
            data = dat)

naive |> tab_model()

	emissions
Predictors	Estimates	CI	p
(Intercept)	11.77	8.81 – 14.73	<0.001
treat	-5.70	-9.46 – -1.93	0.003
counter	0.70	0.24 – 1.17	0.003
treat × counter	-0.53	-1.09 – 0.02	0.061
Observations	744
R² / R² adjusted	0.083 / 0.079

Sem controles, o modelo indica que não existem efeitos da organização social sobre a emissão de gás carbonico. Porém, veremos se esse resultado se mantem após a inserção dos controles apontados pela literatura. O novo modelo possui as seguintes especificações: \[Y = \beta_0 + \beta_1 Counter + \beta_2 Treat + \beta (Counter \times Treat)+ \beta_3 PopDensity+ \beta_4 IncomePerCapita+ \beta_5 PopAge65Plus+ \beta_7 PopCollege + \beta_7 ValueAgrMill + \beta_8 CommerceMill + \beta_9 ValueManufMill +\epsilon\]

Code

PO <- dat %>%
  plm(formula = emissions ~ I(treat * counter) + treat + counter +
        pop_density + income_per_capita + pop_age_65_plus + pop_college +
        value_agr_mill + value_commerce_mill + value_manuf_mill, 
      model = "pooling", index = c("muni_code", "year"))

PO |> tab_model()

	Dependent variable
Predictors	Estimates	CI	p
(Intercept)	13.94	9.39 – 18.49	<0.001
treat * counter	-0.88	-1.21 – -0.56	<0.001
treat	-3.46	-5.60 – -1.32	0.002
counter	1.09	0.80 – 1.37	<0.001
pop density	-0.07	-0.08 – -0.05	<0.001
income per capita	-0.00	-0.00 – -0.00	0.018
pop age 65 plus	-0.36	-0.46 – -0.25	<0.001
pop college	0.39	0.24 – 0.55	<0.001
value agr mill	18.89	12.43 – 25.34	<0.001
value commerce mill	0.08	0.07 – 0.08	<0.001
value manuf mill	0.72	0.47 – 0.98	<0.001
Observations	744
R² / R² adjusted	0.741 / 0.737

O modelo aponta que o ATT da organização social é de uma redução da emissão de 0,88 quilotons de gás carbono por ano, sendo esse resultado estatisticamente significante. Contudo, devemos avaliar se os dados estão corretamente estimados por meio de testes, a fim de avaliar se o modelo deve incorporar efeitos fixos ou aleatórios e o uso de erros robustos devido a heterosedasticidade. Primeiramente, vamos testar se devemos usar um modelo de Efeitos Fixos por meio de um teste F:

Code

FE <- dat %>%
  plm(formula = emissions ~ I(treat * counter) + treat + counter +
        pop_density + income_per_capita + pop_age_65_plus + pop_college +
        value_agr_mill + value_commerce_mill + value_manuf_mill, 
      model = "within", index = c("muni_code", "year"))

pFtest(FE, PO)


    F test for individual effects

data:  emissions ~ I(treat * counter) + treat + counter + pop_density +  ...
F = 84.889, df1 = 61, df2 = 672, p-value < 2.2e-16
alternative hypothesis: significant effects

O teste aponta para a existencia de efeitos específicos, $\alpha_i$, logo, devemos utilizar um modelo de efeitos fixos ou aleatórios. Para definir qual dos dois modelos, realizaremos o teste de Hausman

Code

RE <- dat %>%
  plm(formula = emissions~ I(treat * counter) + treat + counter +
        pop_density + income_per_capita + pop_age_65_plus + pop_college +
        value_agr_mill + value_commerce_mill + value_manuf_mill, 
      model = "random", index = c("muni_code", "year"))

phtest(FE, RE)


    Hausman Test

data:  emissions ~ I(treat * counter) + treat + counter + pop_density +  ...
chisq = 35.357, df = 10, p-value = 0.0001085
alternative hypothesis: one model is inconsistent

O resultado aponta para o uso de Efeitos Fixos no lugar de Efeitos Aleatórios. Não precisamos testar efeitos de tempo, pois o painel é composto por apenas dois pontos no tempo que são capturados pelo tempo.

Code

TWFE <- dat %>%
  plm(formula = emissions ~ I(treat * counter) + treat + counter +
        pop_density + income_per_capita + pop_age_65_plus + pop_college +
        value_agr_mill + value_commerce_mill + value_manuf_mill, 
      model = "within",effect ="twoways" , index = c("muni_code", "year"))

pFtest(TWFE, FE)


    F test for twoways effects

data:  emissions ~ I(treat * counter) + treat + counter + pop_density +  ...
F = 2.3142, df1 = 10, df2 = 662, p-value = 0.01118
alternative hypothesis: significant effects

O teste dá suporte para efeitos específicos de tempo. Assim, estimamos o modelo com Efeitos Fixos em Twoways:

Code

TWFE |> tab_model()

	Dependent variable
Predictors	Estimates	CI	p
treat * counter	-0.48	-0.62 – -0.35	<0.001
treat	0.70	-0.48 – 1.88	0.244
pop density	-0.09	-0.18 – 0.01	0.068
income per capita	0.00	0.00 – 0.00	0.007
pop age 65 plus	-0.06	-0.31 – 0.19	0.655
pop college	-0.60	-1.02 – -0.17	0.006
value agr mill	1.63	-7.57 – 10.84	0.728
value commerce mill	0.04	0.02 – 0.05	<0.001
value manuf mill	0.52	0.20 – 0.83	0.001
Observations	744
R² / R² adjusted	0.104 / -0.005

O modelo apontou para um efeito significativo do tratamento ao longo do tempo, indicando uma redução de 0,48 quilotons de carbono por ano. Após a seleção do modelo adequado, testamos os resultados em busca de heterossedasticidade do modelo.

Code

plot <- tibble(fitted(TWFE),resid(TWFE))

ggplot(data = plot, aes(x = `fitted(TWFE)`, y = `resid(TWFE)`^2)) + 
         geom_point() +
         labs(title = "Resíduos ao quadrado") +
         geom_smooth(method = "loess", se = FALSE)+
  theme_bw()

A análise visual aponta para a existencia de heterocedasticidade nos dados. Vamos conduzir um teste de Breusch-Pagan para avaliar se isso se confirma:

Code

bptest(TWFE, studentize=F)


    Breusch-Pagan test

data:  TWFE
BP = 2504.1, df = 10, p-value < 2.2e-16

Confirmando a análise gráfica, o teste aponta para a existência de heterocedasticidade dos dados. Assim,novamente aplicaremos uma matrix de correção de Arellano para obter os erros padrões robustos.

Code

# Obtendo os erros padrões robustos
ARL  <- vcovHC(TWFE, method = "arellano")

TWFE.ARL <- TWFE; TWFE.ARL$vcov <- ARL

tab_model(TWFE, TWFE.ARL,
          dv.labels = c("Modelo sem erros padrões robustos","Modelo com erros padrões robustos"), 
          show.se = T,
          show.ci = F,
          digits = 4
)

	Modelo sem erros padrões robustos			Modelo com erros padrões robustos
Predictors	Estimates	std. Error	p	Estimates	std. Error	p
treat * counter	-0.4835	0.0699	<0.001	-0.4835	0.2000	0.016
treat	0.7001	0.6006	0.244	0.7001	0.4934	0.156
pop density	-0.0890	0.0487	0.068	-0.0890	0.0440	0.044
income per capita	0.0016	0.0006	0.007	0.0016	0.0010	0.098
pop age 65 plus	-0.0571	0.1277	0.655	-0.0571	0.1051	0.587
pop college	-0.5980	0.2168	0.006	-0.5980	0.3074	0.052
value agr mill	1.6331	4.6892	0.728	1.6331	2.9568	0.581
value commerce mill	0.0374	0.0090	<0.001	0.0374	0.0068	<0.001
value manuf mill	0.5169	0.1601	0.001	0.5169	0.0332	<0.001
Observations	744			744
R² / R² adjusted	0.104 / -0.005			0.104 / -0.005

4 Referências

Card, David, and Alan B. Krueger. 1994. “Minimum Wages and Employment: A Case Study of the Fast-Food Industry in New Jersey and Pennsylvania.” The American Economic Review 84 (4): 772–93.

Fraser, Timothy et al. 2020. “Climate Crisis at City Hall: How Japanese communities mobilize to eliminate emissions” Enviromental Innovation and Social Transitions 37

Fraser, Timothy. 2021. Using Difference-in-Differences Models for Environmental Social Science How Interaction Models can help us work with Data over Time. 20201.RPubs

Leppert, Philipp. 2020. “R Tutorial: Difference-in-Differences” RPubs