Lista2

Lista de Exercícios 02

Curso: Tópicos Especiais de Métodos Quantitativos (UFABC)

Conteúdo: Regressão Multivariada, Metodologia de Surveys, Estudos Longitudinais e Indicadores & Medidas Sintéticas,

Autor: Paula Rafaela Ferreira Barbosa

Data: 05 de agosto 2025

Instalar pacotes

if (!requireNamespace(“dplyr”, quietly = TRUE)) install.packages(“dplyr”) if (!requireNamespace(“tidyverse”, quietly = TRUE)) install.packages(“tidyverse”) if (!requireNamespace(“survey”, quietly = TRUE)) install.packages(“survey”) if (!requireNamespace(“fixest”, quietly = TRUE)) install.packages(“fixest”) if (!requireNamespace(“plm”, quietly = TRUE)) install.packages(“plm”) if (!requireNamespace(“did”, quietly = TRUE)) install.packages(“did”) if (!requireNamespace(“lme4”, quietly = TRUE)) install.packages(“lme4”) if (!requireNamespace(“stargazer”, quietly = TRUE)) install.packages(“stargazer”) if (!requireNamespace(“car”, quietly = TRUE)) install.packages(“car”) install.packages(“lmtest”)

Carregar

library(dplyr) library(tidyverse)
library(survey) library(fixest) library(plm) library(did) library(lme4) library(stargazer) library(car) library(ggplot2) library(lmtest)

Exercício 1 — Regressão Multivariada / Tema-guia: Relação entre capital humano e rendimentos individuais

B1 - (B1) População sintética — Gere 10 000 observações com as variáveis descritas no Roteiro Lab-02. Ajuste um modelo OLS de salário contra escolaridade, experiência e sexo. Comente os coeficientes.

Código:

set.seed(123) # Reprodutibilidade N <- 10000 pop <- tibble( sexo = rbinom(N, 1, 0.5),
escolaridade = rnorm(N, 12, 3), # média de 12 anos de estudo experiencia = pmax(rnorm(N, 10, 5), 0), # média de 10 anos de experiência erro = rnorm(N, 0, 4) # erro aleatório ) |> mutate( salario = 800 + 150 * escolaridade + 60 * experiencia - 120 * sexo + erro ) modelo <- lm(salario ~ escolaridade + experiencia + sexo, data = pop) summary(modelo)

Comentário: Os resultados indicam que todas as variáveis são estatisticamente significativas, com p-valores inferiores a 2e-16 (ou seja, a probabilidade de que os coeficientes sejam iguais a zero, é praticamente nula, o que indica que que essas variáveis influenciam o salário). O coeficiente do intercepto foi estimado em 800, representando o salário médio de um indivíduo com zero anos de escolaridade e experiência, e com sexo igual a 0. A escolaridade apresentou coeficiente de 150, indicando que, a cada ano adicional de estudo, o salário aumenta em média R$150, mantendo-se as demais variáveis constantes. A experiência profissional teve coeficiente de aproximadamente 60, o que significa um acréscimo médio de R$60 no salário por ano de experiência. Já a variável sexo apresentou coeficiente de -120, o que implica que indivíduos com sexo igual a 1 recebem, em média, R$120 a menos do que os com sexo igual a 0. O modelo apresentou um R² de 0,9999 (ou seja, 99,99% da variação nos salários observados é explicada pelas variáveis do modelo).

(B2) Diagnósticos fundamentais—Avalie normalidade dos resíduos (QQ-plot), heterocedasticidade (Breusch–Pagan) e multicolinearidade (VIF). Explique cada teste.

qqnorm(residuals(modelo), main = “QQ-plot dos resíduos”) qqline(residuals(modelo), col = “red”) bptest(modelo) # Teste de Breusch–Pagan - Heterocedasticidade vif(modelo) #VIF plot(modelo$fitted.values, residuals(modelo), xlab = “Valores ajustados”, ylab = “Resíduos”, main = “Resíduos vs Ajustados”) abline(h = 0, col = “red”)

Comentário: Os três diagnósticos indicam que o modelo está adequado. O teste de Breusch–Pagan apresentou p-valor de 0,8338 (valor bem acima de 0,05), o que indica que não há heterocedasticidade (a variância dos erros é constante). Os valores de VIF ficaram próximos de 1 para todas as variáveis, o que mostra ausência de multicolinearidade (as variáveis independentes não estão correlacionadas entre si). Por fim, o QQ-plot mostrou que os resíduos seguem aproximadamente uma distribuição normal (os pontos ficaram próximos da linha reta), confirmando mais um pressuposto importante da regressão linear.

