library(tidyverse) # Conjunto de pacotes e funções para realizar análise de dados no R
library(readxl) # Para fazer leitura e também exportar arquivos em formato Excel.
library(gt) # Para elaboração de tabelas com qualidade de publicação
library(here) # Para facilitar o trabalho com caminhos de pastas
library(sandwich) # vcovHC() para erros-padrão robustos
library(lmtest) # coeftest()Aula Prática - Inferência Causal Baseada em Regressões
Introdução
Antes de prosseguir com a aula prática sobre modelos de regressão, iremos assistir a dois vídeos. O primeiro vídeo apresenta de modo intuitivo a ideia de ceteris paribus em um contexto de regressão. Já o segundo vídeo discute em maior detalhe o problema do viés.
Ceteris Paribus
Viés de Seleção = Viés de Variável Omitida em Regressão
Dale, S., Krueger, A. Estimating the Payoff to Attending a More Selective College: An Application of Selection on Observables and Unobservables. The Quarterly Journal of Economics, Volume 117, Issue 4, November 2002, Pages 1491–1527.
Descrição da atividade
A ideia para a aula prática de hoje será reproduzir a Tabela 7.1 do livro de Stock e Watson, apresentada abaixo.
Observação: na 1a edição do livro, é a Tabela 5.2.
- Fazer download da base de dados e do dicionário de variáveis. Importar dados para o R.
Antes de iniciar a atividade, vamos chamar os pacotes que serão utilizados. Nesta resolução, são utilizados os seguintes pacotes do R. 1
Agora, vamos importar os dados e salvá-los no objeto chamado caschool.
# Carregar os dados
# Ajuste o caminho conforme a localização do arquivo na sua máquina
caschool <- read_excel(here::here("labs","SW_Datasets","caschool.xlsx"))- Estimar os modelos de regressão apresentados nas colunas 1 a 5 da tabela.
Abaixo estimamos as 5 especificações da tabela utilizando o comando lm. Em todos os casos, salvamos o resultado em um objeto para que possamos aplicar outras funções a ele para a obtenção de outros resultados de interesse.
#Especificação 1
espec_1 <- lm(testscr ~ str, data = caschool)
#Especificação 2
espec_2 <- lm(testscr ~ str + el_pct, data = caschool)
#Especificação 3
espec_3 <- lm(testscr ~ str + el_pct + meal_pct, data = caschool)
#Especificação 4
espec_4 <- lm(testscr ~ str + el_pct + calw_pct, data = caschool)
#Especificação 5
espec_5 <- lm(testscr ~ str + el_pct + meal_pct+ calw_pct, data = caschool)Uma vez que o modelo estimado foi salvo, podemos ver os resultados principais com as opções padrão utilizando a função summary. Vejamos o caso da primeira especificação:
#Especificação 1
summary(espec_1)
Call:
lm(formula = testscr ~ str, data = caschool)
Residuals:
Min 1Q Median 3Q Max
-47.727 -14.251 0.483 12.822 48.540
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 698.9330 9.4675 73.825 < 2e-16 ***
str -2.2798 0.4798 -4.751 2.78e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 18.58 on 418 degrees of freedom
Multiple R-squared: 0.05124, Adjusted R-squared: 0.04897
F-statistic: 22.58 on 1 and 418 DF, p-value: 2.783e-06Vocês notam alguma diferença em relação à tabela? Os desvios-padrão dos coeficientes são diferentes! Isso porque a tabela já reporta desvios-padrão robustos, ou seja, que já contemplam a possibilidade de termos um modelo heterocedástico. Para ver o resultado com desvio-padrão robusto, utilizaremos a função coeftest e o pacote sandwich para obter o resultado:2
#Especificação 1
coeftest(espec_1, vcov. = vcovHC(espec_1,"HC1"))
t test of coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 698.93295 10.36436 67.4362 < 2.2e-16 ***
str -2.27981 0.51949 -4.3886 1.447e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Podemos proceder de forma análoga para obter os resultados das demais especificações. Não faremos isso nesta etapa, pois iremos trabalhar no item seguinte para já montar uma tabela de resultados.
- Utilizando o pacote
gte seu estagiário de IA, prepare a tabela para reportar os resultados. Não se preocupe em reproduzir exatamente a formatação, desde que apresente as informações principais (coeficientes, desvios-padrão robustos). Não há necessidade de mostrar na tabela o intervalo de confiança, como é feito apenas para a variável student-teacher ratio no livro.
Abaixo, seguem o código e a Tabela 1 com os resultados gerados a partir dos modelos estimados.3
# ============================================================
# Tabela de Regressões - Desempenho em Testes (testscr) e
# Razão Aluno-Professor (str)
# Dataset: caschool (Stock & Watson)
#
# 5 especificações, erros-padrão robustos (HC1), formato de
# publicação (sem asteriscos de significância).
# ============================================================
# ------------------------------------------------------------
# 1. Especificações
# ------------------------------------------------------------
modelos <- list(espec_1, espec_2, espec_3, espec_4, espec_5)
nomes_colunas <- paste0("(", seq_along(modelos), ")")
# ------------------------------------------------------------
# 2. Ordem e rótulos das variáveis explicativas
# ------------------------------------------------------------
ordem_termos <- c("str", "el_pct", "meal_pct", "calw_pct", "(Intercept)")
rotulos_termos <- c(
str = "Razão aluno-professor (str)",
el_pct = "% de aprendizes de inglês (el_pct)",
meal_pct = "% de alunos com refeição gratuita (meal_pct)",
calw_pct = "% de famílias no CalWORKS (calw_pct)",
`(Intercept)` = "Constante"
)
# ------------------------------------------------------------
# 3. Extrair coeficientes e erro-padrão robusto (HC1) de cada modelo
# ------------------------------------------------------------
extrair_modelo <- function(modelo) {
rob <- coeftest(modelo, vcov = vcovHC(modelo, type = "HC1"))
list(
est = rob[, "Estimate"],
se = rob[, "Std. Error"],
n = nobs(modelo),
r2 = summary(modelo)$r.squared
)
}
resultados <- map(modelos, extrair_modelo)
formata <- function(x, casas = 3) {
if (is.na(x)) "" else sprintf(paste0("%.", casas, "f"), x)
}
# ------------------------------------------------------------
# 4. Montar o corpo da tabela: cada variável em duas linhas
# (estimativa e erro-padrão entre parênteses)
# ------------------------------------------------------------
linhas_coef <- map_dfr(ordem_termos, function(termo) {
estimativas <- map_chr(resultados, function(r) {
if (termo %in% names(r$est)) formata(r$est[[termo]]) else ""
})
erros <- map_chr(resultados, function(r) {
if (termo %in% names(r$se)) paste0("(", formata(r$se[[termo]]), ")") else ""
})
bind_rows(
setNames(as.list(c(rotulos_termos[[termo]], estimativas)), c("variavel", nomes_colunas)),
setNames(as.list(c("", erros)), c("variavel", nomes_colunas))
)
})
# ------------------------------------------------------------
# 5. Linhas de rodapé: número de observações e R²
# ------------------------------------------------------------
linha_n <- setNames(
as.list(c("Observações", map_chr(resultados, ~ as.character(.x$n)))),
c("variavel", nomes_colunas)
)
linha_r2 <- setNames(
as.list(c("R²", map_chr(resultados, ~ formata(.x$r2)))),
c("variavel", nomes_colunas)
)
tabela_final <- bind_rows(linhas_coef, linha_n, linha_r2)
# ------------------------------------------------------------
# 6. Construir a tabela com gt (estilo simples, tipo periódico)
# ------------------------------------------------------------
tab_reg <- tabela_final %>%
gt() %>%
tab_header(
title = md("**Tabela.** Desempenho em Testes e Razão Aluno-Professor")
) %>%
cols_label(variavel = "") %>%
cols_align(align = "center", columns = all_of(nomes_colunas)) %>%
tab_spanner(
label = "Variável dependente: testscr",
columns = all_of(nomes_colunas)
) %>%
tab_style(
style = cell_borders(sides = "top", color = "black", weight = px(1)),
locations = cells_body(rows = variavel == "Observações")
) %>%
tab_options(
table.border.top.style = "solid",
table.border.top.width = px(1.5),
table.border.bottom.style = "solid",
table.border.bottom.width = px(1.5),
column_labels.border.top.style = "solid",
column_labels.border.top.width = px(1),
column_labels.border.bottom.style = "solid",
column_labels.border.bottom.width = px(1),
table_body.hlines.style = "none",
heading.border.bottom.style = "none",
column_labels.font.weight = "bold",
table.font.size = px(11),
heading.title.font.size = px(12),
source_notes.font.size = px(9),
data_row.padding = px(2)
) %>%
tab_source_note(
source_note = "Nota: erros-padrão robustos à heterocedasticidade entre parênteses."
)
tab_reg| Tabela. Desempenho em Testes e Razão Aluno-Professor | |||||
Variável dependente: testscr
|
|||||
|---|---|---|---|---|---|
| (1) | (2) | (3) | (4) | (5) | |
| Razão aluno-professor (str) | -2.280 | -1.101 | -0.998 | -1.308 | -1.014 |
| (0.519) | (0.433) | (0.270) | (0.339) | (0.269) | |
| % de aprendizes de inglês (el_pct) | -0.650 | -0.122 | -0.488 | -0.130 | |
| (0.031) | (0.033) | (0.030) | (0.036) | ||
| % de alunos com refeição gratuita (meal_pct) | -0.547 | -0.529 | |||
| (0.024) | (0.038) | ||||
| % de famílias no CalWORKS (calw_pct) | -0.790 | -0.048 | |||
| (0.068) | (0.059) | ||||
| Constante | 698.933 | 686.032 | 700.150 | 697.999 | 700.392 |
| (10.364) | (8.728) | (5.568) | (6.920) | (5.537) | |
| Observações | 420 | 420 | 420 | 420 | 420 |
| R² | 0.051 | 0.426 | 0.775 | 0.629 | 0.775 |
| Nota: erros-padrão robustos à heterocedasticidade entre parênteses. | |||||
Utilizando a tabela e interpretando os resultados apresentados, responda às seguintes perguntas:
Utilizando os resultados da Tabela 1, seguem sugestões de respostas breveds para as perguntas realizadas em cada item. Para uma análise mais aprofundada de todos os resultados, sugiro fortemente a leitura do capítulo do livro em que os resultados são discutidos detalhadamente.
Qual especificação do modelo você julga ser a mais apropriada?
Sugestão de resposta: A especificação (3) poderia ser escolhida como especificação mais apropriada, uma vez que, além da variável de interesse, controla por características dos alunos que poderiam gerar viés de variável omitida caso não fossem incluídas no modelo. As especificações (4) e (5) apresentam formas alternativas de se controlar pelas condições socioeconômicas dos alunos. Note que os resultados não mudam muito entre as especificações (4) e (5).
Utilizando a especificação escolhida, qual o impacto de reduzir o tamanho médio das turmas em 1 aluno? Interprete o valor do coeficiente e também discuta o resultado em termos de inferência estatística.
Sugestão de resposta: O resultado reportado na coluna (3) indica que, ao reduzir a razão aluno-professor em um aluno por professor, a nota média no teste aumenta aproximadamente 1 ponto, mantendo-se constantes as características dos alunos. Note que o resultado é robusto, ou seja, o coeficiente é relativamente similar nas demais especificações com controles.
Existe relação causal entre tamanho das turmas e notas? Caso sua resposta seja sim, explique qual foi a hipótese de identificação que você utilizou. Se não, diga qual hipótese o modelo não atende e por que ela não seria atendida.
Sugestão de resposta: Caso responda sim, o aluno deve justificar por que o modelo não sofre de viés de variável omitida, ou seja, que não há nenhuma variável não incluída no modelo que seja correlacionada com a razão aluno-professor e com as notas médias do teste.
Caso responda não, o aluno deve apresentar outros mecanismos omitidos que estejam correlacionados com a razão aluno-professor e com as notas médias do teste.
Imagine que você seja prefeito de uma cidade no Brasil. Seu secretário de educação lhe apresenta os resultados desse estudo utilizando os dados da Califórnia e sugere que vocês tentem reduzir o tamanho médio das turmas no município com o intuito de melhorar o desempenho escolar. O que você acharia dessa proposta? Justifique, apresentando argumentos de validade interna ou externa.
Sugestão de resposta: Eu pediria ao secretário para me apresentar por quais motivos ele acreditava que os resultados do estudo feito para Califórnia poderiam ser extrapolados para a população do município, dadas as grandes diferenças socioeconômicas, culturais e institucionais entre o estado americano e o nosso país. Caso a explicação não fosse convincente, eu chegaria à conclusão de que o argumento teria um problema de validade externa e, portanto, não se aplicaria à população residente no município.
Anexo - Prompt
Abaixo, segue o prompt enviado ao Claude para elaboração da tabela em formato de publicação.
Prompt:
Estou trabalhando em um arquivo quarto com o header YAML abaixo. O objetivo é montar uma tabela usando o pacote gt do R. Esta tabela deve reportar o resultados das 5 especificações abaixo em formato para publicação em periódico científico. Note que meu arquivo exporta tanto para o formato html quanto para pdf. As cinco especificações também seguem abaixo. Os pacotes a serem utilizados já foram importados no arquivo. Os desvios-padrão a serem reportados devem ser desvios-padrão robusto a heteroscedásticidade. Além disso, não há necessidade de incluir * para indicar significância estatística.
#Especificação 1 espec_1 <- lm(testscr ~ str, data = caschool)
#Especificação 2 espec_2 <- lm(testscr ~ str + el_pct, data = caschool)
#Especificação 3 espec_3 <- lm(testscr ~ str + el_pct + meal_pct, data = caschool)
#Especificação 4 espec_4 <- lm(testscr ~ str + el_pct + calw_pct, data = caschool)
#Especificação 5 espec_5 <- lm(testscr ~ str + el_pct + meal_pct+ calw_pct, data = caschool)
Header YAML:
---
title: "Aula Prática - Inferência Causal Baseada em Regressões"
reference-location: margin
format:
html:
echo: true
warning: false
error: false
pdf:
echo: false
warning: false
error: false
babel-lang: pt
---Notas de rodapé
Como não quero reportar as mensagens de erro e avisos que o R reproduz após a execução do comando
library, inseri opções específicas para esse bloco de código. Para mais detalhes, veja a página de referência sobre opções de execução.↩︎Alternativamente, é possível utilizar o pacote
estimatrpara obter o resultado de modo mais direto. Veja a documentação do pacote aqui.↩︎Para essa tabela, eu pedi auxílio ao Claude. Veja o prompt inicial no anexo. Depois, revisei o código e fiz alguns ajustes para corrigir erros e detalhes de formatação.↩︎