(I1) Interação & não-linearidade—Inclua o termo de interação escolaridade × sexo e o termo quadrático de experiência. Interprete mudanças nos efeitos marginais.

set.seed(123) N <- 10000 pop <- tibble( sexo = rbinom(N, 1, 0.5), escolaridade = round(rnorm(N, mean = 12, sd = 3)), # anos de estudo experiencia = round(pmax(rnorm(N, mean = 10, sd = 5), 0)), # truncar negativos erro = rnorm(N, 0, 1000) ) %>% mutate( salario = 1500 + 300 * escolaridade + 200 * experiencia - 5 * experiencia^2 + 500 * sexo + 100 * escolaridade * sexo + erro, sexo_fator = factor(sexo, labels = c(“feminino”, “masculino”)), experiencia2 = experiencia^2 )

modelo_i1 <- lm(salario ~ escolaridade * sexo_fator + experiencia + experiencia2, data = pop)

summary(modelo_i1)

Comentário: O modelo indica que mulheres têm um retorno de R$300 por ano de escolaridade, enquanto homens têm um retorno adicional de R$103, totalizando R$434 por ano.

(I2) Erros robustos & ICs — Re-estime o modelo completo usando erros padrão robustos (HC3). Compare IC 95 % com o OLS ingênuo.

coeftest(modelo_i1, vcov = vcovHC(modelo_i1, type = “HC3”))

ic_robustos <- coefci(modelo_i1, vcov. = vcovHC(modelo_i1, type = “HC3”), level = 0.95)

ic_ols <- confint(modelo_i1)

ic_comparacao <- as_tibble(ic_ols, rownames = “variavel”) %>% rename(OLS_lower = 2.5 %, OLS_upper = 97.5 %) %>% left_join(as_tibble(ic_robustos, rownames = “variavel”) %>% rename(Robust_lower = 2.5 %, Robust_upper = 97.5 %), by = “variavel”)

ic_comparacao

Comentário: Ao aplicar erros padrão robustos (HC3), observou-se que os intervalos de confiança de 95% ficaram geralmente mais largos em comparação ao modelo OLS tradicional. Isso indica que, ao corrigir para possíveis problemas de heterocedasticidade, os erros padrão aumentaram — sinal de que o modelo ingênuo estava subestimando a variabilidade das estimativas.

(A) Dados reais—PNAD Contínua

install.packages(“dplyr”) install.packages(“ggplot2”) install.packages(“PNADcIBGE”) library(PNADcIBGE) library(survey) library(dplyr) library(ggplot2)

pnadc_2024t4 <- get_pnadc( year = 2024, quarter = 4, design = TRUE, # retorna como objeto survey vars = c(“UF”, “VD4019”, “V2007”, “V2009”) # Renda, sexo, idade )

dados_sem_pesos <- pnadc_2024t4$variables %>% filter(!is.na(VD4019), VD4019 > 0)

modelo_simples <- lm( VD4019 ~ V2009 + factor(V2007), data = dados_sem_pesos )

summary(modelo_simples)

modelo_pesado <- svyglm( VD4019 ~ V2009 + factor(V2007), design = pnadc_2024t4 )

summary(modelo_pesado)

Comentário: Variáveis segundo dicionario: Renda do trabalho principal (VD4019), sexo (V2007) e idade (V2009). O modelo ponderado com svyglm() mostra que os resultados mudam quando se respeita a estrutura de amostragem da PNAD Contínua.A idade tem efeito positivo e significativo sobre a renda: indivíduos mais velhos tendem a receber mais, embora esse efeito possa não ser linear indefinidamente.O coeficiente da variável factor(V2007)2 (sexo feminino) é negativo, o que indica que mulheres ganham, em média, menos do que homens, mesmo controlando pela idade.

Exercício 2 — Metodologia de Surveys Tema-guia: Satisfação com serviços públicos

B1 - Amostra simulada — Simule 20 000 indivíduos com variáveis renda, idade, sexo e satisfação (Likert 1–5). Sorteie amostra estratificada por sexo (n=1200) e calcule pesos.

set.seed(123)
N <- 20000

pop <- tibble( id = 1:N, sexo = rbinom(N, 1, 0.5), # 0 = feminino, 1 = masculino idade = round(rnorm(N, mean = 40, sd = 15)), # idade média de 40 anos renda = round(rlnorm(N, meanlog = 8, sdlog = 0.6), 2), # renda com distribuição log-normal satisfacao = sample(1:5, size = N, replace = TRUE, prob = c(0.1, 0.15, 0.3, 0.25, 0.2)) # escala Likert )

table(pop$sexo)

prop_sexo <- prop.table(table(pop$sexo)) n_total <- 1200 n_fem <- round(prop_sexo[“0”] * n_total) # feminino n_mas <- n_total - n_fem # masculino (complementa)

pop_mulheres <- pop %>% filter(sexo == 0) pop_homens <- pop %>% filter(sexo == 1)

amostra_mulheres <- pop_mulheres %>% sample_n(n_fem) amostra_homens <- pop_homens %>% sample_n(n_mas)

amostra <- bind_rows(amostra_mulheres, amostra_homens)

tamanho_pop <- pop %>% count(sexo, name = “N_h”) tamanho_amostra <- amostra %>% count(sexo, name = “n_h”) pesos <- left_join(tamanho_pop, tamanho_amostra, by = “sexo”) %>% mutate(peso = N_h / n_h)

amostra_final <- amostra %>% left_join(pesos %>% select(sexo, peso), by = “sexo”)

head(amostra_final)

Comentário: Foi simulada uma população de 20.000 indivíduos com as variáveis sexo, idade, renda e satisfação (escala Likert de 1 a 5), visando representar um cenário de pesquisa de opinião sobre serviços públicos. A amostra de 1.200 indivíduos foi sorteada de forma estratificada por sexo, preservando a proporção de mulheres e homens existente na população. Os resultados mostraram que cada indivíduo da amostra representa aproximadamente 16,7 pessoas da população.

(B2) Estimadores pontuais—Com o pacote survey, estime média ponderada de satisfação e IC 95 %. Compare com estimativa não-ponderada.

desenho <- svydesign(ids = ~1, data = amostra_final, weights = ~peso)

media_ponderada <- svymean(~satisfacao, design = desenho) media_nao_ponderada <- mean(amostra_final$satisfacao)

media_ponderada media_nao_ponderada summary(amostra_final$peso)

Comentário: A média não ponderada apresentou exatamente o mesmo valor (3,265). Isso ocorreu porque todos os pesos amostrais foram iguais (16,67), resultado da amostragem estratificada proporcional à população. Quando a amostra mantém a mesma distribuição dos estratos da população, ponderar ou não gera o mesmo resultado, já que cada observação representa a mesma fração da população.

(I1) Regressão com pesos—Modele satisfacao ~ renda + idade + sexo. Reporte estatísticas de ajuste adequadas

library(survey) desenho <- svydesign(ids = ~1, data = amostra_final, weights = ~peso) modelo_pesos <- svyglm(satisfacao ~ renda + idade + sexo, design = desenho) summary(modelo_pesos)

modelo_sem_pesos <- lm(satisfacao ~ renda + idade + sexo, data = amostra_final) summary(modelo_sem_pesos) library(broom) comp_coef <- bind_rows( tidy(modelo_sem_pesos) %>% mutate(modelo = “Sem pesos”), tidy(modelo_pesos) %>% mutate(modelo = “Com pesos”) )

comp_coef

Comentário: A comparação mostrou que os coeficientes e p-valores foram praticamente idênticos nos dois modelos. Isso ocorre porque, na etapa (B1), a amostra foi obtida por amostragem estratificada proporcional por sexo, resultando em pesos iguais (~16,67) para todos os indivíduos. Nesse cenário, ponderar ou não leva a estimativas quase idênticas. Nenhuma das variáveis explicativas (renda, idade e sexo) apresentou significância estatística ao nível de 5%, indicando que, na amostra simulada, não há evidência de associação linear entre esses fatores e a satisfação. Apenas o intercepto foi significativo, com valor estimado de 3,39, representando a satisfação média para o grupo de referência (mulheres com renda e idade iguais a zero).

(I2) Pós-stratificação — Ajuste pesos via raking para refletir proporção real (55 % mulheres). Re-estime a média de satisfação.

amostra_final <- amostra_final %>% mutate(sexo_fator = factor(sexo, labels = c(“feminino”, “masculino”))) desenho_base <- svydesign(ids = ~1, data = amostra_final, weights = ~peso) marginais <- list(~sexo_fator) proporcoes <- list(data.frame(sexo_fator = c(“feminino”, “masculino”), Freq = c(0.55, 0.45))) desenho_rake <- rake(desenho_base, sample.margins = marginais, population.margins = proporcoes)

media_rake <- svymean(~satisfacao, design = desenho_rake)

media_rake

Comentário: # Comentário: média ponderada de satisfação foi de 3,26

Dados reais — LAPOP 2023

lapop <- read_dta(“C:/Users/paula/Desktop/BRA_2023_LAPOP_AmericasBarometer_v1.0_w.dta”) names(lapop)

lapop <- lapop %>% filter(!is.na(q10inc), !is.na(b2)) %>% # usar q10inc filter(!(q10inc %in% c(888888, 988888))) %>% mutate(quintil_renda = ntile(q10inc, 5)) # criar quintis de renda

names(lapop)[1:50]

lapop <- lapop %>% filter(!is.na(q10inc), !is.na(b13)) %>% # manter casos completos filter(!(q10inc %in% c(888888, 988888))) # remover “não sabe” e “não responde”

lapop <- lapop %>% mutate(quintil_renda = ntile(q10inc, 5))

desenho_lapop <- svydesign( ids = ~upm, strata = ~estratopri, weights = ~wt, data = lapop, nest = TRUE )

conf_por_quintil <- svyby( ~b13, ~quintil_renda, design = desenho_lapop, svymean, na.rm = TRUE )

print(conf_por_quintil)

ggplot(conf_por_quintil, aes(x = factor(quintil_renda), y = b13)) + geom_col(fill = “steelblue”) + geom_errorbar(aes(ymin = b13 - 1.96se, ymax = b13 + 1.96se), width = 0.2) + labs( title = “Confiança no Congresso por quintil de renda - LAPOP Brasil 2023”, x = “Quintil de renda”, y = “Média de confiança (1 = nada, 7 = muito)” )

Comentário: Exercicio complexo, não estou confiando na resposta, mas os estratos de renda mais alta tendem a apresentar menor confiança no Congresso Nacional em comparação aos estratos de renda mais baixa, embora as diferenças entre os três primeiros quintis sejam pouco expressivas.

Exercício 3 — Estudos Longitudinais Tema-guia: Impacto de programa de qualificação profissional

set.seed(123) n_ind <- 1200
n_per <- 6

painel <- expand.grid( id = 1:n_ind, periodo = 1:n_per ) %>% arrange(id, periodo)

caracteristicas <- tibble( id = 1:n_ind, idade_inicial = sample(18:40, n_ind, replace = TRUE), experiencia_inicial = sample(0:10, n_ind, replace = TRUE), tratamento = rbinom(n_ind, 1, 0.5) # 0 ou 1 )

painel <- painel %>% left_join(caracteristicas, by = “id”) %>% mutate( idade = idade_inicial + (periodo - 1), experiencia = experiencia_inicial + (periodo - 1), salario = 1200 + 20 * idade + 50 * experiencia + 200 * tratamento * (periodo >= 4) + # efeito a partir do período 4 rnorm(n(), mean = 0, sd = 300) )

head(painel)

medias <- painel %>% group_by(tratamento, periodo) %>% summarise(media_salario = mean(salario), .groups = “drop”) %>% mutate(grupo = ifelse(tratamento == 1, “Tratamento”, “Controle”))

print(medias)

ggplot(medias, aes(x = periodo, y = media_salario, color = grupo, group = grupo)) + geom_line(size = 1.2) + geom_point(size = 2) + labs( title = “Média salarial por período - Tratamento vs Controle”, x = “Período”, y = “Média salarial” ) + theme_minimal() + scale_color_manual(values = c(“Tratamento” = “blue”, “Controle” = “red”))

Comentário: Ambos os grupos apresentam crescimento salarial ao longo do tempo, explicado pelo aumento da idade e da experiência, mas o grupo Tratamento cresce mais rápido após o início do programa.

(B2) Efeitos fixos — Estime modelo FE e interprete coeficiente de tratamento.

library(plm)

painel <- painel %>% mutate(id = as.factor(id), periodo = as.factor(periodo))

painel <- painel %>% mutate(pos = ifelse(as.numeric(periodo) >= 4, 1, 0), tratamento_pos = tratamento * pos) # interação

modelo_fe <- plm( salario ~ tratamento_pos + idade + experiencia, data = painel, index = c(“id”, “periodo”), model = “within” )

summary(modelo_fe)

O programa teve impacto positivo e significativo nos salários dos participantes, aumentando-os em cerca de R$ 200 por período.

(I1) Hausman — Compare FE versus RE via teste de Hausman.

modelo_fe <- plm( salario ~ tratamento_pos + idade, data = painel, index = c(“id”, “periodo”), model = “within” )

modelo_re <- plm( salario ~ tratamento_pos + idade, data = painel, index = c(“id”, “periodo”), model = “random” )

hausman <- phtest(modelo_fe, modelo_re) hausman

Comentário: Rejeita‑se a hipótese nula de que o modelo de efeitos aleatórios é consistente. Assim, conclui‑se que o modelo de efeitos fixos é o mais adequado.

(I2) DiD escalonado — Use did::att_gt (R) ou python-did. Grafique ATTg,t.

library(did)

painel <- painel %>% mutate( id = as.numeric(as.character(id)), # converter id para numérico periodo_num = as.numeric(as.character(periodo)), first_treat = ifelse(tratamento == 1, 4, 0) # entrada no tratamento )

did_result <- att_gt( yname = “salario”, tname = “periodo_num”, idname = “id”, gname = “first_treat”, data = painel, panel = TRUE )

summary(did_result) ggdid(did_result)

Comentario: O programa de qualificação profissional elevou o salário médio dos participantes em cerca de R$ 200 por período a partir de sua implementação no período 4. Não foram detectadas diferenças significativas antes da implementação, o que reforça a validade da identificação causal pelo método.

Exercicio 4: Exercício 4 — Indicadores & Medidas Sintéticas Tema-guia: Índice Municipal de Vulnerabilidade Social (IMVS)

(B1) Dados simulados — Gere 500 municípios com indicadores padronizados: renda per capita, analfabetismo, saneamento, mortalidade infantil e desemprego.

set.seed(123)

n_mun <- 500

dados_imvs <- tibble( municipio = paste0(“Município_”, 1:n_mun), renda_pc = rnorm(n_mun, mean = 1200, sd = 300), # Renda per capita (R$) analfabetismo = rnorm(n_mun, mean = 8, sd = 3), # % saneamento = rnorm(n_mun, mean = 75, sd = 10), # % mortalidade_infantil = rnorm(n_mun, mean = 15, sd = 5), # por mil nascidos vivos desemprego = rnorm(n_mun, mean = 10, sd = 4) # % )

dados_padronizados <- dados_imvs %>% mutate( across( renda_pc:desemprego, ~ scale(.)[,1], # padronização .names = “{.col}_z” ) )

head(dados_padronizados)

(B2) Índice simples — Calcule IMVS pela média aritmética e classifique em quintis.

dados_imvs <- dados_padronizados %>% mutate( imvs = rowMeans(select(., renda_pc_z, analfabetismo_z, saneamento_z, mortalidade_infantil_z, desemprego_z), na.rm = TRUE), quintil_imvs = ntile(imvs, 5) # 1 = menor índice, 5 = maior índice )

head(dados_imvs %>% select(municipio, imvs, quintil_imvs))

Comentário: Quanto maior o IMVS, mais altos os valores médios dos indicadores padronizados.Quanto menor o IMVS, menores os valores médios dos indicadores padronizado

(I1) PCA — Aplique PCA (autovalor ≥ 1) e construa índice ponderado pelas cargas.

indicadores <- dados_padronizados %>% select(renda_pc_z, analfabetismo_z, saneamento_z, mortalidade_infantil_z, desemprego_z)

pca_result <- prcomp(indicadores, scale. = FALSE)

autovalores <- pca_result$sdev^2 print(autovalores)

componentes_selecionados <- which(autovalores >= 1) print(componentes_selecionados)

cargas <- pca_result$rotation[, componentes_selecionados] print(cargas)

var_exp <- autovalores[componentes_selecionados] / sum(autovalores[componentes_selecionados]) print(var_exp)

indice_pca <- as.matrix(indicadores) %% as.matrix(cargas %% diag(var_exp))

indice_pca <- rowSums(indice_pca)

indice_pca <- as.numeric(scale(indice_pca))

dados_pca <- dados_padronizados %>% mutate(indice_pca = indice_pca)

head(dados_pca %>% select(municipio, indice_pca))

Comentário: Valores positivos = municípios mais vulneráveis, considerando a combinação ponderada dos indicadores. Valores negativos = municípios menos vulneráveis no conjunto.