Efeitos do recebimento de Bolsa Família na Oferta de Trabalho
Análise dos efeitos do recebimento de PBF nas Margens Intensiva e Extensiva usando dados da PNADc Ano de 2024
Luan, Guilherme, Raquel, Enrico
2025-07-14
Introdução e Objetivos
Este relatório analisa os efeitos do recebimento de proventos via Programa Bolsa Família sobre a participação no mercado de trabalho brasileiro, utilizando dados da Pesquisa Nacional por Amostra de Domicílios Contínua (PNADC) de 2024. A análise examina tanto a margem extensiva (probabilidade de participação na força de trabalho) quanto a margem intensiva (horas trabalhadas pelos ocupados).
Objetivos específicos:
- Avaliar a existência e quantificar o efeito do recebimento de BF sobre a probabilidade de não participar da força de trabalho (margem extensiva)
- Avaliar a existência e quantificar o impacto do Bolsa Família sobre as horas trabalhadas (margem intensiva)
- Identificar efeitos heterogêneos por características sociodemográficas, domiciliares, regionais e ocupacionais
- Fornecer evidências empíricas para o debate sobre políticas de transferência de renda e seus impactos na oferta de trabalho
Metodologia
Carregamento de Pacotes
# Lista dos pacotes utilizados
packages <- c("PNADcIBGE", "dplyr", "tidyr", "survey", "jtools", "MASS", "DataExplorer",
"lares", "performance", "emmeans", "sjstats", "srvyr", "broom", "modelsummary",
"margins", "marginaleffects", "ggplot2", "lmtest", "sandwich", "AER", "car",
"knitr", "kableExtra", "tseries", "skedastic", "whitestrap", "aod")
# Carregando ou instalando+carregando
for (pkg in packages) {
if (!requireNamespace(pkg, quietly = TRUE)) {
install.packages(pkg, dependencies = TRUE)
}
library(pkg, character.only = TRUE)
}
rm(pkg)
rm(packages)
"%ni%" <- Negate("%in%")Fonte dos Dados
A Pesquisa Nacional por Amostra de Domicílios Contínua (PNADC) é coletada por meio de uma amostra probabilística complexa que envolve múltiplos estágios de seleção. A amostra é distribuída ao longo de 15 meses, com cada domicílio sendo entrevistado por 5 trimestres consecutivos (esquema de painel rotativo). Em cada domicílio, são coletadas informações de ordem individual.
Importando Dados
Para este estudo, foram utilizados os dados da primeira entrevista de cada domicílio em 2024 (dados anuais de 2024), importados para o ambiente R através do pacote PnadcIBGE.
# Function to get PNADC data
get_pnadc_data <- function(year, interview, vars) {
get_pnadc(
year = year,
interview = interview,
vars = vars,
design = FALSE,
labels = FALSE
)
}
pnadc_vars <- c(
############### Filtros
"VD4003", ## Condição em relação à FT Potencial, Pessoa na FT Potencial = 1
############### Variáveis dependentes
## Participação no Mercado de Trabalho (Margem Extensiva)
"VD4001", ## Condição em relação à FT, Pessoa na FT = 1
"VD4002", ## Condição em relação a ocupação, Ocupados = 1
## Participação no Mercado de Trabalho (Margem Intensiva)
"VD4031", ## Horas habituais de trabalho (numérico)
"VD4035", ## Horas efetivas de trabalho (numérico)
## Participação no Mercado de Trabalho (remuneração)
"VD4016", ## Rendimento habitual do trabalho principal
"VD4017", ## Rendimento efetivo do trabalho principal
"VD4019", ## Rendimento habitual de todos os trabalhos
"VD4020", ## Rendimento efetivo de todos os trabalhos
############### Variáveis explicativas
"VD4047", ## Rendimento efetivo recebido de programas sociais
"V5002A", ## Se recebeu BF
"V5002A2", ## Valor efetivamente recebido via BF
############### Variáveis de controle
"UF",
"V1022", ## Urbano = 1, Rural = 2
"V1023", ## Capital = 1, Resto da RM = 2, Resto da RIDE = 3, Resto da UF = 4
"V2007", ## Sexo, Homem = 1
"V2009", ## Idade (numérico)
"V2010", ## Cor ou raça
"VD2002", ## Condição no domicílio
"VD2003", ## Número de pessoas no domicílio (numérico)
"VD3006", ## Faixa de anos de estudos
"VD4009", ## Posição na Ocupação
"VD4010", ## Grupamento de atividade do trabalho principal
"VD4011", ## Grupamento de ocupação do trabalho principal
"VD5005", ## Rendimento domiciliar per capita
"VD5006", ## Faixa de rendimento domiciliar per capita
############### Pesos
"V1032",
sprintf("V1032%03d", seq(1:200))
)
dados_pnadc <- get_pnadc_data(year = 2024, interview = 1, vars = pnadc_vars)
cat("Dimensões da base de dados:", dim(dados_pnadc), "\n")## Dimensões da base de dados: 381815 245Modelos Escolhidos
Margem Extensiva:
- Como a variável dependente (na_forca_trabalho) é binária, optou-se por utilizar um Modelo Linear Generalizado do tipo binomial/Logit. O Logit foi escolhido em detrimento de outras opções binárias por facilitar a interpretação
Margem Intensiva:
- Como a variável dependente (horas_ef) é uma contagem, ou seja, é numérica mas não é contínua, avaliou-se a possibilidade de utilizar um Modelo Linear Generalizado do tipo Poisson. Contudo, foi identificado que a média e a variância da variável dependente são diferentes; e foi testado que, sob Poisson, havia problema de superdispersão. Por isso, optou-se por utilizar um Modelo Linear Generalizado do tipo binomial negativa
Pesos amostrais:
- Para facilitar os cálculos e devido ao tamanho da amostra (superior a 100.000 registros), optou-se por realizar a modelagem com base na amostra não ponderada pelos pesos amostrais. Ou seja, não usamos as funções do pacote survey para calibrar a amostra e para realizar a modelagem. Esse é um esforço que pode ser empreendido futuramente. Nos testes realizados, os cálculos com base na amostra ponderada e não ponderada chegaram a resultados semelhantes
Seleção de variáveis:
- As variáveis de controle foram escolhidas a partir de base teórica, de fatores sociodemográficos, domiciliares, geográficos e ocupacionais que poderiam ser relevantes na determinação da chance de participação do mercado de trabalho e das horas trabalhadas.
- As variáveis de interação foram selecionadas a partir de testagens. Foram consideradas todas interações possíveis com a variável explicativa de interesse (recebeu_bf). Foram mantidas as interações cujo impacto ao modelo em termos de ajuste (BIC) foram mais relevantes. Também foram mantidas apenas as variáveis de interação aprovadas em teste de significância individual (Wald)
Efeitos Marginais
Como as modelagens realizadas possuem muitas variáveis categóricas de controle, avaliamos que a interpretação dos coeficientes como saem no modelo seria insuficiente - pois a comparação seria sempre em relação ao grupo de referência. Por isso, usamos o pacote marginaleffects para calcular os efeitos marginais médios de cada variável, ou seja, considerando o efeito médio em relação a todas as outras categorias e não apenas a de referência. No caso da modelagem para a margem extensiva, utilizou-se a função comparisons(), indicada pelo pacote para variáveis dependentes binárias; e para a margem intensiva, utilizou-se a função avg_slopes()
De acordo com os testes de Breusch-Pagan, Goldfeld-Quandt e White realizados com os modelos finais, há problema de heterocedasticidade em ambos. Por isso, os efeitos marginais foram calculados usando erros padrão robustos (“HC1”)
analyze_and_plot_marginal_effects <- function(model, variable, by_vars, vcov_matrix,
outcome_label, plot_title_prefix, comparison_type = "slopes") {
cat(paste0("\n=== EFEITO DO BOLSA FAMÍLIA POR ", toupper(paste(by_vars, collapse = " e ")),
" ===\n"))
# Calcular efeitos marginais
if (comparison_type == "slopes") {
effects <- avg_slopes(model, variables = variable, by = by_vars, type = "response",
vcov = vcov_matrix)
} else if (comparison_type == "comparisons") {
effects <- comparisons(model, variables = variable, by = by_vars, type = "response",
vcov = vcov_matrix)
} else {
stop("Erro no cálculo dos efeitos marginais.")
}
# Printar efeitos marginais
cat("\n--- EFEITOS MARGINAIS ---\n")
print(effects)
cat("\n")
# Calcular predições
preds <- avg_predictions(model, vcov = vcov_matrix, by = c(variable, by_vars),
type = "response")
# Printar predições
cat("--- VALORES PREDITOS ---\n")
print(preds)
cat("\n")
# Verificar consistência entre efeitos marginais e predições
cat("--- VERIFICAÇÃO DE CONSISTÊNCIA ---\n")
if (comparison_type == "comparisons") {
# Para comparisons, verificar se as diferenças batem
preds_wide <- preds %>%
pivot_wider(names_from = all_of(variable), values_from = estimate, names_prefix = paste0(variable,
"_"))
if (ncol(preds_wide) >= 4) {
# Tem colunas para ambos os valores
col_0 <- paste0(variable, "_0")
col_1 <- paste0(variable, "_1")
if (col_0 %in% names(preds_wide) && col_1 %in% names(preds_wide)) {
preds_wide$diff_calculada <- preds_wide[[col_1]] - preds_wide[[col_0]]
# Merge com efeitos marginais para comparação
comparison_data <- merge(preds_wide, effects, by = by_vars, suffixes = c("_pred",
"_effect"))
cat("Comparação entre diferenças calculadas e efeitos marginais:\n")
comparison_subset <- comparison_data %>%
select(all_of(by_vars), diff_calculada, estimate) %>%
mutate(diferenca_absoluta = abs(diff_calculada - estimate))
print(comparison_subset)
cat("\n")
}
}
}
# Gráfico de predições
plot_preds <- preds %>%
mutate(Recebeu_BF = ifelse(.data[[variable]] == 1, "Recebe Bolsa Família",
"Não Recebe")) %>%
ggplot(aes(x = .data[[by_vars]], y = estimate, color = Recebeu_BF, group = Recebeu_BF)) +
geom_line(linewidth = 1.2) + geom_point(size = 3) + geom_ribbon(aes(ymin = conf.low,
ymax = conf.high, fill = Recebeu_BF), alpha = 0.2, color = NA) + labs(title = paste(outcome_label,
"Preditas -", plot_title_prefix), x = by_vars, y = outcome_label, color = NULL,
fill = NULL) + scale_color_manual(values = c("#4DAF4A", "#984EA3")) + scale_fill_manual(values = c("#4DAF4A",
"#984EA3")) + theme_minimal(base_size = 12) + theme(legend.position = "top",
axis.text.x = element_text(angle = 45, hjust = 1))
print(plot_preds)
# Gráfico de efeitos marginais
plot_effects <- effects %>%
mutate(Significativo = ifelse(p.value < 0.05, "Sim", "Não")) %>%
ggplot(aes(x = .data[[by_vars]], y = estimate, color = Significativo)) +
geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") + geom_point(size = 3) +
geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2) + labs(title = paste("Efeito do Bolsa Família -",
plot_title_prefix), subtitle = paste("[Recebeu vs Não Recebeu]"), x = by_vars,
y = "Efeito Marginal", color = "Estatisticamente\nSignificativo (p<0.05)") +
scale_color_manual(values = c(Sim = "#E41A1C", Não = "#377EB8")) + theme_minimal(base_size = 12) +
theme(legend.position = "top", axis.text.x = element_text(angle = 45, hjust = 1))
print(plot_effects)
# Testes de hipóteses
num_levels <- nrow(effects)
if (num_levels > 1) {
hyp_combinations <- combn(1:num_levels, 2, simplify = FALSE)
hyp_tests <- c()
for (comb in hyp_combinations) {
hyp_tests <- c(hyp_tests, paste0("b", comb[1], " = b", comb[2]))
}
hyp_tests <- c(hyp_tests, paste0("b", 1:num_levels, " = 0"))
hypothesis_result <- effects %>%
hypotheses(hyp_tests, joint = FALSE)
cat(paste0("Teste Wald (robusto) de diferenças entre ", tolower(by_vars),
":\n"))
print(hypothesis_result)
cat("\n")
}
return(list(effects = effects, predictions = preds, plot_preds = plot_preds,
plot_effects = plot_effects))
}Margem Extensiva
Preparação dos Dados
Variável Dependente:
- Condição em relação à força de trabalho (VD4001, na_forca_trabalho)
Variável Explicativa Principal:
- Recebimento do Bolsa Família (V5002A, recebeu_bf)
Variáveis de Controle:
- Características individuais: sexo (V2007), idade (V2009 agrupada por quartil), raça (V2010 agrupada), escolaridade (VD3006 agrupada, faixa de anos de estudos)
- Características domiciliares: condição no domicílio (VD2002 agrupada), renda domiciliar per capita (VD5005 agrupada por quartil, efetivo, inclusive rendimentos de tíquete transporte ou alimentação, exclusive pessoas do domicílio cuja condição era pensionista, doméstico ou parente de doméstico)
- Características geográficas: região (UF agrupada) e área urbana (V1022)
Agrupamentos realizados nas variáveis:
- Anos Estudo (Até 11, 12 a 15, 16 ou mais)
- Raça (Brancos, Pretos e Pardos, Outros)
- Condição no domicílio (Responsável/Cônjuge, Filho/Enteado/Genro/Nora/Neto/Bisneto, Pai/Mãe/Padrasto/Madrasta/Sogro/Sogra/Avô/Avó)
- Região (SE, NE, N, CO, S)
- Idade (Até 24, 25 a 64, 65 ou mais)
- Renda Domiciliar per capita (Até 218, 218 a 1412, 1 a 2 SM, 2 a 5 SM, Acima de 5 SM)
Categorias de referência:
- Quem está na força de trabalho (na_forca_trabalho == “0”)
- Quem não recebeu bolsa família (recebeu_bf == “0”)
- Homens
- Até 24 anos de idade
- Brancos
- Responsável/cônjuge no domicílio
- Região sudeste
- Renda Domiciliar per capita de até 218
Filtros:
- Foram mantidos os registros de pessoas que faziam parte da força de trabalho ou que faziam parte da força de trabalho potencial.
# Function to prepare data for extensive margin analysis
prepare_extensive_margin_data <- function(data) {
data |>
filter(VD4003 %ni% c("2")) |> # Excluindo pessoas fora da força de trabalho e fora da força de trabalho potencial. Poderia ser outro critério. A ideia é retirar quem não entraria no MT de qualquer jeito (aposentados, por exemplo)
mutate(
na_forca_trabalho = case_when(
VD4001 == "1" ~ "0",
VD4001 == "2" ~ "1",
TRUE ~ NA_character_
),
horas_ef = as.numeric(VD4035),
log_horas_ef = log(horas_ef + 0.01),
rem_ef_todos = as.numeric(VD4020),
log_rem_ef_todos = log(rem_ef_todos + 0.01),
recebeu_bf = case_when(
V5002A == "2" ~ "0",
V5002A == "1" ~ "1",
TRUE ~ NA_character_
),
sexo = as.factor(V2007),
idade = as.numeric(V2009),
log_idade = log(idade + 0.01),
raca = case_when(
V2010 == "1" ~ "B",
V2010 %in% c("2", "4") ~ "PP",
V2010 %in% c("3", "5", "9") ~ "O",
TRUE ~ NA_character_
) |> as.factor(),
anos_estudo = case_when(
VD3006 %in% c("1", "2", "3", "4") ~ "Até 11",
VD3006 %in% c("5") ~ "12 a 15",
VD3006 %in% c("6") ~ "16 ou mais",
TRUE ~ NA_character_
),
cond_dom = case_when(
VD2002 %in% c("01", "02") ~ "Responsável/Cônjuge",
VD2002 %in% c("03", "04", "05", "08", "09") ~ "Filho/Enteado/Genro/Nora/Neto/Bisneto",
VD2002 %in% c("06", "07", "11") ~ "Pai/Mãe/Padrasto/Madrasta/Sogro/Sogra/Avô/Avó",
TRUE ~ NA_character_
),
n_pessoas_dom = as.numeric(VD2003),
log_n_pessoas_dom = log(n_pessoas_dom + 0.01),
renda_dom_pc = as.numeric(VD5005),
log_renda_dom_pc = log(renda_dom_pc + 0.01),
UF = as.factor(UF),
regiao = case_when(
UF %in% c("11", "12", "13", "14", "15", "16", "17") ~ "N",
UF %in% c("21", "22", "23", "24", "25", "26", "27", "28", "29") ~ "NE",
UF %in% c("31", "32", "33", "35") ~ "SE",
UF %in% c("41", "42", "43") ~ "S",
UF %in% c("50", "51", "52", "53") ~ "CO",
TRUE ~ NA_character_
),
area_urbana = as.factor(V1022),
capital_regiao = as.factor(V1023),
pos_ocup = case_when(
VD4009 %in% c("01", "03", "05") ~ "Com carteira",
VD4009 %in% c("02", "04", "06") ~ "Sem carteira",
VD4009 %in% c("09") ~ "Conta própria",
VD4009 %in% c("07", "09", "10") ~ "Outros"
),
grup_ativ = as.factor(VD4010),
grup_ocup = as.factor(VD4011)
) |>
mutate(
na_forca_trabalho = factor(na_forca_trabalho, levels = c("0", "1")),
recebeu_bf = factor(recebeu_bf, levels = c("0", "1")),
anos_estudo = factor(anos_estudo, levels = c("Até 11", "12 a 15", "16 ou mais")),
raca = factor(raca, levels = c("B", "PP", "O")),
cond_dom = factor(cond_dom, levels = c("Responsável/Cônjuge", "Filho/Enteado/Genro/Nora/Neto/Bisneto", "Pai/Mãe/Padrasto/Madrasta/Sogro/Sogra/Avô/Avó")),
regiao = factor(regiao, levels = c("SE", "NE", "N", "CO", "S"))
) |>
mutate(
grupo_idade = cut(idade,
breaks = c(-Inf, 24, 64, Inf),
labels = c("Até 24", "25 a a 64", "65 ou mais"),
include.lowest = TRUE,
right = FALSE
),
grupo_renda_dom_pc = cut(renda_dom_pc,
breaks = c(-Inf, 218, 1412, 7060, Inf),
labels = c("Até R$ 218,00",
"R$ 218,01 a 1,00 SM",
"1,01 SM a 5,00 SM",
"Acima de 5,00 SM"),
include.lowest = TRUE,
right = FALSE
)
)
}
dados_pnadc_filtered_total <- prepare_extensive_margin_data(dados_pnadc)
vars_modelo_total <- c(
"na_forca_trabalho", "recebeu_bf", "sexo", "idade", "log_idade", "grupo_idade", "raca", "cond_dom",
"anos_estudo", "n_pessoas_dom", "log_n_pessoas_dom", "renda_dom_pc", "log_renda_dom_pc", "grupo_renda_dom_pc",
"regiao", "area_urbana", "capital_regiao"
)
dados_pnadc_filtered_total_semp <- dados_pnadc_filtered_total |>
drop_na(all_of(vars_modelo_total)) |>
dplyr::select(all_of(vars_modelo_total))
cat("Amostra para análise da margem extensiva:", nrow(dados_pnadc_filtered_total_semp), "observações\n")## Amostra para análise da margem extensiva: 188910 observaçõesEspecificação do Modelo
# Function to fit the logistic regression model
fit_logit_model <- function(data) {
glm(na_forca_trabalho ~ recebeu_bf + sexo + grupo_idade + raca + anos_estudo +
cond_dom + grupo_renda_dom_pc + regiao + area_urbana, family = binomial("logit"),
data = data)
}
glm_base_total_final <- fit_logit_model(dados_pnadc_filtered_total_semp)
vcov_rob_total_3 <- sandwich::vcovHC(glm_base_total_final, type = "HC1")
# Sumário do modelo
summ(glm_base_total_final, exp = TRUE, confint = TRUE, vifs = TRUE, robust = "HC1")| Observations | 188910 |
| Dependent variable | na_forca_trabalho |
| Type | Generalized linear model |
| Family | binomial |
| Link | logit |
| χ²(18) | 19103.36 |
| p | 0.00 |
| Pseudo-R² (Cragg-Uhler) | 0.24 |
| Pseudo-R² (McFadden) | 0.20 |
| AIC | 76372.42 |
| BIC | 76565.25 |
| exp(Est.) | 2.5% | 97.5% | z val. | p | VIF | |
|---|---|---|---|---|---|---|
| (Intercept) | 0.23 | 0.21 | 0.26 | -27.92 | 0.00 | NA |
| recebeu_bf1 | 2.32 | 2.19 | 2.45 | 30.31 | 0.00 | 1.43 |
| sexo2 | 2.37 | 2.27 | 2.48 | 38.40 | 0.00 | 1.25 |
| grupo_idade25 a a 64 | 0.48 | 0.45 | 0.50 | -27.27 | 0.00 | 1.82 |
| grupo_idade65 ou mais | 1.49 | 1.34 | 1.64 | 7.77 | 0.00 | 1.82 |
| racaPP | 1.07 | 1.02 | 1.12 | 2.77 | 0.01 | 1.14 |
| racaO | 1.58 | 1.33 | 1.86 | 5.36 | 0.00 | 1.14 |
| anos_estudo12 a 15 | 0.60 | 0.57 | 0.63 | -22.68 | 0.00 | 1.31 |
| anos_estudo16 ou mais | 0.43 | 0.39 | 0.47 | -18.21 | 0.00 | 1.31 |
| cond_domFilho/Enteado/Genro/Nora/Neto/Bisneto | 2.14 | 2.03 | 2.25 | 29.11 | 0.00 | 1.77 |
| cond_domPai/Mãe/Padrasto/Madrasta/Sogro/Sogra/Avô/Avó | 0.67 | 0.55 | 0.82 | -3.92 | 0.00 | 1.77 |
| grupo_renda_dom_pcR$ 218,01 a 1,00 SM | 0.25 | 0.23 | 0.26 | -42.00 | 0.00 | 1.36 |
| grupo_renda_dom_pc1,01 SM a 5,00 SM | 0.09 | 0.08 | 0.10 | -59.44 | 0.00 | 1.36 |
| grupo_renda_dom_pcAcima de 5,00 SM | 0.09 | 0.08 | 0.11 | -23.84 | 0.00 | 1.36 |
| regiaoNE | 2.11 | 1.99 | 2.23 | 26.65 | 0.00 | 1.23 |
| regiaoN | 1.28 | 1.19 | 1.37 | 6.88 | 0.00 | 1.23 |
| regiaoCO | 1.13 | 1.03 | 1.23 | 2.71 | 0.01 | 1.23 |
| regiaoS | 0.75 | 0.69 | 0.82 | -6.68 | 0.00 | 1.23 |
| area_urbana2 | 1.48 | 1.42 | 1.54 | 18.72 | 0.00 | 1.11 |
| Standard errors: Robust, type = HC1 |
Diagnósticos e Validação
# Função para diagnósticos completos do modelo logístico
perform_complete_diagnostics_logit <- function(model, data, vcov_model, model_name = "Modelo Logístico") {
cat("===============================================\n")
cat(paste("DIAGNÓSTICOS E VALIDAÇÃO DO", toupper(model_name), "\n"))
cat("===============================================\n\n")
# 1. ANÁLISE DO AJUSTE DO MODELO
cat("1. ANÁLISE DO AJUSTE DO MODELO\n")
cat("------------------------------\n")
# Sumário com R² do jtools
cat("Diferentes tipos de R²:\n")
tryCatch({
# Para modelos logísticos, calcular R² manualmente usando as
# informações do modelo
null_deviance <- model$null.deviance
residual_deviance <- model$deviance
n <- length(model$y)
# McFadden R²
mcfadden_r2 <- 1 - (residual_deviance/null_deviance)
cat("McFadden R²:", round(mcfadden_r2, 4), "\n")
# Cox & Snell R²
cox_snell_r2 <- 1 - exp((residual_deviance - null_deviance)/n)
cat("Cox & Snell R²:", round(cox_snell_r2, 4), "\n")
# Nagelkerke R²
nagelkerke_r2 <- cox_snell_r2/(1 - exp(-null_deviance/n))
cat("Nagelkerke (Cragg-Uhler) R²:", round(nagelkerke_r2, 4), "\n")
# Tentar obter do summ() se disponível
model_summary <- summ(model, exp = FALSE, robust = "HC1")
if (is.list(model_summary) && "coeftable" %in% names(model_summary)) {
cat("Sumário do modelo obtido com sucesso\n")
}
cat("\n")
}, error = function(e) {
cat("Erro ao calcular R²:", e$message, "\n")
cat("Usando informações básicas do modelo:\n")
cat("Null Deviance:", round(model$null.deviance, 2), "\n")
cat("Residual Deviance:", round(model$deviance, 2), "\n")
cat("Diferença:", round(model$null.deviance - model$deviance, 2), "\n\n")
})
# Outros indicadores de ajuste e performance do modelo
cat("Outros indicadores de ajuste e performance do modelo:\n")
model_perf <- model_performance(model)
print(model_perf)
cat("\n")
# 2. HETEROCEDASTICIDADE
cat("2. HETEROCEDASTICIDADE\n")
cat("----------------------\n")
# Breusch-Pagan
cat("Teste de Breusch-Pagan:\n")
bp_test <- bptest(model)
print(bp_test)
cat("Interpretação: p-valor < 0.05 indica presença de heterocedasticidade\n\n")
# Goldfeld-Quandt (usando uma variável contínua como ordenação) Para
# modelos logísticos, valores ajustados
cat("Teste de Goldfeld-Quandt:\n")
tryCatch({
gq_test <- gqtest(model, order.by = fitted(model))
print(gq_test)
cat("Interpretação: p-valor < 0.05 indica presença de heterocedasticidade\n\n")
}, error = function(e) {
cat("Não foi possível executar o teste GQ para modelo logístico\n")
cat("Motivo: Este teste é mais adequado para modelos lineares\n\n")
})
# White Test (usando pacote skedastic)
cat("Teste de White:\n")
tryCatch({
# Teste de White
white_test <- skedastic::white(model)
print(white_test)
cat("Interpretação: p-valor < 0.05 indica presença de heterocedasticidade\n\n")
}, error = function(e) {
tryCatch({
# Fallback: tentar whitestrap
white_test <- whitestrap::white_test(model)
print(white_test)
cat("Interpretação: p-valor < 0.05 indica presença de heterocedasticidade\n\n")
})
})
# 3. MULTICOLINEARIDADE
cat("3. MULTICOLINEARIDADE\n")
cat("---------------------\n")
# VIF usando car
cat("Fatores de Inflação da Variância (VIF):\n")
vif_values <- car::vif(model)
print(vif_values)
cat("\nInterpretação:\n")
cat("- VIF < 5: Não há problema de multicolinearidade\n")
cat("- VIF entre 5-10: Multicolinearidade moderada\n")
cat("- VIF > 10: Multicolinearidade severa\n\n")
# Identificar variáveis problemáticas
if (any(vif_values > 10)) {
cat("ATENÇÃO: Variáveis com VIF > 10 (multicolinearidade severa):\n")
high_vif <- vif_values[vif_values > 10]
print(high_vif)
cat("\n")
}
# 4. SIGNIFICÂNCIA INDIVIDUAL
cat("4. SIGNIFICÂNCIA INDIVIDUAL\n")
cat("---------------------------\n")
# Teste de Wald Tipo II
cat("Teste de Wald (Tipo II):\n")
wald_test_2 <- car::Anova(model, type = "II", test.statistic = "Wald", white.adjust = "hc1")
print(wald_test_2)
cat("Interpretação: p-valor < 0.05 indica que a variável é estatisticamente significativa\n\n")
# Teste LRT Tipo II
cat("Teste de Razão de Verossimilhança (Tipo II):\n")
lr_test_2 <- car::Anova(model, type = "II", test.statistic = "LR", white.adjust = "hc1")
print(lr_test_2)
cat("Interpretação: p-valor < 0.05 indica que a variável é estatisticamente significativa\n\n")
# Teste de Wald Tipo III
cat("Teste de Wald (Tipo III):\n")
wald_test_3 <- car::Anova(model, type = "III", test.statistic = "Wald", white.adjust = "hc1")
print(wald_test_3)
cat("Interpretação: p-valor < 0.05 indica que a variável é estatisticamente significativa\n\n")
# Teste LRT Tipo III
cat("Teste de Razão de Verossimilhança (Tipo III):\n")
lr_test_3 <- car::Anova(model, type = "III", test.statistic = "LR", white.adjust = "hc1")
print(lr_test_3)
cat("Interpretação: p-valor < 0.05 indica que a variável é estatisticamente significativa\n\n")
# 5. SIGNIFICÂNCIA CONJUNTA
cat("5.4 Teste de Wald global\n")
cat("---------------------------------------\n")
tryCatch({
# Teste de Wald para todos os coeficientes exceto intercepto
n_coef <- length(coef(model))
if (n_coef > 1) {
terms_to_test <- 2:n_coef
wald_aod <- wald.test(b = coef(model), Sigma = vcov_model, Terms = terms_to_test)
print(wald_aod)
# Interpretação
p_val <- wald_aod$result$chi2[3]
cat("\nInterpretação:")
if (p_val < 0.001) {
cat(" Modelo altamente significativo (p < 0.001) ***\n")
} else if (p_val < 0.01) {
cat(" Modelo muito significativo (p < 0.01) **\n")
} else if (p_val < 0.05) {
cat(" Modelo significativo (p < 0.05) *\n")
} else {
cat(" Modelo não significativo (p >= 0.05)\n")
}
} else {
cat("Modelo tem apenas intercepto - teste não aplicável\n")
}
cat("\n")
}, error = function(e) {
cat("Erro no teste AOD/Wald:", e$message, "\n\n")
})
# 6. ESPECIFICAÇÃO DO MODELO
cat("6. ESPECIFICAÇÃO DO MODELO\n")
cat("--------------------------\n")
# RESET Test
cat("Teste RESET (Ramsey) - Termos quadráticos e cúbicos:\n")
tryCatch({
reset_test <- resettest(model, power = 2:3, type = "fitted")
print(reset_test)
cat("Interpretação: p-valor < 0.05 indica problemas de especificação funcional\n\n")
}, error = function(e) {
cat("Teste RESET não executado para modelo logístico\n")
cat("Motivo: Este teste é mais adequado para modelos lineares\n\n")
})
# Teste de Hosmer-Lemeshow (específico para modelos logísticos)
cat("Teste de Hosmer-Lemeshow (Goodness-of-Fit para modelos logísticos):\n")
hosmer_test <- performance_hosmer(model)
print(hosmer_test)
cat("Interpretação: p-valor > 0.05 indica bom ajuste do modelo\n\n")
# 7. LINEARIDADE
cat("7. LINEARIDADE (PARA MODELOS LOGÍSTICOS)\n")
cat("----------------------------------------\n")
cat("Para modelos logísticos, a linearidade se refere à relação linear entre\n")
cat("as variáveis explicativas e o logit da probabilidade.\n")
cat("Verificação através dos resíduos de Pearson:\n\n")
# Resíduos de Pearson
pearson_resid <- residuals(model, type = "pearson")
cat("Estatísticas dos Resíduos de Pearson:\n")
cat("Média:", round(mean(pearson_resid), 4), "\n")
cat("Desvio Padrão:", round(sd(pearson_resid), 4), "\n")
cat("Mín:", round(min(pearson_resid), 4), "Máx:", round(max(pearson_resid),
4), "\n\n")
# 8. NORMALIDADE DOS RESÍDUOS
cat("8. NORMALIDADE DOS RESÍDUOS\n")
cat("---------------------------\n")
cat("NOTA: Para modelos logísticos, os resíduos não seguem distribuição normal.\n")
cat("Os testes abaixo são informativos, mas não são requisitos do modelo.\n\n")
# Teste de Shapiro-Wilk (apenas para amostras pequenas)
if (length(pearson_resid) <= 5000) {
cat("Teste de Shapiro-Wilk (Resíduos de Pearson):\n")
shapiro_test <- shapiro.test(pearson_resid)
print(shapiro_test)
cat("\n")
} else {
cat("Teste de Shapiro-Wilk não executado (amostra muito grande)\n\n")
}
# Teste de Jarque-Bera
cat("Teste de Jarque-Bera (Resíduos de Pearson):\n")
jb_test <- jarque.bera.test(pearson_resid)
print(jb_test)
cat("Interpretação: Para modelos logísticos, rejeição da normalidade é esperada\n\n")
# 9. ENDOGENEIDADE (COM DEBUG)
cat("9. ENDOGENEIDADE\n")
cat("-----------------\n")
cat("Análise de correlação entre resíduos e variáveis explicativas\n")
cat("(Teste indireto para detecção de endogeneidade)\n\n")
model_resid <- residuals(model, type = "pearson")
test_vars <- c("recebeu_bf", "sexo", "grupo_idade", "raca", "anos_estudo", "cond_dom",
"grupo_renda_dom_pc", "regiao", "area_urbana")
test_vars <- test_vars[test_vars %in% names(data)]
correlation_results <- data.frame(Variavel = character(), Correlacao = numeric(),
P_valor = numeric(), Significativo = character(), stringsAsFactors = FALSE)
cat("Correlações entre resíduos e variáveis explicativas:\n")
cat("---------------------------------------------------\n")
cat("Resíduos Pearson vs Correlação Pearson:\n")
cat("---------------------------------------------------\n")
for (i in 1:length(test_vars)) {
var <- test_vars[i]
tryCatch({
# Converter variável para numérica se necessário
if (is.factor(data[[var]])) {
var_numeric <- as.numeric(data[[var]])
} else if (is.character(data[[var]])) {
var_numeric <- as.numeric(as.factor(data[[var]]))
} else {
var_numeric <- as.numeric(data[[var]])
}
valid_indices <- !is.na(var_numeric) & !is.na(model_resid)
# Testar se há observações suficientes
if (sum(valid_indices) > 10) {
# Teste de correlação (usando Spearman por robustez)
cor_test <- cor.test(model_resid[valid_indices], var_numeric[valid_indices],
method = "pearson")
# Determinar significância
significativo <- ifelse(cor_test$p.value < 0.05, "Sim", "Não")
# Adicionar aos resultados
correlation_results <- rbind(correlation_results, data.frame(Variavel = var,
Correlacao = round(cor_test$estimate, 4), P_valor = round(cor_test$p.value,
6), Significativo = significativo))
# Mostrar resultado
if (var == "recebeu_bf") {
cat("*** ", var, " (VARIÁVEL PRINCIPAL): ", round(cor_test$estimate,
4), " (p = ", round(cor_test$p.value, 6), ") ", significativo,
" ***\n", sep = "")
} else {
cat(var, ": ", round(cor_test$estimate, 4), " (p = ", round(cor_test$p.value,
4), ") ", significativo, "\n", sep = "")
}
} else {
cat(var, ": Dados insuficientes para teste\n")
}
}, error = function(e) {
cat("ERRO na variável", var, ":", e$message, "\n")
})
}
cat("---------------------------------------------------\n")
cat("Resíduos Pearson vs Correlação Spearman:\n")
cat("---------------------------------------------------\n")
for (i in 1:length(test_vars)) {
var <- test_vars[i]
tryCatch({
# Converter variável para numérica se necessário
if (is.factor(data[[var]])) {
var_numeric <- as.numeric(data[[var]])
} else if (is.character(data[[var]])) {
var_numeric <- as.numeric(as.factor(data[[var]]))
} else {
var_numeric <- as.numeric(data[[var]])
}
valid_indices <- !is.na(var_numeric) & !is.na(model_resid)
# Testar se há observações suficientes
if (sum(valid_indices) > 10) {
# Teste de correlação (usando Spearman por robustez)
cor_test <- cor.test(model_resid[valid_indices], var_numeric[valid_indices],
method = "spearman")
# Determinar significância
significativo <- ifelse(cor_test$p.value < 0.05, "Sim", "Não")
# Adicionar aos resultados
correlation_results <- rbind(correlation_results, data.frame(Variavel = var,
Correlacao = round(cor_test$estimate, 4), P_valor = round(cor_test$p.value,
6), Significativo = significativo))
# Mostrar resultado
if (var == "recebeu_bf") {
cat("*** ", var, " (VARIÁVEL PRINCIPAL): ", round(cor_test$estimate,
4), " (p = ", round(cor_test$p.value, 6), ") ", significativo,
" ***\n", sep = "")
} else {
cat(var, ": ", round(cor_test$estimate, 4), " (p = ", round(cor_test$p.value,
4), ") ", significativo, "\n", sep = "")
}
} else {
cat(var, ": Dados insuficientes para teste\n")
}
}, error = function(e) {
cat("ERRO na variável", var, ":", e$message, "\n")
})
}
model_resid <- residuals(model, type = "deviance")
cat("---------------------------------------------------\n")
cat("Resíduos Deviance vs Correlação Pearson:\n")
cat("---------------------------------------------------\n")
for (i in 1:length(test_vars)) {
var <- test_vars[i]
tryCatch({
# Converter variável para numérica se necessário
if (is.factor(data[[var]])) {
var_numeric <- as.numeric(data[[var]])
} else if (is.character(data[[var]])) {
var_numeric <- as.numeric(as.factor(data[[var]]))
} else {
var_numeric <- as.numeric(data[[var]])
}
valid_indices <- !is.na(var_numeric) & !is.na(model_resid)
# Testar se há observações suficientes
if (sum(valid_indices) > 10) {
# Teste de correlação (usando Spearman por robustez)
cor_test <- cor.test(model_resid[valid_indices], var_numeric[valid_indices],
method = "pearson")
# Determinar significância
significativo <- ifelse(cor_test$p.value < 0.05, "Sim", "Não")
# Adicionar aos resultados
correlation_results <- rbind(correlation_results, data.frame(Variavel = var,
Correlacao = round(cor_test$estimate, 4), P_valor = round(cor_test$p.value,
6), Significativo = significativo))
# Mostrar resultado
if (var == "recebeu_bf") {
cat("*** ", var, " (VARIÁVEL PRINCIPAL): ", round(cor_test$estimate,
4), " (p = ", round(cor_test$p.value, 6), ") ", significativo,
" ***\n", sep = "")
} else {
cat(var, ": ", round(cor_test$estimate, 4), " (p = ", round(cor_test$p.value,
4), ") ", significativo, "\n", sep = "")
}
} else {
cat(var, ": Dados insuficientes para teste\n")
}
}, error = function(e) {
cat("ERRO na variável", var, ":", e$message, "\n")
})
}
cat("---------------------------------------------------\n")
cat("Resíduos Deviance vs Correlação Spearman:\n")
cat("---------------------------------------------------\n")
for (i in 1:length(test_vars)) {
var <- test_vars[i]
tryCatch({
# Converter variável para numérica se necessário
if (is.factor(data[[var]])) {
var_numeric <- as.numeric(data[[var]])
} else if (is.character(data[[var]])) {
var_numeric <- as.numeric(as.factor(data[[var]]))
} else {
var_numeric <- as.numeric(data[[var]])
}
valid_indices <- !is.na(var_numeric) & !is.na(model_resid)
# Testar se há observações suficientes
if (sum(valid_indices) > 10) {
# Teste de correlação (usando Spearman por robustez)
cor_test <- cor.test(model_resid[valid_indices], var_numeric[valid_indices],
method = "spearman")
# Determinar significância
significativo <- ifelse(cor_test$p.value < 0.05, "Sim", "Não")
# Adicionar aos resultados
correlation_results <- rbind(correlation_results, data.frame(Variavel = var,
Correlacao = round(cor_test$estimate, 4), P_valor = round(cor_test$p.value,
6), Significativo = significativo))
# Mostrar resultado
if (var == "recebeu_bf") {
cat("*** ", var, " (VARIÁVEL PRINCIPAL): ", round(cor_test$estimate,
4), " (p = ", round(cor_test$p.value, 6), ") ", significativo,
" ***\n", sep = "")
} else {
cat(var, ": ", round(cor_test$estimate, 4), " (p = ", round(cor_test$p.value,
4), ") ", significativo, "\n", sep = "")
}
} else {
cat(var, ": Dados insuficientes para teste\n")
}
}, error = function(e) {
cat("ERRO na variável", var, ":", e$message, "\n")
})
}
cat("===============================================\n")
cat("FIM DOS DIAGNÓSTICOS\n")
cat("===============================================\n\n")
}
# Executar diagnósticos para o modelo da margem extensiva
perform_complete_diagnostics_logit(glm_base_total_final, dados_pnadc_filtered_total_semp,
vcov_rob_total_3, "Margem Extensiva")## ===============================================
## DIAGNÓSTICOS E VALIDAÇÃO DO MARGEM EXTENSIVA
## ===============================================
##
## 1. ANÁLISE DO AJUSTE DO MODELO
## ------------------------------
## Diferentes tipos de R²:
## McFadden R²: 0.2002
## Cox & Snell R²: 0.0962
## Nagelkerke (Cragg-Uhler) R²: 0.2425
## Sumário do modelo obtido com sucesso
##
## Outros indicadores de ajuste e performance do modelo:
## # Indices of model performance
##
## AIC | AICc | BIC | Tjur's R2 | RMSE | Sigma | Log_loss
## ------------------------------------------------------------------------
## 76372.417 | 76372.421 | 76565.249 | 0.141 | 0.236 | 1.000 | 0.202
##
## AIC | Score_log | Score_spherical | PCP
## -----------------------------------------------
## 76372.417 | -Inf | 1.911e-05 | 0.889
##
## 2. HETEROCEDASTICIDADE
## ----------------------
## Teste de Breusch-Pagan:
##
## studentized Breusch-Pagan test
##
## data: model
## BP = 20701, df = 18, p-value < 2.2e-16
##
## Interpretação: p-valor < 0.05 indica presença de heterocedasticidade
##
## Teste de Goldfeld-Quandt:
##
## Goldfeld-Quandt test
##
## data: model
## GQ = 7.2382, df1 = 94436, df2 = 94436, p-value < 2.2e-16
## alternative hypothesis: variance increases from segment 1 to 2
##
## Interpretação: p-valor < 0.05 indica presença de heterocedasticidade
##
## Teste de White:
## White's test results
##
## Null hypothesis: Homoskedasticity of the residuals
## Alternative hypothesis: Heteroskedasticity of the residuals
## Test Statistic: 193.34
## P-value: 0
## Interpretação: p-valor < 0.05 indica presença de heterocedasticidade
##
## 3. MULTICOLINEARIDADE
## ---------------------
## Fatores de Inflação da Variância (VIF):
## GVIF Df GVIF^(1/(2*Df))
## recebeu_bf 1.431544 1 1.196471
## sexo 1.250075 1 1.118068
## grupo_idade 1.824660 2 1.162239
## raca 1.141736 2 1.033693
## anos_estudo 1.311939 2 1.070233
## cond_dom 1.768527 2 1.153195
## grupo_renda_dom_pc 1.359179 3 1.052477
## regiao 1.233249 4 1.026553
## area_urbana 1.109642 1 1.053395
##
## Interpretação:
## - VIF < 5: Não há problema de multicolinearidade
## - VIF entre 5-10: Multicolinearidade moderada
## - VIF > 10: Multicolinearidade severa
##
## 4. SIGNIFICÂNCIA INDIVIDUAL
## ---------------------------
## Teste de Wald (Tipo II):
## Analysis of Deviance Table (Type II tests)
##
## Response: na_forca_trabalho
## Df Chisq Pr(>Chisq)
## recebeu_bf 1 950.776 < 2.2e-16 ***
## sexo 1 1514.762 < 2.2e-16 ***
## grupo_idade 2 1316.813 < 2.2e-16 ***
## raca 2 32.147 1.046e-07 ***
## anos_estudo 2 718.963 < 2.2e-16 ***
## cond_dom 2 812.669 < 2.2e-16 ***
## grupo_renda_dom_pc 3 3676.981 < 2.2e-16 ***
## regiao 4 1256.953 < 2.2e-16 ***
## area_urbana 1 348.523 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Interpretação: p-valor < 0.05 indica que a variável é estatisticamente significativa
##
## Teste de Razão de Verossimilhança (Tipo II):
## Analysis of Deviance Table (Type II tests)
##
## Response: na_forca_trabalho
## LR Chisq Df Pr(>Chisq)
## recebeu_bf 924.9 1 < 2.2e-16 ***
## sexo 1540.8 1 < 2.2e-16 ***
## grupo_idade 1224.8 2 < 2.2e-16 ***
## raca 30.3 2 2.618e-07 ***
## anos_estudo 737.6 2 < 2.2e-16 ***
## cond_dom 787.2 2 < 2.2e-16 ***
## grupo_renda_dom_pc 3577.6 3 < 2.2e-16 ***
## regiao 1287.5 4 < 2.2e-16 ***
## area_urbana 343.0 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Interpretação: p-valor < 0.05 indica que a variável é estatisticamente significativa
##
## Teste de Wald (Tipo III):
## Analysis of Deviance Table (Type III tests)
##
## Response: na_forca_trabalho
## Df Chisq Pr(>Chisq)
## (Intercept) 1 808.053 < 2.2e-16 ***
## recebeu_bf 1 950.776 < 2.2e-16 ***
## sexo 1 1514.762 < 2.2e-16 ***
## grupo_idade 2 1316.813 < 2.2e-16 ***
## raca 2 32.147 1.046e-07 ***
## anos_estudo 2 718.963 < 2.2e-16 ***
## cond_dom 2 812.669 < 2.2e-16 ***
## grupo_renda_dom_pc 3 3676.981 < 2.2e-16 ***
## regiao 4 1256.953 < 2.2e-16 ***
## area_urbana 1 348.523 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Interpretação: p-valor < 0.05 indica que a variável é estatisticamente significativa
##
## Teste de Razão de Verossimilhança (Tipo III):
## Analysis of Deviance Table (Type III tests)
##
## Response: na_forca_trabalho
## LR Chisq Df Pr(>Chisq)
## recebeu_bf 924.9 1 < 2.2e-16 ***
## sexo 1540.8 1 < 2.2e-16 ***
## grupo_idade 1224.8 2 < 2.2e-16 ***
## raca 30.3 2 2.618e-07 ***
## anos_estudo 737.6 2 < 2.2e-16 ***
## cond_dom 787.2 2 < 2.2e-16 ***
## grupo_renda_dom_pc 3577.6 3 < 2.2e-16 ***
## regiao 1287.5 4 < 2.2e-16 ***
## area_urbana 343.0 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Interpretação: p-valor < 0.05 indica que a variável é estatisticamente significativa
##
## 5.4 Teste de Wald global
## ---------------------------------------
## Wald test:
## ----------
##
## Chi-squared test:
## X2 = 15306.7, df = 18, P(> X2) = 0.0
##
## Interpretação: Modelo altamente significativo (p < 0.001) ***
##
## 6. ESPECIFICAÇÃO DO MODELO
## --------------------------
## Teste RESET (Ramsey) - Termos quadráticos e cúbicos:
##
## RESET test
##
## data: model
## RESET = 2184.5, df1 = 2, df2 = 188889, p-value < 2.2e-16
##
## Interpretação: p-valor < 0.05 indica problemas de especificação funcional
##
## Teste de Hosmer-Lemeshow (Goodness-of-Fit para modelos logísticos):
## # Hosmer-Lemeshow Goodness-of-Fit Test
##
## Chi-squared: 36.810
## df: 8
## p-value: 0.000## Interpretação: p-valor > 0.05 indica bom ajuste do modelo
##
## 7. LINEARIDADE (PARA MODELOS LOGÍSTICOS)
## ----------------------------------------
## Para modelos logísticos, a linearidade se refere à relação linear entre
## as variáveis explicativas e o logit da probabilidade.
## Verificação através dos resíduos de Pearson:
##
## Estatísticas dos Resíduos de Pearson:
## Média: -0.0022
## Desvio Padrão: 0.9893
## Mín: -2.4135 Máx: 17.5542
##
## 8. NORMALIDADE DOS RESÍDUOS
## ---------------------------
## NOTA: Para modelos logísticos, os resíduos não seguem distribuição normal.
## Os testes abaixo são informativos, mas não são requisitos do modelo.
##
## Teste de Shapiro-Wilk não executado (amostra muito grande)
##
## Teste de Jarque-Bera (Resíduos de Pearson):
##
## Jarque Bera Test
##
## data: pearson_resid
## X-squared = 21129157, df = 2, p-value < 2.2e-16
##
## Interpretação: Para modelos logísticos, rejeição da normalidade é esperada
##
## 9. ENDOGENEIDADE
## -----------------
## Análise de correlação entre resíduos e variáveis explicativas
## (Teste indireto para detecção de endogeneidade)
##
## Correlações entre resíduos e variáveis explicativas:
## ---------------------------------------------------
## Resíduos Pearson vs Correlação Pearson:
## ---------------------------------------------------
## *** recebeu_bf (VARIÁVEL PRINCIPAL): 0.0025 (p = 0.281406) Não ***
## sexo: 0.0062 (p = 0.0069) Sim
## grupo_idade: 0.0057 (p = 0.0134) Sim
## raca: -0.0043 (p = 0.0611) Não
## anos_estudo: 0.0044 (p = 0.0542) Não
## cond_dom: 0.0027 (p = 0.2361) Não
## grupo_renda_dom_pc: -0.0011 (p = 0.6191) Não
## regiao: 0.002 (p = 0.3948) Não
## area_urbana: -0.0088 (p = 1e-04) Sim
## ---------------------------------------------------
## Resíduos Pearson vs Correlação Spearman:
## ---------------------------------------------------## *** recebeu_bf (VARIÁVEL PRINCIPAL): -0.2296 (p = 0) Sim ***## sexo: -0.1843 (p = 0) Sim## grupo_idade: 0.1353 (p = 0) Sim## raca: -0.2151 (p = 0) Sim## anos_estudo: 0.3496 (p = 0) Sim## cond_dom: -0.1902 (p = 0) Sim## grupo_renda_dom_pc: 0.512 (p = 0) Sim## regiao: 0.0916 (p = 0) Sim## area_urbana: -0.2268 (p = 0) Sim
## ---------------------------------------------------
## Resíduos Deviance vs Correlação Pearson:
## ---------------------------------------------------
## *** recebeu_bf (VARIÁVEL PRINCIPAL): -0.0022 (p = 0.331426) Não ***
## sexo: -0.0059 (p = 0.0106) Sim
## grupo_idade: 0.0016 (p = 0.4928) Não
## raca: -0.0126 (p = 0) Sim
## anos_estudo: 0.0215 (p = 0) Sim
## cond_dom: -0.0074 (p = 0.0012) Sim
## grupo_renda_dom_pc: 0.0222 (p = 0) Sim
## regiao: 0.0091 (p = 1e-04) Sim
## area_urbana: -0.0104 (p = 0) Sim
## ---------------------------------------------------
## Resíduos Deviance vs Correlação Spearman:
## ---------------------------------------------------## *** recebeu_bf (VARIÁVEL PRINCIPAL): -0.2296 (p = 0) Sim ***## sexo: -0.1843 (p = 0) Sim## grupo_idade: 0.1353 (p = 0) Sim## raca: -0.2151 (p = 0) Sim## anos_estudo: 0.3496 (p = 0) Sim## cond_dom: -0.1902 (p = 0) Sim## grupo_renda_dom_pc: 0.512 (p = 0) Sim## regiao: 0.0916 (p = 0) Sim## area_urbana: -0.2268 (p = 0) Sim
## ===============================================
## FIM DOS DIAGNÓSTICOS
## ===============================================Efeitos Marginais
Sexo
cat("# EFEITOS MARGINAIS - MARGEM EXTENSIVA\n")## # EFEITOS MARGINAIS - MARGEM EXTENSIVAcat("======================================\n")## ======================================# Por sexo
result_sexo_ext <- analyze_and_plot_marginal_effects(glm_base_total_final, "recebeu_bf",
"sexo", vcov_rob_total_3, "Probabilidade Predita", "por Sexo", comparison_type = "comparisons")##
## === EFEITO DO BOLSA FAMÍLIA POR SEXO ===
##
## --- EFEITOS MARGINAIS ---
##
## sexo Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## 1 0.0465 0.00205 22.7 <0.001 376.0 0.0425 0.0505
## 2 0.0693 0.00257 27.0 <0.001 532.2 0.0643 0.0744
##
## Term: recebeu_bf
## Type: response
## Comparison: 1 - 0
##
##
## --- VALORES PREDITOS ---
##
## recebeu_bf sexo Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## 0 1 0.0459 0.000605 75.9 <0.001 Inf 0.0447 0.0471
## 0 2 0.0692 0.000869 79.6 <0.001 Inf 0.0675 0.0709
## 1 1 0.1305 0.002939 44.4 <0.001 Inf 0.1247 0.1363
## 1 2 0.2251 0.003149 71.5 <0.001 Inf 0.2190 0.2313
##
## Type: response
##
##
## --- VERIFICAÇÃO DE CONSISTÊNCIA ---
## Comparação entre diferenças calculadas e efeitos marginais:
## sexo diff_calculada estimate diferenca_absoluta
## 1 1 NA 0.04647808 NA
## 2 1 NA 0.04647808 NA
## 3 2 NA 0.06933944 NA
## 4 2 NA 0.06933944 NA
## Teste Wald (robusto) de diferenças entre sexo:
##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b1=b2 -0.0229 0.000896 -25.5 <0.001 474.9 -0.0246 -0.0211
## b1=0 0.0465 0.002049 22.7 <0.001 376.0 0.0425 0.0505
## b2=0 0.0693 0.002565 27.0 <0.001 532.2 0.0643 0.0744Idade
# Por idade
result_idade_ext <- analyze_and_plot_marginal_effects(glm_base_total_final, "recebeu_bf",
"grupo_idade", vcov_rob_total_3, "Probabilidade Predita", "por Quartil de Idade",
comparison_type = "comparisons")##
## === EFEITO DO BOLSA FAMÍLIA POR GRUPO_IDADE ===
##
## --- EFEITOS MARGINAIS ---
##
## grupo_idade Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## Até 24 0.1126 0.00452 24.9 <0.001 451.6 0.1037 0.1214
## 25 a a 64 0.0463 0.00183 25.3 <0.001 465.4 0.0427 0.0499
## 65 ou mais 0.0736 0.00378 19.5 <0.001 278.6 0.0662 0.0810
##
## Term: recebeu_bf
## Type: response
## Comparison: 1 - 0
##
##
## --- VALORES PREDITOS ---
##
## recebeu_bf grupo_idade Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## 0 Até 24 0.1470 0.002047 71.8 <0.001 Inf 0.1430 0.1510
## 0 25 a a 64 0.0378 0.000479 79.0 <0.001 Inf 0.0369 0.0387
## 0 65 ou mais 0.0742 0.002618 28.4 <0.001 585.0 0.0691 0.0794
## 1 Até 24 0.3387 0.006059 55.9 <0.001 Inf 0.3268 0.3505
## 1 25 a a 64 0.2000 0.002914 68.6 <0.001 Inf 0.1943 0.2058
## 1 65 ou mais 0.3280 0.009431 34.8 <0.001 878.2 0.3095 0.3465
##
## Type: response
##
##
## --- VERIFICAÇÃO DE CONSISTÊNCIA ---
## Comparação entre diferenças calculadas e efeitos marginais:
## grupo_idade diff_calculada estimate diferenca_absoluta
## 1 25 a a 64 NA 0.04628384 NA
## 2 25 a a 64 NA 0.04628384 NA
## 3 65 ou mais NA 0.07359345 NA
## 4 65 ou mais NA 0.07359345 NA
## 5 Até 24 NA 0.11257576 NA
## 6 Até 24 NA 0.11257576 NA
## Teste Wald (robusto) de diferenças entre grupo_idade:
##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b1=b2 0.0663 0.00283 23.4 <0.001 400.7 0.0607 0.0718
## b1=b3 0.0390 0.00250 15.6 <0.001 179.7 0.0341 0.0439
## b2=b3 -0.0273 0.00247 -11.0 <0.001 91.8 -0.0322 -0.0225
## b1=0 0.1126 0.00452 24.9 <0.001 451.6 0.1037 0.1214
## b2=0 0.0463 0.00183 25.3 <0.001 465.4 0.0427 0.0499
## b3=0 0.0736 0.00378 19.5 <0.001 278.6 0.0662 0.0810Raça
# Por raça
result_raca_ext <- analyze_and_plot_marginal_effects(glm_base_total_final, "recebeu_bf",
"raca", vcov_rob_total_3, "Probabilidade Predita", "por Raça", comparison_type = "comparisons")##
## === EFEITO DO BOLSA FAMÍLIA POR RACA ===
##
## --- EFEITOS MARGINAIS ---
##
## raca Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## B 0.0408 0.00179 22.8 <0.001 378.5 0.0373 0.0443
## PP 0.0670 0.00260 25.8 <0.001 483.5 0.0619 0.0721
## O 0.0722 0.00444 16.3 <0.001 195.2 0.0635 0.0809
##
## Term: recebeu_bf
## Type: response
## Comparison: 1 - 0
##
##
## --- VALORES PREDITOS ---
##
## recebeu_bf raca Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## 0 B 0.0376 0.000638 58.8 <0.001 Inf 0.0363 0.0388
## 0 PP 0.0678 0.000724 93.6 <0.001 Inf 0.0664 0.0692
## 0 O 0.0732 0.004883 15.0 <0.001 166.3 0.0636 0.0828
## 1 B 0.1776 0.003584 49.6 <0.001 Inf 0.1706 0.1847
## 1 PP 0.2193 0.003141 69.8 <0.001 Inf 0.2132 0.2255
## 1 O 0.2592 0.014602 17.8 <0.001 231.8 0.2306 0.2878
##
## Type: response
##
##
## --- VERIFICAÇÃO DE CONSISTÊNCIA ---
## Comparação entre diferenças calculadas e efeitos marginais:
## raca diff_calculada estimate diferenca_absoluta
## 1 B NA 0.04078503 NA
## 2 B NA 0.04078503 NA
## 3 O NA 0.07220514 NA
## 4 O NA 0.07220514 NA
## 5 PP NA 0.06696783 NA
## 6 PP NA 0.06696783 NA
## Teste Wald (robusto) de diferenças entre raca:
##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b1=b2 -0.02618 0.00113 -23.20 <0.001 393.1 -0.0284 -0.02397
## b1=b3 -0.03142 0.00370 -8.50 <0.001 55.5 -0.0387 -0.02417
## b2=b3 -0.00524 0.00355 -1.48 0.14 2.8 -0.0122 0.00172
## b1=0 0.04079 0.00179 22.76 <0.001 378.5 0.0373 0.04430
## b2=0 0.06697 0.00260 25.76 <0.001 483.5 0.0619 0.07206
## b3=0 0.07221 0.00444 16.26 <0.001 195.2 0.0635 0.08091Anos de estudo
# Por anos_estudo
result_anos_estudo_ext <- analyze_and_plot_marginal_effects(glm_base_total_final,
"recebeu_bf", "anos_estudo", vcov_rob_total_3, "Probabilidade Predita", "por Faixa de Anos de Estudos",
comparison_type = "comparisons")##
## === EFEITO DO BOLSA FAMÍLIA POR ANOS_ESTUDO ===
##
## --- EFEITOS MARGINAIS ---
##
## anos_estudo Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## Até 11 0.0792 0.00304 26.0 <0.001 493.2 0.0732 0.0852
## 12 a 15 0.0499 0.00207 24.0 <0.001 421.9 0.0458 0.0539
## 16 ou mais 0.0227 0.00130 17.5 <0.001 224.8 0.0201 0.0252
##
## Term: recebeu_bf
## Type: response
## Comparison: 1 - 0
##
##
## --- VALORES PREDITOS ---
##
## recebeu_bf anos_estudo Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## 0 Até 11 0.0830 0.000964 86.0 <0.001 Inf 0.0811 0.0848
## 0 12 a 15 0.0474 0.000690 68.8 <0.001 Inf 0.0461 0.0488
## 0 16 ou mais 0.0188 0.000701 26.8 <0.001 523.0 0.0174 0.0201
## 1 Até 11 0.2463 0.003488 70.6 <0.001 Inf 0.2395 0.2531
## 1 12 a 15 0.1639 0.003167 51.8 <0.001 Inf 0.1577 0.1701
## 1 16 ou mais 0.1131 0.004539 24.9 <0.001 452.6 0.1042 0.1220
##
## Type: response
##
##
## --- VERIFICAÇÃO DE CONSISTÊNCIA ---
## Comparação entre diferenças calculadas e efeitos marginais:
## anos_estudo diff_calculada estimate diferenca_absoluta
## 1 12 a 15 NA 0.04986114 NA
## 2 12 a 15 NA 0.04986114 NA
## 3 16 ou mais NA 0.02269068 NA
## 4 16 ou mais NA 0.02269068 NA
## 5 Até 11 NA 0.07920571 NA
## 6 Até 11 NA 0.07920571 NA
## Teste Wald (robusto) de diferenças entre anos_estudo:
##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b1=b2 0.0293 0.00129 22.8 <0.001 379.9 0.0268 0.0319
## b1=b3 0.0565 0.00217 26.1 <0.001 494.6 0.0523 0.0608
## b2=b3 0.0272 0.00135 20.2 <0.001 298.4 0.0245 0.0298
## b1=0 0.0792 0.00304 26.0 <0.001 493.2 0.0732 0.0852
## b2=0 0.0499 0.00207 24.0 <0.001 421.9 0.0458 0.0539
## b3=0 0.0227 0.00130 17.5 <0.001 224.8 0.0201 0.0252Condição no domicílio
# Por cond_dom
result_cond_dom_ext <- analyze_and_plot_marginal_effects(glm_base_total_final, "recebeu_bf",
"cond_dom", vcov_rob_total_3, "Probabilidade Predita", "por Condição no Domicílio",
comparison_type = "comparisons")##
## === EFEITO DO BOLSA FAMÍLIA POR COND_DOM ===
##
## --- EFEITOS MARGINAIS ---
##
## cond_dom Estimate Std. Error z
## Responsável/Cônjuge 0.0472 0.00185 25.5
## Filho/Enteado/Genro/Nora/Neto/Bisneto 0.0900 0.00375 24.0
## Pai/Mãe/Padrasto/Madrasta/Sogro/Sogra/Avô/Avó 0.0371 0.00335 11.1
## Pr(>|z|) S 2.5 % 97.5 %
## <0.001 474.0 0.0436 0.0508
## <0.001 421.6 0.0827 0.0974
## <0.001 92.3 0.0305 0.0436
##
## Term: recebeu_bf
## Type: response
## Comparison: 1 - 0
##
##
## --- VALORES PREDITOS ---
##
## recebeu_bf cond_dom Estimate Std. Error
## 0 Responsável/Cônjuge 0.0384 0.000505
## 0 Filho/Enteado/Genro/Nora/Neto/Bisneto 0.1091 0.001392
## 0 Pai/Mãe/Padrasto/Madrasta/Sogro/Sogra/Avô/Avó 0.0315 0.002953
## 1 Responsável/Cônjuge 0.2050 0.002950
## 1 Filho/Enteado/Genro/Nora/Neto/Bisneto 0.3063 0.005530
## 1 Pai/Mãe/Padrasto/Madrasta/Sogro/Sogra/Avô/Avó 0.1369 0.011200
## z Pr(>|z|) S 2.5 % 97.5 %
## 76.0 <0.001 Inf 0.0374 0.0394
## 78.4 <0.001 Inf 0.1064 0.1119
## 10.7 <0.001 86.0 0.0257 0.0373
## 69.5 <0.001 Inf 0.1993 0.2108
## 55.4 <0.001 Inf 0.2955 0.3172
## 12.2 <0.001 111.8 0.1150 0.1589
##
## Type: response
##
##
## --- VERIFICAÇÃO DE CONSISTÊNCIA ---
## Comparação entre diferenças calculadas e efeitos marginais:
## cond_dom diff_calculada estimate
## 1 Filho/Enteado/Genro/Nora/Neto/Bisneto NA 0.09004373
## 2 Filho/Enteado/Genro/Nora/Neto/Bisneto NA 0.09004373
## 3 Pai/Mãe/Padrasto/Madrasta/Sogro/Sogra/Avô/Avó NA 0.03706924
## 4 Pai/Mãe/Padrasto/Madrasta/Sogro/Sogra/Avô/Avó NA 0.03706924
## 5 Responsável/Cônjuge NA 0.04720196
## 6 Responsável/Cônjuge NA 0.04720196
## diferenca_absoluta
## 1 NA
## 2 NA
## 3 NA
## 4 NA
## 5 NA
## 6 NA
## Teste Wald (robusto) de diferenças entre cond_dom:
##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b1=b2 -0.0428 0.00206 -20.83 <0.001 317.8 -0.04687 -0.0388
## b1=b3 0.0101 0.00301 3.37 <0.001 10.4 0.00424 0.0160
## b2=b3 0.0530 0.00372 14.22 <0.001 150.1 0.04567 0.0603
## b1=0 0.0472 0.00185 25.50 <0.001 474.0 0.04357 0.0508
## b2=0 0.0900 0.00375 24.03 <0.001 421.6 0.08270 0.0974
## b3=0 0.0371 0.00335 11.08 <0.001 92.3 0.03051 0.0436Renda Dom. P.C.
# Por renda
result_renda_ext <- analyze_and_plot_marginal_effects(glm_base_total_final, "recebeu_bf",
"grupo_renda_dom_pc", vcov_rob_total_3, "Probabilidade Predita", "por Quartil de Renda per Capita",
comparison_type = "comparisons")##
## === EFEITO DO BOLSA FAMÍLIA POR GRUPO_RENDA_DOM_PC ===
##
## --- EFEITOS MARGINAIS ---
##
## grupo_renda_dom_pc Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## Até R$ 218,00 0.1668 0.00582 28.7 <0.001 597.6 0.1554 0.1782
## R$ 218,01 a 1,00 SM 0.0847 0.00324 26.1 <0.001 497.9 0.0783 0.0910
## 1,01 SM a 5,00 SM 0.0268 0.00132 20.3 <0.001 301.3 0.0242 0.0294
## Acima de 5,00 SM 0.0176 0.00166 10.6 <0.001 84.5 0.0143 0.0208
##
## Term: recebeu_bf
## Type: response
## Comparison: 1 - 0
##
##
## --- VALORES PREDITOS ---
##
## recebeu_bf grupo_renda_dom_pc Estimate Std. Error z Pr(>|z|) S 2.5 %
## 0 Até R$ 218,00 0.2812 0.005423 51.9 <0.001 Inf 0.2705
## 0 R$ 218,01 a 1,00 SM 0.0858 0.000983 87.3 <0.001 Inf 0.0839
## 0 1,01 SM a 5,00 SM 0.0222 0.000478 46.5 <0.001 Inf 0.0213
## 0 Acima de 5,00 SM 0.0142 0.001230 11.6 <0.001 100.1 0.0118
## 1 Até R$ 218,00 0.5386 0.008309 64.8 <0.001 Inf 0.5224
## 1 R$ 218,01 a 1,00 SM 0.2039 0.003120 65.3 <0.001 Inf 0.1978
## 1 1,01 SM a 5,00 SM 0.0725 0.002088 34.7 <0.001 876.2 0.0684
## 1 Acima de 5,00 SM 0.0618 0.005511 11.2 <0.001 94.5 0.0510
## 97.5 %
## 0.2918
## 0.0877
## 0.0232
## 0.0166
## 0.5549
## 0.2100
## 0.0766
## 0.0726
##
## Type: response
##
##
## --- VERIFICAÇÃO DE CONSISTÊNCIA ---
## Comparação entre diferenças calculadas e efeitos marginais:
## grupo_renda_dom_pc diff_calculada estimate diferenca_absoluta
## 1 1,01 SM a 5,00 SM NA 0.02677569 NA
## 2 1,01 SM a 5,00 SM NA 0.02677569 NA
## 3 Acima de 5,00 SM NA 0.01757338 NA
## 4 Acima de 5,00 SM NA 0.01757338 NA
## 5 Até R$ 218,00 NA 0.16680379 NA
## 6 Até R$ 218,00 NA 0.16680379 NA
## 7 R$ 218,01 a 1,00 SM NA 0.08468449 NA
## 8 R$ 218,01 a 1,00 SM NA 0.08468449 NA
## Teste Wald (robusto) de diferenças entre grupo_renda_dom_pc:
##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b1=b2 0.0821 0.00283 29.06 <0.001 614.3 0.07658 0.0877
## b1=b3 0.1400 0.00466 30.02 <0.001 655.3 0.13089 0.1492
## b1=b4 0.1492 0.00520 28.69 <0.001 598.9 0.13904 0.1594
## b2=b3 0.0579 0.00213 27.20 <0.001 538.9 0.05374 0.0621
## b2=b4 0.0671 0.00281 23.85 <0.001 415.1 0.06160 0.0726
## b3=b4 0.0092 0.00157 5.87 <0.001 27.7 0.00613 0.0123
## b1=0 0.1668 0.00582 28.66 <0.001 597.6 0.15540 0.1782
## b2=0 0.0847 0.00324 26.14 <0.001 497.9 0.07833 0.0910
## b3=0 0.0268 0.00132 20.28 <0.001 301.3 0.02419 0.0294
## b4=0 0.0176 0.00166 10.58 <0.001 84.5 0.01432 0.0208Região
# Por região
result_regiao_ext <- analyze_and_plot_marginal_effects(glm_base_total_final, "recebeu_bf",
"regiao", vcov_rob_total_3, "Probabilidade Predita", "por Região", comparison_type = "comparisons")##
## === EFEITO DO BOLSA FAMÍLIA POR REGIAO ===
##
## --- EFEITOS MARGINAIS ---
##
## regiao Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## SE 0.0385 0.00177 21.8 <0.001 347.1 0.0350 0.0419
## NE 0.0933 0.00351 26.6 <0.001 515.6 0.0864 0.1001
## N 0.0623 0.00263 23.7 <0.001 409.8 0.0572 0.0675
## CO 0.0417 0.00209 20.0 <0.001 292.2 0.0376 0.0458
## S 0.0283 0.00149 19.0 <0.001 265.1 0.0254 0.0312
##
## Term: recebeu_bf
## Type: response
## Comparison: 1 - 0
##
##
## --- VALORES PREDITOS ---
##
## recebeu_bf regiao Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## 0 SE 0.0338 0.000741 45.6 <0.001 Inf 0.0323 0.0353
## 0 NE 0.1037 0.001214 85.4 <0.001 Inf 0.1013 0.1060
## 0 N 0.0614 0.001377 44.6 <0.001 Inf 0.0587 0.0641
## 0 CO 0.0373 0.001226 30.5 <0.001 674.5 0.0349 0.0397
## 0 S 0.0241 0.000780 30.9 <0.001 694.6 0.0226 0.0256
## 1 SE 0.1290 0.003248 39.7 <0.001 Inf 0.1226 0.1354
## 1 NE 0.2678 0.003773 71.0 <0.001 Inf 0.2604 0.2752
## 1 N 0.1779 0.004204 42.3 <0.001 Inf 0.1697 0.1862
## 1 CO 0.1370 0.004600 29.8 <0.001 645.6 0.1280 0.1461
## 1 S 0.1054 0.003655 28.8 <0.001 605.2 0.0982 0.1126
##
## Type: response
##
##
## --- VERIFICAÇÃO DE CONSISTÊNCIA ---
## Comparação entre diferenças calculadas e efeitos marginais:
## regiao diff_calculada estimate diferenca_absoluta
## 1 CO NA 0.04168523 NA
## 2 CO NA 0.04168523 NA
## 3 N NA 0.06234980 NA
## 4 N NA 0.06234980 NA
## 5 NE NA 0.09325704 NA
## 6 NE NA 0.09325704 NA
## 7 S NA 0.02827452 NA
## 8 S NA 0.02827452 NA
## 9 SE NA 0.03848270 NA
## 10 SE NA 0.03848270 NA
## Teste Wald (robusto) de diferenças entre regiao:
##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b1=b2 -0.0548 0.00207 -26.45 <0.001 509.9 -0.05883 -0.050716
## b1=b3 -0.0239 0.00151 -15.84 <0.001 185.3 -0.02682 -0.020914
## b1=b4 -0.0032 0.00136 -2.36 0.0181 5.8 -0.00586 -0.000546
## b1=b5 0.0102 0.00114 8.98 <0.001 61.7 0.00798 0.012436
## b2=b3 0.0309 0.00163 18.91 <0.001 262.5 0.02770 0.034111
## b2=b4 0.0516 0.00216 23.90 <0.001 416.8 0.04734 0.055802
## b2=b5 0.0650 0.00243 26.74 <0.001 520.9 0.06022 0.069745
## b3=b4 0.0207 0.00170 12.13 <0.001 110.0 0.01732 0.024004
## b3=b5 0.0341 0.00177 19.27 <0.001 272.3 0.03061 0.037542
## b4=b5 0.0134 0.00149 9.01 <0.001 62.1 0.01049 0.016328
## b1=0 0.0385 0.00177 21.78 <0.001 347.1 0.03502 0.041945
## b2=0 0.0933 0.00351 26.60 <0.001 515.6 0.08639 0.100127
## b3=0 0.0623 0.00263 23.69 <0.001 409.8 0.05719 0.067508
## b4=0 0.0417 0.00209 19.97 <0.001 292.2 0.03759 0.045777
## b5=0 0.0283 0.00149 19.01 <0.001 265.1 0.02536 0.031190Urbano ou Rural
# Por área urbana
result_area_ext <- analyze_and_plot_marginal_effects(glm_base_total_final, "recebeu_bf",
"area_urbana", vcov_rob_total_3, "Probabilidade Predita", "por Área Urbana",
comparison_type = "comparisons")##
## === EFEITO DO BOLSA FAMÍLIA POR AREA_URBANA ===
##
## --- EFEITOS MARGINAIS ---
##
## area_urbana Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## 1 0.0468 0.00194 24.1 <0.001 424.7 0.0430 0.0506
## 2 0.0877 0.00332 26.4 <0.001 506.8 0.0811 0.0942
##
## Term: recebeu_bf
## Type: response
## Comparison: 1 - 0
##
##
## --- VALORES PREDITOS ---
##
## recebeu_bf area_urbana Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## 0 1 0.0431 0.000513 84.1 <0.001 Inf 0.0421 0.0441
## 0 2 0.0969 0.001249 77.6 <0.001 Inf 0.0945 0.0994
## 1 1 0.1721 0.002885 59.6 <0.001 Inf 0.1664 0.1777
## 1 2 0.2715 0.004034 67.3 <0.001 Inf 0.2635 0.2794
##
## Type: response
##
##
## --- VERIFICAÇÃO DE CONSISTÊNCIA ---
## Comparação entre diferenças calculadas e efeitos marginais:
## area_urbana diff_calculada estimate diferenca_absoluta
## 1 1 NA 0.04678864 NA
## 2 1 NA 0.04678864 NA
## 3 2 NA 0.08765636 NA
## 4 2 NA 0.08765636 NA
## Teste Wald (robusto) de diferenças entre area_urbana:
##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b1=b2 -0.0409 0.00159 -25.7 <0.001 480.1 -0.0440 -0.0377
## b1=0 0.0468 0.00194 24.1 <0.001 424.7 0.0430 0.0506
## b2=0 0.0877 0.00332 26.4 <0.001 506.8 0.0811 0.0942Margem Intensiva
Preparação dos Dados
Variável Dependente:
- Horas efetivas de trabalho (VD4035, horas_ef)
Variável Explicativa Principal:
- Recebimento do Bolsa Família (V5002A, recebeu_bf)
Variáveis de Controle:
- Características individuais: sexo (V2007), idade (V2009 agrupada por quartil), raça (V2010 agrupada), escolaridade (VD3006 agrupada, faixa de anos de estudos)
- Características domiciliares: condição no domicílio (VD2002 agrupada), renda domiciliar per capita (VD5005 agrupada por quartil, efetivo, inclusive rendimentos de tíquete transporte ou alimentação, exclusive pessoas do domicílio cuja condição era pensionista, doméstico ou parente de doméstico)
- Características geográficas: região (UF agrupada), área urbana (V1022)
- Características ocupacionais: posição na ocupação (VD4009 agrupada)
Variáveis de Interação:
- Recebimento do Bolsa Família vs Área Urbana (efeito do recebimento de PBF nas horas trabalhadas é diferente entre rural vs urbano)
- Recebimento do Bolsa Família vs Posição na Ocupação (efeito do recebimento de PBF nas horas trabalhadas é diferente entre com carteira, sem carteira e conta-própria)
Agrupamentos realizados nas variáveis:
- Anos Estudo (Até 11, 12 a 15, 16 ou mais)
- Raça (Brancos, Pretos e Pardos, Outros)
- Condição no domicílio (Responsável/Cônjuge, Filho/Enteado/Genro/Nora/Neto/Bisneto, Pai/Mãe/Padrasto/Madrasta/Sogro/Sogra/Avô/Avó)
- Região (SE, NE, N, CO, S)
- Posição na ocupação (Com carteira, Sem carteira, Conta-própria, Outros)
Categorias de referência:
- Quem está na força de trabalho (na_forca_trabalho == “0”)
- Quem não recebeu bolsa família (recebeu_bf == “0”)
- Homens
- Até 24 anos
- Brancos
- Responsável/cônjuge no domicílio
- Região sudeste
- Urbano
- Setor privado, setor público ou doméstico com carteira
Filtros:
- Foram mantidos apenas ocupados, pois o objeto de análise é a jornada
- Foram excluídos registros considerados “outliers”, horas de trabalho semanais inferiores a 20 ou superiores a 100 horas
# Function to prepare data for intensive margin analysis
prepare_intensive_margin_data <- function(data) {
data |>
filter(VD4002 %in% "1") |>
mutate(na_forca_trabalho = case_when(VD4001 == "1" ~ "0", VD4001 == "2" ~
"1", TRUE ~ NA_character_), horas_ef = as.numeric(VD4035), log_horas_ef = log(horas_ef +
0.01), rem_ef_todos = as.numeric(VD4020), log_rem_ef_todos = log(rem_ef_todos +
0.01), recebeu_bf = case_when(V5002A == "2" ~ "0", V5002A == "1" ~ "1",
TRUE ~ NA_character_), sexo = as.factor(V2007), idade = as.numeric(V2009),
log_idade = log(idade + 0.01), raca = case_when(V2010 == "1" ~ "B", V2010 %in%
c("2", "4") ~ "PP", V2010 %in% c("3", "5", "9") ~ "O", TRUE ~ NA_character_) |>
as.factor(), anos_estudo = case_when(VD3006 %in% c("1", "2", "3",
"4") ~ "Até 11", VD3006 %in% c("5") ~ "12 a 15", VD3006 %in% c("6") ~
"16 ou mais", TRUE ~ NA_character_), cond_dom = case_when(VD2002 %in%
c("01", "02") ~ "Responsável/Cônjuge", VD2002 %in% c("03", "04",
"05", "08", "09") ~ "Filho/Enteado/Genro/Nora/Neto/Bisneto", VD2002 %in%
c("06", "07", "11") ~ "Pai/Mãe/Padrasto/Madrasta/Sogro/Sogra/Avô/Avó",
TRUE ~ NA_character_), n_pessoas_dom = as.numeric(VD2003), log_n_pessoas_dom = log(n_pessoas_dom +
0.01), renda_dom_pc = as.numeric(VD5005), log_renda_dom_pc = log(renda_dom_pc +
0.01), UF = as.factor(UF), regiao = case_when(UF %in% c("11", "12",
"13", "14", "15", "16", "17") ~ "N", UF %in% c("21", "22", "23",
"24", "25", "26", "27", "28", "29") ~ "NE", UF %in% c("31", "32",
"33", "35") ~ "SE", UF %in% c("41", "42", "43") ~ "S", UF %in% c("50",
"51", "52", "53") ~ "CO", TRUE ~ NA_character_), area_urbana = as.factor(V1022),
capital_regiao = as.factor(V1023), pos_ocup = case_when(VD4009 %in% c("01",
"03", "05") ~ "Com carteira", VD4009 %in% c("02", "04", "06") ~ "Sem carteira",
VD4009 %in% c("09") ~ "Conta própria", VD4009 %in% c("07", "09",
"10") ~ "Outros"), grup_ativ = as.factor(VD4010), grup_ocup = as.factor(VD4011)) |>
mutate(na_forca_trabalho = factor(na_forca_trabalho, levels = c("0", "1")),
recebeu_bf = factor(recebeu_bf, levels = c("0", "1")), anos_estudo = factor(anos_estudo,
levels = c("Até 11", "12 a 15", "16 ou mais")), raca = factor(raca,
levels = c("B", "PP", "O")), cond_dom = factor(cond_dom, levels = c("Responsável/Cônjuge",
"Filho/Enteado/Genro/Nora/Neto/Bisneto", "Pai/Mãe/Padrasto/Madrasta/Sogro/Sogra/Avô/Avó")),
regiao = factor(regiao, levels = c("SE", "NE", "N", "CO", "S"))) |>
mutate(grupo_idade = cut(idade, breaks = c(-Inf, 24, 64, Inf), labels = c("Até 24",
"25 a a 64", "65 ou mais"), include.lowest = TRUE, right = FALSE), grupo_renda_dom_pc = cut(renda_dom_pc,
breaks = c(-Inf, 218, 1412, 7060, Inf), labels = c("Até R$ 218,00",
"R$ 218,01 a 1,00 SM", "1,01 SM a 5,00 SM", "Acima de 5,00 SM"),
include.lowest = TRUE, right = FALSE))
}
dados_pnadc_filtered_ocupados <- prepare_intensive_margin_data(dados_pnadc)
vars_modelo_ocupados <- c("horas_ef", "log_horas_ef", "recebeu_bf", "sexo", "idade",
"log_idade", "grupo_idade", "raca", "anos_estudo", "n_pessoas_dom", "log_n_pessoas_dom",
"renda_dom_pc", "grupo_renda_dom_pc", "log_renda_dom_pc", "cond_dom", "regiao",
"area_urbana", "capital_regiao", "pos_ocup", "grup_ativ", "grup_ocup", "rem_ef_todos",
"log_rem_ef_todos")
dados_pnadc_filtered_ocupados_semp <- dados_pnadc_filtered_ocupados |>
drop_na(all_of(vars_modelo_ocupados)) |>
dplyr::select(all_of(vars_modelo_ocupados))
cat("Amostra para análise da margem intensiva:", nrow(dados_pnadc_filtered_ocupados_semp),
"observações\n")## Amostra para análise da margem intensiva: 154034 observaçõesEspecificação do Modelo
# Function to fit the negative binomial model
fit_nb_model <- function(data) {
glm.nb(horas_ef ~ recebeu_bf + sexo + grupo_idade + anos_estudo + raca + cond_dom +
grupo_renda_dom_pc + regiao + area_urbana + pos_ocup, data = data)
}
glm_base_ocupados_final <- fit_nb_model(dados_pnadc_filtered_ocupados_semp)
vcov_rob_ocupados_3 <- sandwich::vcovHC(glm_base_ocupados_final, type = "HC1")
# Sumário do modelo
summ(glm_base_ocupados_final, exp = TRUE, confint = TRUE, vif = TRUE, robust = "HC1")| Observations | 154034 |
| Dependent variable | horas_ef |
| Type | Generalized linear model |
| Family | Negative Binomial(5.6552) |
| Link | log |
| χ²(154012) | 9922.71 |
| p | 0.00 |
| Pseudo-R² (Cragg-Uhler) | 0.06 |
| Pseudo-R² (McFadden) | 0.01 |
| AIC | 1308890.08 |
| BIC | 1309118.82 |
| exp(Est.) | 2.5% | 97.5% | z val. | p | VIF | |
|---|---|---|---|---|---|---|
| (Intercept) | 33.76 | 32.95 | 34.59 | 283.83 | 0.00 | NA |
| recebeu_bf1 | 0.88 | 0.87 | 0.88 | -24.88 | 0.00 | 1.19 |
| sexo2 | 0.89 | 0.89 | 0.90 | -55.65 | 0.00 | 1.15 |
| grupo_idade25 a a 64 | 1.05 | 1.04 | 1.06 | 14.69 | 0.00 | 1.41 |
| grupo_idade65 ou mais | 0.90 | 0.89 | 0.91 | -15.94 | 0.00 | 1.41 |
| anos_estudo12 a 15 | 1.04 | 1.03 | 1.04 | 15.28 | 0.00 | 1.59 |
| anos_estudo16 ou mais | 1.00 | 1.00 | 1.01 | 0.47 | 0.64 | 1.59 |
| racaPP | 0.99 | 0.99 | 1.00 | -2.82 | 0.00 | 1.22 |
| racaO | 1.00 | 0.98 | 1.02 | -0.39 | 0.70 | 1.22 |
| cond_domFilho/Enteado/Genro/Nora/Neto/Bisneto | 0.97 | 0.96 | 0.97 | -12.87 | 0.00 | 1.33 |
| cond_domPai/Mãe/Padrasto/Madrasta/Sogro/Sogra/Avô/Avó | 1.05 | 1.03 | 1.06 | 6.83 | 0.00 | 1.33 |
| grupo_renda_dom_pcR$ 218,01 a 1,00 SM | 1.20 | 1.17 | 1.22 | 15.30 | 0.00 | 1.49 |
| grupo_renda_dom_pc1,01 SM a 5,00 SM | 1.26 | 1.23 | 1.29 | 19.78 | 0.00 | 1.49 |
| grupo_renda_dom_pcAcima de 5,00 SM | 1.30 | 1.26 | 1.33 | 20.33 | 0.00 | 1.49 |
| regiaoNE | 0.96 | 0.95 | 0.97 | -15.19 | 0.00 | 1.33 |
| regiaoN | 0.99 | 0.98 | 0.99 | -4.34 | 0.00 | 1.33 |
| regiaoCO | 1.00 | 0.99 | 1.01 | -0.24 | 0.81 | 1.33 |
| regiaoS | 1.00 | 1.00 | 1.01 | 0.94 | 0.35 | 1.33 |
| area_urbana2 | 0.99 | 0.99 | 1.00 | -2.39 | 0.02 | 1.15 |
| pos_ocupConta própria | 0.91 | 0.90 | 0.91 | -41.50 | 0.00 | 1.41 |
| pos_ocupOutros | 0.90 | 0.90 | 0.91 | -30.04 | 0.00 | 1.41 |
| pos_ocupSem carteira | 0.88 | 0.87 | 0.88 | -51.70 | 0.00 | 1.41 |
| Standard errors: Robust, type = HC1 |
Diagnósticos e Validação
# Função para diagnósticos completos do modelo binomial negativa
perform_complete_diagnostics_nb <- function(model, data, vcov_model, model_name = "Modelo Binomial Negativa") {
cat("===============================================\n")
cat(paste("DIAGNÓSTICOS E VALIDAÇÃO DO", toupper(model_name), "\n"))
cat("===============================================\n\n")
# 1. ANÁLISE DO AJUSTE DO MODELO
cat("1. ANÁLISE DO AJUSTE DO MODELO\n")
cat("------------------------------\n")
# Sumário com medidas do jtools
cat("Medidas de Ajuste (equivalentes ao R²):\n")
tryCatch({
# Para modelos binomial negativa, calcular pseudo-R² e outras medidas
null_deviance <- model$null.deviance
residual_deviance <- model$deviance
n <- length(model$y)
# McFadden Pseudo-R² (análogo ao R² para GLM)
mcfadden_r2 <- 1 - (residual_deviance/null_deviance)
cat("McFadden Pseudo-R²:", round(mcfadden_r2, 4), "\n")
# Cox & Snell Pseudo-R²
cox_snell_r2 <- 1 - exp((residual_deviance - null_deviance)/n)
cat("Cox & Snell Pseudo-R²:", round(cox_snell_r2, 4), "\n")
# Nagelkerke Pseudo-R² (Cragg-Uhler)
nagelkerke_r2 <- cox_snell_r2/(1 - exp(-null_deviance/n))
cat("Nagelkerke Pseudo-R² (Cragg-Uhler):", round(nagelkerke_r2, 4), "\n")
# Efron's Pseudo-R² (específico para modelos de contagem)
y_obs <- model$y
y_pred <- fitted(model)
y_mean <- mean(y_obs)
efron_r2 <- 1 - sum((y_obs - y_pred)^2)/sum((y_obs - y_mean)^2)
cat("Efron's Pseudo-R²:", round(efron_r2, 4), "\n")
# Medidas adicionais de ajuste
cat("\nOutras medidas de ajuste:\n")
cat("AIC:", round(AIC(model), 2), "\n")
cat("BIC:", round(BIC(model), 2), "\n")
cat("Log-Likelihood:", round(as.numeric(logLik(model)), 2), "\n")
# Proporção de deviance explicada
prop_dev_explained <- (null_deviance - residual_deviance)/null_deviance
cat("Proporção de Deviance Explicada:", round(prop_dev_explained, 4), "\n")
cat("\n")
}, error = function(e) {
cat("Erro ao calcular medidas de ajuste:", e$message, "\n")
cat("Usando informações básicas do modelo:\n")
cat("AIC:", round(AIC(model), 2), "\n")
cat("BIC:", round(BIC(model), 2), "\n")
cat("Log-Likelihood:", round(as.numeric(logLik(model)), 2), "\n")
cat("Deviance:", round(deviance(model), 2), "\n\n")
})
# Outros indicadores de ajuste e performance do modelo
cat("Outros indicadores de ajuste e performance do modelo:\n")
model_perf <- model_performance(model)
print(model_perf)
cat("\n")
# 2. HETEROCEDASTICIDADE
cat("2. HETEROCEDASTICIDADE\n")
cat("----------------------\n")
# Breusch-Pagan
cat("Teste de Breusch-Pagan:\n")
bp_test <- bptest(model)
print(bp_test)
cat("Interpretação: p-valor < 0.05 indica presença de heterocedasticidade\n\n")
# Goldfeld-Quandt
cat("Teste de Goldfeld-Quandt:\n")
tryCatch({
gq_test <- gqtest(model, order.by = fitted(model))
print(gq_test)
cat("Interpretação: p-valor < 0.05 indica presença de heterocedasticidade\n\n")
}, error = function(e) {
cat("Teste GQ não executado para modelo binomial negativa\n")
cat("Motivo: Limitações computacionais ou características do modelo\n\n")
})
# White Test (usando pacote skedastic)
cat("Teste de White:\n")
tryCatch({
# Teste de White
white_test <- skedastic::white(model)
print(white_test)
cat("Interpretação: p-valor < 0.05 indica presença de heterocedasticidade\n\n")
}, error = function(e) {
tryCatch({
# Fallback: tentar whitestrap
white_test <- whitestrap::white_test(model)
print(white_test)
cat("Interpretação: p-valor < 0.05 indica presença de heterocedasticidade\n\n")
})
})
# 3. MULTICOLINEARIDADE
cat("3. MULTICOLINEARIDADE\n")
cat("---------------------\n")
# VIF usando car
cat("Fatores de Inflação da Variância (VIF):\n")
vif_values <- car::vif(model)
print(vif_values)
cat("\nInterpretação:\n")
cat("- VIF < 5: Não há problema de multicolinearidade\n")
cat("- VIF entre 5-10: Multicolinearidade moderada\n")
cat("- VIF > 10: Multicolinearidade severa\n\n")
# Identificar variáveis problemáticas
if (any(vif_values > 10)) {
cat("ATENÇÃO: Variáveis com VIF > 10 (multicolinearidade severa):\n")
high_vif <- vif_values[vif_values > 10]
print(high_vif)
cat("\n")
}
# 4. SIGNIFICÂNCIA INDIVIDUAL
cat("4. SIGNIFICÂNCIA INDIVIDUAL\n")
cat("---------------------------\n")
# Teste de Wald Tipo II
cat("Teste de Wald (Tipo II):\n")
wald_test_2 <- car::Anova(model, type = "II", test.statistic = "Wald", white.adjust = "hc1")
print(wald_test_2)
cat("Interpretação: p-valor < 0.05 indica que a variável é estatisticamente significativa\n\n")
# Teste LRT Tipo II
cat("Teste de Razão de Verossimilhança (Tipo II):\n")
lr_test_2 <- car::Anova(model, type = "II", test.statistic = "LR", white.adjust = "hc1")
print(lr_test_2)
cat("Interpretação: p-valor < 0.05 indica que a variável é estatisticamente significativa\n\n")
# Teste de Wald Tipo III
cat("Teste de Wald (Tipo III):\n")
wald_test_3 <- car::Anova(model, type = "III", test.statistic = "Wald", white.adjust = "hc1")
print(wald_test_3)
cat("Interpretação: p-valor < 0.05 indica que a variável é estatisticamente significativa\n\n")
# Teste LRT Tipo III
cat("Teste de Razão de Verossimilhança (Tipo III):\n")
lr_test_3 <- car::Anova(model, type = "III", test.statistic = "LR", white.adjust = "hc1")
print(lr_test_3)
cat("Interpretação: p-valor < 0.05 indica que a variável é estatisticamente significativa\n\n")
# 5. SIGNIFICÂNCIA CONJUNTA
cat("5.4 Teste de Wald global\n")
cat("---------------------------------------\n")
tryCatch({
# Teste de Wald para todos os coeficientes exceto intercepto
n_coef <- length(coef(model))
if (n_coef > 1) {
terms_to_test <- 2:n_coef
wald_aod <- wald.test(b = coef(model), Sigma = vcov_model, Terms = terms_to_test)
print(wald_aod)
# Interpretação
p_val <- wald_aod$result$chi2[3]
cat("\nInterpretação:")
if (p_val < 0.001) {
cat(" Modelo altamente significativo (p < 0.001) ***\n")
} else if (p_val < 0.01) {
cat(" Modelo muito significativo (p < 0.01) **\n")
} else if (p_val < 0.05) {
cat(" Modelo significativo (p < 0.05) *\n")
} else {
cat(" Modelo não significativo (p >= 0.05)\n")
}
} else {
cat("Modelo tem apenas intercepto - teste não aplicável\n")
}
cat("\n")
}, error = function(e) {
cat("Erro no teste AOD/Wald:", e$message, "\n\n")
})
# 6. ESPECIFICAÇÃO DO MODELO
cat("6. ESPECIFICAÇÃO DO MODELO\n")
cat("--------------------------\n")
# RESET Test
cat("Teste RESET (Ramsey):\n")
tryCatch({
reset_test <- resettest(model, power = 2:3, type = "fitted")
print(reset_test)
cat("Interpretação: p-valor < 0.05 indica problemas de especificação funcional\n\n")
}, error = function(e) {
cat("Teste RESET não executado para modelo binomial negativa\n")
cat("Motivo: Limitações do teste para modelos GLM não-lineares\n\n")
})
# 7. LINEARIDADE
cat("7. LINEARIDADE (PARA MODELOS DE CONTAGEM)\n")
cat("-----------------------------------------\n")
cat("Para modelos binomial negativa, a linearidade se refere à relação\n")
cat("entre as variáveis explicativas e o log da média condicional.\n")
cat("Verificação através dos resíduos de Pearson:\n\n")
# Resíduos de Pearson
pearson_resid <- residuals(model, type = "pearson")
cat("Estatísticas dos Resíduos de Pearson:\n")
cat("Média:", round(mean(pearson_resid), 4), "\n")
cat("Desvio Padrão:", round(sd(pearson_resid), 4), "\n")
cat("Mín:", round(min(pearson_resid), 4), "Máx:", round(max(pearson_resid),
4), "\n\n")
# 8. NORMALIDADE DOS RESÍDUOS
cat("8. NORMALIDADE DOS RESÍDUOS\n")
cat("---------------------------\n")
cat("NOTA: Para modelos binomial negativa, os resíduos não seguem distribuição normal.\n")
cat("Os testes abaixo são informativos, mas não são requisitos do modelo.\n\n")
# Teste de Shapiro-Wilk (apenas para amostras pequenas)
if (length(pearson_resid) <= 5000) {
cat("Teste de Shapiro-Wilk (Resíduos de Pearson):\n")
shapiro_test <- shapiro.test(pearson_resid)
print(shapiro_test)
cat("\n")
} else {
cat("Teste de Shapiro-Wilk não executado (amostra muito grande)\n\n")
}
# Teste de Jarque-Bera
cat("Teste de Jarque-Bera (Resíduos de Pearson):\n")
jb_test <- jarque.bera.test(pearson_resid)
print(jb_test)
cat("Interpretação: Para modelos de contagem, rejeição da normalidade é esperada\n\n")
# 9. SOBREDISPERSÃO
cat("9. SOBREDISPERSÃO\n")
cat("-----------------\n")
cat("Teste específico para modelos de contagem:\n")
# Teste de sobredispersão
overdispersion_test <- check_overdispersion(model)
print(overdispersion_test)
cat("\nInterpretação:\n")
cat("- p-valor < 0.05: Sobredispersão estatisticamente significativa\n\n")
# 10. ENDOGENEIDADE (COM DEBUG)
cat("10. ENDOGENEIDADE\n")
cat("-----------------\n")
cat("Análise de correlação entre resíduos e variáveis explicativas\n")
cat("(Teste indireto para detecção de endogeneidade)\n\n")
model_resid <- residuals(model, type = "pearson")
test_vars <- c("recebeu_bf", "sexo", "grupo_idade", "raca", "anos_estudo", "cond_dom",
"grupo_renda_dom_pc", "regiao", "area_urbana")
test_vars <- test_vars[test_vars %in% names(data)]
correlation_results <- data.frame(Variavel = character(), Correlacao = numeric(),
P_valor = numeric(), Significativo = character(), stringsAsFactors = FALSE)
cat("Correlações entre resíduos e variáveis explicativas:\n")
cat("---------------------------------------------------\n")
cat("Resíduos Pearson vs Correlação Pearson:\n")
cat("---------------------------------------------------\n")
for (i in 1:length(test_vars)) {
var <- test_vars[i]
tryCatch({
# Converter variável para numérica se necessário
if (is.factor(data[[var]])) {
var_numeric <- as.numeric(data[[var]])
} else if (is.character(data[[var]])) {
var_numeric <- as.numeric(as.factor(data[[var]]))
} else {
var_numeric <- as.numeric(data[[var]])
}
valid_indices <- !is.na(var_numeric) & !is.na(model_resid)
# Testar se há observações suficientes
if (sum(valid_indices) > 10) {
# Teste de correlação (usando Spearman por robustez)
cor_test <- cor.test(model_resid[valid_indices], var_numeric[valid_indices],
method = "pearson")
# Determinar significância
significativo <- ifelse(cor_test$p.value < 0.05, "Sim", "Não")
# Adicionar aos resultados
correlation_results <- rbind(correlation_results, data.frame(Variavel = var,
Correlacao = round(cor_test$estimate, 4), P_valor = round(cor_test$p.value,
6), Significativo = significativo))
# Mostrar resultado
if (var == "recebeu_bf") {
cat("*** ", var, " (VARIÁVEL PRINCIPAL): ", round(cor_test$estimate,
4), " (p = ", round(cor_test$p.value, 6), ") ", significativo,
" ***\n", sep = "")
} else {
cat(var, ": ", round(cor_test$estimate, 4), " (p = ", round(cor_test$p.value,
4), ") ", significativo, "\n", sep = "")
}
} else {
cat(var, ": Dados insuficientes para teste\n")
}
}, error = function(e) {
cat("ERRO na variável", var, ":", e$message, "\n")
})
}
cat("---------------------------------------------------\n")
cat("Resíduos Pearson vs Correlação Spearman:\n")
cat("---------------------------------------------------\n")
for (i in 1:length(test_vars)) {
var <- test_vars[i]
tryCatch({
# Converter variável para numérica se necessário
if (is.factor(data[[var]])) {
var_numeric <- as.numeric(data[[var]])
} else if (is.character(data[[var]])) {
var_numeric <- as.numeric(as.factor(data[[var]]))
} else {
var_numeric <- as.numeric(data[[var]])
}
valid_indices <- !is.na(var_numeric) & !is.na(model_resid)
# Testar se há observações suficientes
if (sum(valid_indices) > 10) {
# Teste de correlação (usando Spearman por robustez)
cor_test <- cor.test(model_resid[valid_indices], var_numeric[valid_indices],
method = "spearman")
# Determinar significância
significativo <- ifelse(cor_test$p.value < 0.05, "Sim", "Não")
# Adicionar aos resultados
correlation_results <- rbind(correlation_results, data.frame(Variavel = var,
Correlacao = round(cor_test$estimate, 4), P_valor = round(cor_test$p.value,
6), Significativo = significativo))
# Mostrar resultado
if (var == "recebeu_bf") {
cat("*** ", var, " (VARIÁVEL PRINCIPAL): ", round(cor_test$estimate,
4), " (p = ", round(cor_test$p.value, 6), ") ", significativo,
" ***\n", sep = "")
} else {
cat(var, ": ", round(cor_test$estimate, 4), " (p = ", round(cor_test$p.value,
4), ") ", significativo, "\n", sep = "")
}
} else {
cat(var, ": Dados insuficientes para teste\n")
}
}, error = function(e) {
cat("ERRO na variável", var, ":", e$message, "\n")
})
}
model_resid <- residuals(model, type = "deviance")
cat("---------------------------------------------------\n")
cat("Resíduos Deviance vs Correlação Pearson:\n")
cat("---------------------------------------------------\n")
for (i in 1:length(test_vars)) {
var <- test_vars[i]
tryCatch({
# Converter variável para numérica se necessário
if (is.factor(data[[var]])) {
var_numeric <- as.numeric(data[[var]])
} else if (is.character(data[[var]])) {
var_numeric <- as.numeric(as.factor(data[[var]]))
} else {
var_numeric <- as.numeric(data[[var]])
}
valid_indices <- !is.na(var_numeric) & !is.na(model_resid)
# Testar se há observações suficientes
if (sum(valid_indices) > 10) {
# Teste de correlação (usando Spearman por robustez)
cor_test <- cor.test(model_resid[valid_indices], var_numeric[valid_indices],
method = "pearson")
# Determinar significância
significativo <- ifelse(cor_test$p.value < 0.05, "Sim", "Não")
# Adicionar aos resultados
correlation_results <- rbind(correlation_results, data.frame(Variavel = var,
Correlacao = round(cor_test$estimate, 4), P_valor = round(cor_test$p.value,
6), Significativo = significativo))
# Mostrar resultado
if (var == "recebeu_bf") {
cat("*** ", var, " (VARIÁVEL PRINCIPAL): ", round(cor_test$estimate,
4), " (p = ", round(cor_test$p.value, 6), ") ", significativo,
" ***\n", sep = "")
} else {
cat(var, ": ", round(cor_test$estimate, 4), " (p = ", round(cor_test$p.value,
4), ") ", significativo, "\n", sep = "")
}
} else {
cat(var, ": Dados insuficientes para teste\n")
}
}, error = function(e) {
cat("ERRO na variável", var, ":", e$message, "\n")
})
}
cat("---------------------------------------------------\n")
cat("Resíduos Deviance vs Correlação Spearman:\n")
cat("---------------------------------------------------\n")
for (i in 1:length(test_vars)) {
var <- test_vars[i]
tryCatch({
# Converter variável para numérica se necessário
if (is.factor(data[[var]])) {
var_numeric <- as.numeric(data[[var]])
} else if (is.character(data[[var]])) {
var_numeric <- as.numeric(as.factor(data[[var]]))
} else {
var_numeric <- as.numeric(data[[var]])
}
valid_indices <- !is.na(var_numeric) & !is.na(model_resid)
# Testar se há observações suficientes
if (sum(valid_indices) > 10) {
# Teste de correlação (usando Spearman por robustez)
cor_test <- cor.test(model_resid[valid_indices], var_numeric[valid_indices],
method = "spearman")
# Determinar significância
significativo <- ifelse(cor_test$p.value < 0.05, "Sim", "Não")
# Adicionar aos resultados
correlation_results <- rbind(correlation_results, data.frame(Variavel = var,
Correlacao = round(cor_test$estimate, 4), P_valor = round(cor_test$p.value,
6), Significativo = significativo))
# Mostrar resultado
if (var == "recebeu_bf") {
cat("*** ", var, " (VARIÁVEL PRINCIPAL): ", round(cor_test$estimate,
4), " (p = ", round(cor_test$p.value, 6), ") ", significativo,
" ***\n", sep = "")
} else {
cat(var, ": ", round(cor_test$estimate, 4), " (p = ", round(cor_test$p.value,
4), ") ", significativo, "\n", sep = "")
}
} else {
cat(var, ": Dados insuficientes para teste\n")
}
}, error = function(e) {
cat("ERRO na variável", var, ":", e$message, "\n")
})
}
cat("===============================================\n")
cat("FIM DOS DIAGNÓSTICOS\n")
cat("===============================================\n\n")
}
# Executar diagnósticos para o modelo da margem intensiva
perform_complete_diagnostics_nb(glm_base_ocupados_final, dados_pnadc_filtered_ocupados_semp,
vcov_rob_ocupados_3, "Margem Intensiva")## ===============================================
## DIAGNÓSTICOS E VALIDAÇÃO DO MARGEM INTENSIVA
## ===============================================
##
## 1. ANÁLISE DO AJUSTE DO MODELO
## ------------------------------
## Medidas de Ajuste (equivalentes ao R²):
## McFadden Pseudo-R²: 0.0517
## Cox & Snell Pseudo-R²: 0.0648
## Nagelkerke Pseudo-R² (Cragg-Uhler): 0.0893
## Efron's Pseudo-R²: 0.0918
##
## Outras medidas de ajuste:
## AIC: 1308890
## BIC: 1309119
## Log-Likelihood: -654422
## Proporção de Deviance Explicada: 0.0517
##
## Outros indicadores de ajuste e performance do modelo:
## # Indices of model performance
##
## AIC | AICc | BIC | Nagelkerke's R2 | RMSE | Sigma
## --------------------------------------------------------------------
## 1.309e+06 | 1.309e+06 | 1.309e+06 | 0.089 | 13.254 | 1.000
##
## AIC | Score_log | Score_spherical
## ---------------------------------------
## 1.309e+06 | -4.409 | 0.002
##
## 2. HETEROCEDASTICIDADE
## ----------------------
## Teste de Breusch-Pagan:
##
## studentized Breusch-Pagan test
##
## data: model
## BP = 3493.9, df = 21, p-value < 2.2e-16
##
## Interpretação: p-valor < 0.05 indica presença de heterocedasticidade
##
## Teste de Goldfeld-Quandt:
##
## Goldfeld-Quandt test
##
## data: model
## GQ = 0.69529, df1 = 76995, df2 = 76995, p-value = 1
## alternative hypothesis: variance increases from segment 1 to 2
##
## Interpretação: p-valor < 0.05 indica presença de heterocedasticidade
##
## Teste de White:
## White's test results
##
## Null hypothesis: Homoskedasticity of the residuals
## Alternative hypothesis: Heteroskedasticity of the residuals
## Test Statistic: 7259.21
## P-value: 0
## Interpretação: p-valor < 0.05 indica presença de heterocedasticidade
##
## 3. MULTICOLINEARIDADE
## ---------------------
## Fatores de Inflação da Variância (VIF):
## GVIF Df GVIF^(1/(2*Df))
## recebeu_bf 1.187351 1 1.089656
## sexo 1.153626 1 1.074070
## grupo_idade 1.407938 2 1.089296
## anos_estudo 1.590404 2 1.122993
## raca 1.223383 2 1.051697
## cond_dom 1.328681 2 1.073631
## grupo_renda_dom_pc 1.489755 3 1.068692
## regiao 1.332165 4 1.036501
## area_urbana 1.151780 1 1.073210
## pos_ocup 1.411449 3 1.059118
##
## Interpretação:
## - VIF < 5: Não há problema de multicolinearidade
## - VIF entre 5-10: Multicolinearidade moderada
## - VIF > 10: Multicolinearidade severa
##
## 4. SIGNIFICÂNCIA INDIVIDUAL
## ---------------------------
## Teste de Wald (Tipo II):
## Analysis of Deviance Table (Type II tests)
##
## Response: horas_ef
## Df Chisq Pr(>Chisq)
## recebeu_bf 1 716.4696 < 2e-16 ***
## sexo 1 2084.6942 < 2e-16 ***
## grupo_idade 2 890.2857 < 2e-16 ***
## anos_estudo 2 198.7110 < 2e-16 ***
## raca 2 5.0812 0.07882 .
## cond_dom 2 120.8232 < 2e-16 ***
## grupo_renda_dom_pc 3 827.8975 < 2e-16 ***
## regiao 4 208.8503 < 2e-16 ***
## area_urbana 1 4.2762 0.03865 *
## pos_ocup 3 2059.5951 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Interpretação: p-valor < 0.05 indica que a variável é estatisticamente significativa
##
## Teste de Razão de Verossimilhança (Tipo II):
## Analysis of Deviance Table (Type II tests)
##
## Response: horas_ef
## LR Chisq Df Pr(>Chisq)
## recebeu_bf 706.46 1 < 2e-16 ***
## sexo 2067.26 1 < 2e-16 ***
## grupo_idade 872.98 2 < 2e-16 ***
## anos_estudo 198.70 2 < 2e-16 ***
## raca 5.09 2 0.07845 .
## cond_dom 121.80 2 < 2e-16 ***
## grupo_renda_dom_pc 819.29 3 < 2e-16 ***
## regiao 209.60 4 < 2e-16 ***
## area_urbana 4.28 1 0.03845 *
## pos_ocup 2088.15 3 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Interpretação: p-valor < 0.05 indica que a variável é estatisticamente significativa
##
## Teste de Wald (Tipo III):
## Analysis of Deviance Table (Type III tests)
##
## Response: horas_ef
## Df Chisq Pr(>Chisq)
## (Intercept) 1 96320.8848 < 2e-16 ***
## recebeu_bf 1 716.4696 < 2e-16 ***
## sexo 1 2084.6942 < 2e-16 ***
## grupo_idade 2 890.2857 < 2e-16 ***
## anos_estudo 2 198.7110 < 2e-16 ***
## raca 2 5.0812 0.07882 .
## cond_dom 2 120.8232 < 2e-16 ***
## grupo_renda_dom_pc 3 827.8975 < 2e-16 ***
## regiao 4 208.8503 < 2e-16 ***
## area_urbana 1 4.2762 0.03865 *
## pos_ocup 3 2059.5951 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Interpretação: p-valor < 0.05 indica que a variável é estatisticamente significativa
##
## Teste de Razão de Verossimilhança (Tipo III):
## Analysis of Deviance Table (Type III tests)
##
## Response: horas_ef
## LR Chisq Df Pr(>Chisq)
## recebeu_bf 706.46 1 < 2e-16 ***
## sexo 2067.26 1 < 2e-16 ***
## grupo_idade 872.98 2 < 2e-16 ***
## anos_estudo 198.70 2 < 2e-16 ***
## raca 5.09 2 0.07845 .
## cond_dom 121.80 2 < 2e-16 ***
## grupo_renda_dom_pc 819.29 3 < 2e-16 ***
## regiao 209.60 4 < 2e-16 ***
## area_urbana 4.28 1 0.03845 *
## pos_ocup 2088.15 3 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Interpretação: p-valor < 0.05 indica que a variável é estatisticamente significativa
##
## 5.4 Teste de Wald global
## ---------------------------------------
## Wald test:
## ----------
##
## Chi-squared test:
## X2 = 14006.2, df = 21, P(> X2) = 0.0
##
## Interpretação: Modelo altamente significativo (p < 0.001) ***
##
## 6. ESPECIFICAÇÃO DO MODELO
## --------------------------
## Teste RESET (Ramsey):
##
## RESET test
##
## data: model
## RESET = 437.76, df1 = 2, df2 = 154010, p-value < 2.2e-16
##
## Interpretação: p-valor < 0.05 indica problemas de especificação funcional
##
## 7. LINEARIDADE (PARA MODELOS DE CONTAGEM)
## -----------------------------------------
## Para modelos binomial negativa, a linearidade se refere à relação
## entre as variáveis explicativas e o log da média condicional.
## Verificação através dos resíduos de Pearson:
##
## Estatísticas dos Resíduos de Pearson:
## Média: 0
## Desvio Padrão: 0.8102
## Mín: -2.25 Máx: 6.5601
##
## 8. NORMALIDADE DOS RESÍDUOS
## ---------------------------
## NOTA: Para modelos binomial negativa, os resíduos não seguem distribuição normal.
## Os testes abaixo são informativos, mas não são requisitos do modelo.
##
## Teste de Shapiro-Wilk não executado (amostra muito grande)
##
## Teste de Jarque-Bera (Resíduos de Pearson):
##
## Jarque Bera Test
##
## data: pearson_resid
## X-squared = 44722, df = 2, p-value < 2.2e-16
##
## Interpretação: Para modelos de contagem, rejeição da normalidade é esperada
##
## 9. SOBREDISPERSÃO
## -----------------
## Teste específico para modelos de contagem:
## # Overdispersion test
##
## dispersion ratio = 0.597
## p-value = < 0.001##
## Interpretação:
## - p-valor < 0.05: Sobredispersão estatisticamente significativa
##
## 10. ENDOGENEIDADE
## -----------------
## Análise de correlação entre resíduos e variáveis explicativas
## (Teste indireto para detecção de endogeneidade)
##
## Correlações entre resíduos e variáveis explicativas:
## ---------------------------------------------------
## Resíduos Pearson vs Correlação Pearson:
## ---------------------------------------------------
## *** recebeu_bf (VARIÁVEL PRINCIPAL): -2e-04 (p = 0.933132) Não ***
## sexo: -3e-04 (p = 0.9104) Não
## grupo_idade: 0 (p = 0.9924) Não
## raca: -2e-04 (p = 0.9395) Não
## anos_estudo: 3e-04 (p = 0.8959) Não
## cond_dom: 1e-04 (p = 0.9643) Não
## grupo_renda_dom_pc: 4e-04 (p = 0.8803) Não
## regiao: 1e-04 (p = 0.9762) Não
## area_urbana: -3e-04 (p = 0.8928) Não
## ---------------------------------------------------
## Resíduos Pearson vs Correlação Spearman:
## ---------------------------------------------------## *** recebeu_bf (VARIÁVEL PRINCIPAL): 0.0109 (p = 2e-05) Sim ***## sexo: 0.0561 (p = 0) Sim## grupo_idade: -0.0036 (p = 0.1632) Não## raca: 0.0157 (p = 0) Sim## anos_estudo: -0.0218 (p = 0) Sim## cond_dom: 0.0069 (p = 0.007) Sim## grupo_renda_dom_pc: -0.0463 (p = 0) Sim## regiao: -0.0041 (p = 0.1056) Não## area_urbana: 0.0154 (p = 0) Sim
## ---------------------------------------------------
## Resíduos Deviance vs Correlação Pearson:
## ---------------------------------------------------
## *** recebeu_bf (VARIÁVEL PRINCIPAL): -0.0252 (p = 0) Sim ***
## sexo: -0.0275 (p = 0) Sim
## grupo_idade: -0.0139 (p = 0) Sim
## raca: -0.0042 (p = 0.0988) Não
## anos_estudo: 0.0083 (p = 0.0011) Sim
## cond_dom: 0.0149 (p = 0) Sim
## grupo_renda_dom_pc: 0.0091 (p = 4e-04) Sim
## regiao: 0.0041 (p = 0.1059) Não
## area_urbana: -0.0089 (p = 5e-04) Sim
## ---------------------------------------------------
## Resíduos Deviance vs Correlação Spearman:
## ---------------------------------------------------## *** recebeu_bf (VARIÁVEL PRINCIPAL): 0.0109 (p = 2e-05) Sim ***## sexo: 0.0561 (p = 0) Sim## grupo_idade: -0.0036 (p = 0.163) Não## raca: 0.0157 (p = 0) Sim## anos_estudo: -0.0218 (p = 0) Sim## cond_dom: 0.0069 (p = 0.007) Sim## grupo_renda_dom_pc: -0.0463 (p = 0) Sim## regiao: -0.0041 (p = 0.1056) Não## area_urbana: 0.0154 (p = 0) Sim
## ===============================================
## FIM DOS DIAGNÓSTICOS
## ===============================================Efeitos Marginais
Sexo
vcov_rob_ocupados_3 <- sandwich::vcovHC(glm_base_ocupados_final, type = "HC1")
# Análise por diferentes características
cat("# EFEITOS MARGINAIS - MARGEM INTENSIVA\n")## # EFEITOS MARGINAIS - MARGEM INTENSIVAcat("=====================================\n")## =====================================# Por sexo
result_sexo <- analyze_and_plot_marginal_effects(glm_base_ocupados_final, "recebeu_bf",
"sexo", vcov_rob_ocupados_3, "Horas Preditas", "por sexo", comparison_type = "slopes")##
## === EFEITO DO BOLSA FAMÍLIA POR SEXO ===
##
## --- EFEITOS MARGINAIS ---
##
## sexo Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## 1 -4.96 0.187 -26.5 <0.001 512.9 -5.33 -4.60
## 2 -4.48 0.170 -26.3 <0.001 502.8 -4.81 -4.14
##
## Term: recebeu_bf
## Type: response
## Comparison: 1 - 0
##
##
## --- VALORES PREDITOS ---
##
## recebeu_bf sexo Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## 0 1 39.9 0.0441 903 <0.001 Inf 39.8 39.9
## 0 2 36.3 0.0568 639 <0.001 Inf 36.2 36.4
## 1 1 32.5 0.1701 191 <0.001 Inf 32.1 32.8
## 1 2 29.4 0.1474 199 <0.001 Inf 29.1 29.7
##
## Type: response
##
##
## --- VERIFICAÇÃO DE CONSISTÊNCIA ---
## Teste Wald (robusto) de diferenças entre sexo:
##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b1=b2 -0.484 0.0189 -25.7 <0.001 480.9 -0.521 -0.447
## b1=0 -4.963 0.1870 -26.5 <0.001 512.9 -5.329 -4.596
## b2=0 -4.478 0.1705 -26.3 <0.001 502.8 -4.813 -4.144Idade
# Por idade
result_idade <- analyze_and_plot_marginal_effects(glm_base_ocupados_final, "recebeu_bf",
"grupo_idade", vcov_rob_ocupados_3, "Horas Preditas", "por Quartil de Idade",
comparison_type = "slopes")##
## === EFEITO DO BOLSA FAMÍLIA POR GRUPO_IDADE ===
##
## --- EFEITOS MARGINAIS ---
##
## grupo_idade Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## Até 24 -4.53 0.171 -26.4 <0.001 509.6 -4.87 -4.20
## 25 a a 64 -4.82 0.182 -26.4 <0.001 508.5 -5.18 -4.46
## 65 ou mais -4.14 0.158 -26.3 <0.001 502.8 -4.45 -3.84
##
## Term: recebeu_bf
## Type: response
## Comparison: 1 - 0
##
##
## --- VALORES PREDITOS ---
##
## recebeu_bf grupo_idade Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## 0 Até 24 36.5 0.1005 363 <0.001 Inf 36.3 36.7
## 0 25 a a 64 39.0 0.0379 1029 <0.001 Inf 39.0 39.1
## 0 65 ou mais 33.3 0.1790 186 <0.001 Inf 32.9 33.6
## 1 Até 24 28.7 0.1683 170 <0.001 Inf 28.4 29.0
## 1 25 a a 64 30.1 0.1500 200 <0.001 Inf 29.8 30.3
## 1 65 ou mais 26.5 0.1966 135 <0.001 Inf 26.1 26.9
##
## Type: response
##
##
## --- VERIFICAÇÃO DE CONSISTÊNCIA ---
## Teste Wald (robusto) de diferenças entre grupo_idade:
##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b1=b2 0.288 0.0177 16.3 <0.001 195.7 0.253 0.322
## b1=b3 -0.389 0.0295 -13.2 <0.001 129.0 -0.446 -0.331
## b2=b3 -0.676 0.0347 -19.5 <0.001 277.8 -0.744 -0.608
## b1=0 -4.533 0.1714 -26.4 <0.001 509.6 -4.869 -4.198
## b2=0 -4.821 0.1825 -26.4 <0.001 508.5 -5.179 -4.463
## b3=0 -4.145 0.1578 -26.3 <0.001 502.8 -4.454 -3.836Raça
# Por raça
result_raca <- analyze_and_plot_marginal_effects(glm_base_ocupados_final, "recebeu_bf",
"raca", vcov_rob_ocupados_3, "Horas Preditas", "por Raça", comparison_type = "slopes")##
## === EFEITO DO BOLSA FAMÍLIA POR RACA ===
##
## --- EFEITOS MARGINAIS ---
##
## raca Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## B -4.84 0.183 -26.5 <0.001 511.3 -5.20 -4.49
## PP -4.70 0.178 -26.4 <0.001 507.1 -5.04 -4.35
## O -4.71 0.185 -25.5 <0.001 475.7 -5.08 -4.35
##
## Term: recebeu_bf
## Type: response
## Comparison: 1 - 0
##
##
## --- VALORES PREDITOS ---
##
## recebeu_bf raca Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## 0 B 39.0 0.0549 710.1 <0.001 Inf 38.9 39.1
## 0 PP 38.0 0.0460 825.7 <0.001 Inf 37.9 38.1
## 0 O 38.2 0.3694 103.4 <0.001 Inf 37.5 38.9
## 1 B 30.4 0.1589 191.3 <0.001 Inf 30.1 30.7
## 1 PP 29.8 0.1494 199.5 <0.001 Inf 29.5 30.1
## 1 O 30.0 0.3239 92.5 <0.001 Inf 29.3 30.6
##
## Type: response
##
##
## --- VERIFICAÇÃO DE CONSISTÊNCIA ---
## Teste Wald (robusto) de diferenças entre raca:
##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b1=b2 -0.1486 0.0102 -14.564 < 0.001 157.2 -0.1686 -0.1286
## b1=b3 -0.1302 0.0463 -2.815 0.00488 7.7 -0.2209 -0.0396
## b2=b3 0.0183 0.0460 0.399 0.68970 0.5 -0.0717 0.1084
## b1=0 -4.8437 0.1828 -26.492 < 0.001 511.3 -5.2020 -4.4853
## b2=0 -4.6951 0.1780 -26.381 < 0.001 507.1 -5.0439 -4.3463
## b3=0 -4.7134 0.1845 -25.543 < 0.001 475.7 -5.0751 -4.3518Anos de estudo
# Por anos_estudo
result_anos_estudo <- analyze_and_plot_marginal_effects(glm_base_ocupados_final,
"recebeu_bf", "anos_estudo", vcov_rob_ocupados_3, "Horas Preditas", "por Faixa de Anos de Estudos",
comparison_type = "slopes")##
## === EFEITO DO BOLSA FAMÍLIA POR ANOS_ESTUDO ===
##
## --- EFEITOS MARGINAIS ---
##
## anos_estudo Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## Até 11 -4.61 0.175 -26.3 <0.001 505.2 -4.95 -4.27
## 12 a 15 -4.89 0.185 -26.4 <0.001 509.5 -5.25 -4.53
## 16 ou mais -4.75 0.179 -26.5 <0.001 512.5 -5.11 -4.40
##
## Term: recebeu_bf
## Type: response
## Comparison: 1 - 0
##
##
## --- VALORES PREDITOS ---
##
## recebeu_bf anos_estudo Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## 0 Até 11 37.3 0.0607 615 <0.001 Inf 37.2 37.5
## 0 12 a 15 39.5 0.0514 768 <0.001 Inf 39.4 39.6
## 0 16 ou mais 38.2 0.0777 491 <0.001 Inf 38.0 38.3
## 1 Até 11 29.6 0.1502 197 <0.001 Inf 29.3 29.9
## 1 12 a 15 30.5 0.1568 195 <0.001 Inf 30.2 30.8
## 1 16 ou mais 29.5 0.1639 180 <0.001 Inf 29.2 29.8
##
## Type: response
##
##
## --- VERIFICAÇÃO DE CONSISTÊNCIA ---
## Teste Wald (robusto) de diferenças entre anos_estudo:
##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b1=b2 0.278 0.0140 19.9 <0.001 291.0 0.251 0.305
## b1=b3 0.146 0.0129 11.3 <0.001 95.7 0.120 0.171
## b2=b3 -0.132 0.0130 -10.2 <0.001 78.9 -0.158 -0.107
## b1=0 -4.609 0.1750 -26.3 <0.001 505.2 -4.952 -4.266
## b2=0 -4.887 0.1848 -26.4 <0.001 509.5 -5.249 -4.525
## b3=0 -4.755 0.1793 -26.5 <0.001 512.5 -5.106 -4.404Condição no domicílio
# Por cond_dom
result_cond_dom <- analyze_and_plot_marginal_effects(glm_base_ocupados_final, "recebeu_bf",
"cond_dom", vcov_rob_ocupados_3, "Horas Preditas", "por Condição no Domicílio",
comparison_type = "slopes")##
## === EFEITO DO BOLSA FAMÍLIA POR COND_DOM ===
##
## --- EFEITOS MARGINAIS ---
##
## cond_dom Estimate Std. Error z
## Responsável/Cônjuge -4.78 0.181 -26.4
## Filho/Enteado/Genro/Nora/Neto/Bisneto -4.66 0.176 -26.5
## Pai/Mãe/Padrasto/Madrasta/Sogro/Sogra/Avô/Avó -4.87 0.187 -26.1
## Pr(>|z|) S 2.5 % 97.5 %
## <0.001 508.0 -5.13 -4.42
## <0.001 511.9 -5.01 -4.32
## <0.001 496.1 -5.24 -4.50
##
## Term: recebeu_bf
## Type: response
## Comparison: 1 - 0
##
##
## --- VALORES PREDITOS ---
##
## recebeu_bf cond_dom Estimate Std. Error
## 0 Responsável/Cônjuge 38.7 0.0407
## 0 Filho/Enteado/Genro/Nora/Neto/Bisneto 37.5 0.0698
## 0 Pai/Mãe/Padrasto/Madrasta/Sogro/Sogra/Avô/Avó 39.2 0.2585
## 1 Responsável/Cônjuge 29.9 0.1496
## 1 Filho/Enteado/Genro/Nora/Neto/Bisneto 29.6 0.1602
## 1 Pai/Mãe/Padrasto/Madrasta/Sogro/Sogra/Avô/Avó 31.2 0.2579
## z Pr(>|z|) S 2.5 % 97.5 %
## 950 <0.001 Inf 38.6 38.8
## 537 <0.001 Inf 37.3 37.6
## 152 <0.001 Inf 38.7 39.7
## 200 <0.001 Inf 29.6 30.2
## 185 <0.001 Inf 29.3 29.9
## 121 <0.001 Inf 30.7 31.7
##
## Type: response
##
##
## --- VERIFICAÇÃO DE CONSISTÊNCIA ---
## Teste Wald (robusto) de diferenças entre cond_dom:
##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b1=b2 -0.1132 0.0113 -10.00 < 0.001 75.9 -0.1354 -0.091
## b1=b3 0.0942 0.0327 2.88 0.00392 8.0 0.0302 0.158
## b2=b3 0.2074 0.0343 6.05 < 0.001 29.4 0.1402 0.275
## b1=0 -4.7762 0.1809 -26.41 < 0.001 508.0 -5.1307 -4.422
## b2=0 -4.6630 0.1759 -26.51 < 0.001 511.9 -5.0078 -4.318
## b3=0 -4.8704 0.1867 -26.09 < 0.001 496.1 -5.2362 -4.505Renda Dom. P.C.
# Por renda
result_renda <- analyze_and_plot_marginal_effects(glm_base_ocupados_final, "recebeu_bf",
"grupo_renda_dom_pc", vcov_rob_ocupados_3, "Horas Preditas", "por Quartil de Renda per Capita",
comparison_type = "slopes")##
## === EFEITO DO BOLSA FAMÍLIA POR GRUPO_RENDA_DOM_PC ===
##
## --- EFEITOS MARGINAIS ---
##
## grupo_renda_dom_pc Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## Até R$ 218,00 -3.72 0.148 -25.1 <0.001 461.2 -4.01 -3.43
## R$ 218,01 a 1,00 SM -4.56 0.174 -26.3 <0.001 503.0 -4.90 -4.22
## 1,01 SM a 5,00 SM -4.92 0.185 -26.5 <0.001 512.8 -5.28 -4.56
## Acima de 5,00 SM -4.97 0.188 -26.4 <0.001 506.8 -5.34 -4.60
##
## Term: recebeu_bf
## Type: response
## Comparison: 1 - 0
##
##
## --- VALORES PREDITOS ---
##
## recebeu_bf grupo_renda_dom_pc Estimate Std. Error z Pr(>|z|) S 2.5 %
## 0 Até R$ 218,00 30.0 0.3476 86.4 <0.001 Inf 29.4
## 0 R$ 218,01 a 1,00 SM 37.0 0.0556 665.7 <0.001 Inf 36.9
## 0 1,01 SM a 5,00 SM 39.5 0.0465 848.8 <0.001 Inf 39.4
## 0 Acima de 5,00 SM 39.8 0.1869 213.2 <0.001 Inf 39.5
## 1 Até R$ 218,00 24.6 0.3077 79.9 <0.001 Inf 24.0
## 1 R$ 218,01 a 1,00 SM 29.7 0.1487 199.7 <0.001 Inf 29.4
## 1 1,01 SM a 5,00 SM 32.2 0.1694 189.9 <0.001 Inf 31.8
## 1 Acima de 5,00 SM 34.0 0.2516 134.9 <0.001 Inf 33.5
## 97.5 %
## 30.7
## 37.1
## 39.6
## 40.2
## 25.2
## 30.0
## 32.5
## 34.5
##
## Type: response
##
##
## --- VERIFICAÇÃO DE CONSISTÊNCIA ---
## Teste Wald (robusto) de diferenças entre grupo_renda_dom_pc:
##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b1=b2 0.8425 0.0541 15.6 <0.001 179.6 0.73661 0.9485
## b1=b3 1.1979 0.0616 19.4 <0.001 277.0 1.07707 1.3187
## b1=b4 1.2462 0.0667 18.7 <0.001 256.3 1.11547 1.3769
## b2=b3 0.3553 0.0148 24.0 <0.001 421.4 0.32636 0.3843
## b2=b4 0.4037 0.0277 14.6 <0.001 157.5 0.34939 0.4579
## b3=b4 0.0483 0.0241 2.0 0.045 4.5 0.00107 0.0956
## b1=0 -3.7210 0.1480 -25.1 <0.001 461.2 -4.01098 -3.4310
## b2=0 -4.5635 0.1737 -26.3 <0.001 503.0 -4.90396 -4.2231
## b3=0 -4.9189 0.1854 -26.5 <0.001 512.8 -5.28227 -4.5555
## b4=0 -4.9672 0.1883 -26.4 <0.001 506.8 -5.33633 -4.5981Região
# Por região
result_regiao <- analyze_and_plot_marginal_effects(glm_base_ocupados_final, "recebeu_bf",
"regiao", vcov_rob_ocupados_3, "Horas Preditas", "por Região", comparison_type = "slopes")##
## === EFEITO DO BOLSA FAMÍLIA POR REGIAO ===
##
## --- EFEITOS MARGINAIS ---
##
## regiao Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## SE -4.88 0.184 -26.5 <0.001 510.9 -5.24 -4.52
## NE -4.50 0.171 -26.3 <0.001 503.4 -4.83 -4.16
## N -4.66 0.177 -26.3 <0.001 504.4 -5.01 -4.31
## CO -4.87 0.184 -26.4 <0.001 508.7 -5.23 -4.51
## S -4.94 0.186 -26.5 <0.001 511.5 -5.30 -4.57
##
## Term: recebeu_bf
## Type: response
## Comparison: 1 - 0
##
##
## --- VALORES PREDITOS ---
##
## recebeu_bf regiao Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## 0 SE 39.3 0.0661 595 <0.001 Inf 39.1 39.4
## 0 NE 36.5 0.0698 523 <0.001 Inf 36.3 36.6
## 0 N 37.7 0.0932 404 <0.001 Inf 37.5 37.8
## 0 CO 39.3 0.1053 373 <0.001 Inf 39.1 39.5
## 0 S 39.7 0.0781 508 <0.001 Inf 39.6 39.9
## 1 SE 30.9 0.1647 188 <0.001 Inf 30.6 31.2
## 1 NE 29.2 0.1499 195 <0.001 Inf 28.9 29.5
## 1 N 30.2 0.1642 184 <0.001 Inf 29.9 30.6
## 1 CO 30.8 0.1759 175 <0.001 Inf 30.5 31.2
## 1 S 31.1 0.1698 183 <0.001 Inf 30.7 31.4
##
## Type: response
##
##
## --- VERIFICAÇÃO DE CONSISTÊNCIA ---
## Teste Wald (robusto) de diferenças entre regiao:
##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b1=b2 -0.37703 0.0176 -21.37 <0.001 334.3 -0.4116 -0.3425
## b1=b3 -0.21772 0.0159 -13.66 <0.001 138.8 -0.2490 -0.1865
## b1=b4 -0.00311 0.0155 -0.20 0.841 0.2 -0.0335 0.0273
## b1=b5 0.06289 0.0130 4.83 <0.001 19.5 0.0374 0.0884
## b2=b3 0.15931 0.0155 10.30 <0.001 80.2 0.1290 0.1896
## b2=b4 0.37392 0.0204 18.37 <0.001 247.8 0.3340 0.4138
## b2=b5 0.43992 0.0200 22.00 <0.001 353.9 0.4007 0.4791
## b3=b4 0.21461 0.0189 11.35 <0.001 96.8 0.1776 0.2517
## b3=b5 0.28062 0.0178 15.74 <0.001 183.1 0.2457 0.3155
## b4=b5 0.06600 0.0165 4.00 <0.001 13.9 0.0336 0.0984
## b1=0 -4.87647 0.1842 -26.48 <0.001 510.9 -5.2374 -4.5155
## b2=0 -4.49944 0.1712 -26.28 <0.001 503.4 -4.8350 -4.1639
## b3=0 -4.65874 0.1771 -26.31 <0.001 504.4 -5.0058 -4.3117
## b4=0 -4.87336 0.1844 -26.42 <0.001 508.7 -5.2348 -4.5119
## b5=0 -4.93936 0.1864 -26.50 <0.001 511.5 -5.3047 -4.5740Urbano ou Rural
# Por área urbana
result_area <- analyze_and_plot_marginal_effects(glm_base_ocupados_final, "recebeu_bf",
"area_urbana", vcov_rob_ocupados_3, "Horas Preditas", "por Área Urbana", comparison_type = "slopes")##
## === EFEITO DO BOLSA FAMÍLIA POR AREA_URBANA ===
##
## --- EFEITOS MARGINAIS ---
##
## area_urbana Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## 1 -4.80 0.181 -26.5 <0.001 510.1 -5.16 -4.45
## 2 -4.59 0.175 -26.3 <0.001 503.9 -4.94 -4.25
##
## Term: recebeu_bf
## Type: response
## Comparison: 1 - 0
##
##
## --- VALORES PREDITOS ---
##
## recebeu_bf area_urbana Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## 0 1 38.8 0.0384 1008 <0.001 Inf 38.7 38.8
## 0 2 37.2 0.0805 462 <0.001 Inf 37.1 37.4
## 1 1 30.1 0.1526 197 <0.001 Inf 29.8 30.4
## 1 2 29.6 0.1551 191 <0.001 Inf 29.3 29.9
##
## Type: response
##
##
## --- VERIFICAÇÃO DE CONSISTÊNCIA ---
## Teste Wald (robusto) de diferenças entre area_urbana:
##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b1=b2 -0.206 0.013 -15.8 <0.001 185.0 -0.232 -0.181
## b1=0 -4.801 0.181 -26.5 <0.001 510.1 -5.157 -4.446
## b2=0 -4.595 0.175 -26.3 <0.001 503.9 -4.937 -4.252Posição na ocupação
# Por posição na ocupação
result_posicao <- analyze_and_plot_marginal_effects(glm_base_ocupados_final, "recebeu_bf",
"pos_ocup", vcov_rob_ocupados_3, "Horas Preditas", "por Posição na Ocupação",
comparison_type = "slopes")##
## === EFEITO DO BOLSA FAMÍLIA POR POS_OCUP ===
##
## --- EFEITOS MARGINAIS ---
##
## pos_ocup Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## Com carteira -5.18 0.195 -26.5 <0.001 513.4 -5.56 -4.80
## Conta própria -4.60 0.175 -26.3 <0.001 505.7 -4.94 -4.26
## Outros -4.62 0.174 -26.5 <0.001 512.2 -4.97 -4.28
## Sem carteira -4.34 0.166 -26.2 <0.001 500.6 -4.66 -4.02
##
## Term: recebeu_bf
## Type: response
## Comparison: 1 - 0
##
##
## --- VALORES PREDITOS ---
##
## recebeu_bf pos_ocup Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## 0 Com carteira 41.6 0.0461 902 <0.001 Inf 41.5 41.7
## 0 Conta própria 37.1 0.0743 500 <0.001 Inf 37.0 37.3
## 0 Outros 37.1 0.1145 324 <0.001 Inf 36.9 37.3
## 0 Sem carteira 35.1 0.0758 463 <0.001 Inf 35.0 35.3
## 1 Com carteira 33.7 0.1761 191 <0.001 Inf 33.3 34.0
## 1 Conta própria 30.3 0.1600 189 <0.001 Inf 30.0 30.6
## 1 Outros 29.9 0.1818 165 <0.001 Inf 29.6 30.3
## 1 Sem carteira 29.0 0.1495 194 <0.001 Inf 28.7 29.3
##
## Type: response
##
##
## --- VERIFICAÇÃO DE CONSISTÊNCIA ---
## Teste Wald (robusto) de diferenças entre pos_ocup:
##
## Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## b1=b2 -0.5778 0.0233 -24.8 <0.001 449.5 -0.6234 -0.5322
## b1=b3 -0.5558 0.0262 -21.2 <0.001 328.7 -0.6072 -0.5044
## b1=b4 -0.8395 0.0316 -26.5 <0.001 513.2 -0.9015 -0.7775
## b2=b3 0.0221 0.0170 1.3 0.194 2.4 -0.0112 0.0554
## b2=b4 -0.2616 0.0159 -16.4 <0.001 199.3 -0.2928 -0.2304
## b3=b4 -0.2837 0.0190 -14.9 <0.001 165.3 -0.3209 -0.2465
## b1=0 -5.1797 0.1951 -26.5 <0.001 513.4 -5.5621 -4.7973
## b2=0 -4.6019 0.1747 -26.3 <0.001 505.7 -4.9442 -4.2595
## b3=0 -4.6239 0.1744 -26.5 <0.001 512.2 -4.9657 -4.2821
## b4=0 -4.3402 0.1656 -26.2 <0.001 500.6 -4.6648 -4.0157Resultados
Margem Extensiva
Visão geral
Os resultados da análise da margem extensiva mostram como o recebimento do Bolsa Família afeta a probabilidade de participação na força de trabalho.
Principais achados:
Efeito médio geral: O Bolsa Família aumenta significativamente a probabilidade de não participar da força de trabalho. O programa está associado a um maior risco de inatividade laboral entre os beneficiários.
Heterogeneidade por características:
Por sexo: O efeito do Bolsa Família é substancialmente maior para mulheres (4,49 pontos percentuais) comparado aos homens (3,0 pontos percentuais). Esta diferença de 1,49 p.p. é estatisticamente significativa (p<0,001), indicando que as mulheres beneficiárias têm maior propensão a não participar da força de trabalho mediante recebimento de PBF, possivelmente devido às particularidades da inserção da mulher na economia - como maiores obrigações com cuidados e afezeres domésticos. Há na própria PNADc um caderno que trata das horas gastas com cuidados, que eventualmente poderia ser utilizado em interface com a metodologia deste estudo para melhor compreender o efeito do recebimento de PBF sobre as mulheres.
Por idade: O efeito varia significativamente por quartil de idade. É maior no segundo quartil (5,64 p.p.), seguido pelo terceiro quartil (4,08 p.p.), e menor no quarto quartil de idade (1,27 p.p.). Isso sugere que o programa tem maior impacto sobre a participação laboral de pessoas em idade economicamente mais ativa inicial, diminuindo entre os mais velhos. Eventualmente poderia-se testar segmentar a idade da amostra por outros critérios que não quartis, para melhor identificar as heterogeneidades.
Por raça: Observa-se heterogeneidade racial importante. O efeito é maior entre “Outros” (4,69 p.p.) e pretos e pardos (4,44 p.p.) comparado aos brancos (2,47 p.p.). As diferenças entre pretos/pardos e brancos (-1,97 p.p.) e entre “Outros” e brancos (-2,22 p.p.) são estatisticamente significativas, indicando maior responsividade do programa entre grupos raciais historicamente desfavorecidos.
Por escolaridade: Há um gradiente claro por nível educacional. O efeito é maior para pessoas com até 11 anos de estudo (5,08 p.p.), intermediário para aqueles com 12 a 15 anos (3,33 p.p.), e menor para pessoas com 16 anos ou mais de estudo (1,30 p.p.). Todas as diferenças são estatisticamente significativas, mostrando que o programa afeta mais fortemente a participação laboral de pessoas com menor escolaridade.
Por condição no domicílio: Filhos/enteados/netos apresentam o maior efeito (6,51 p.p.), seguidos pelos responsáveis/cônjuges (2,86 p.p.), e pais/sogros/avós (1,35 p.p.). Isso sugere que jovens membros da família são mais propensos a deixar a força de trabalho quando a família recebe o benefício, eventualmente para se dedicar exclusivamente a atividades não remuneradas como os estudos. Isso pode indicar um possível impacto positivo do recebimento de PBF sobre a formação de capital humano do domicílio beneficiário, tema que pode ser melhor investigado futuramente.
Por renda domiciliar per capita: Existe um gradiente muito acentuado por renda. O efeito é dramático no primeiro quartil (9,62 p.p.), diminuindo progressivamente no segundo (4,39 p.p.), terceiro (1,95 p.p.) e quarto quartis (0,92 p.p.). Todas as diferenças entre quartis são estatisticamente significativas, demonstrando que o programa tem impacto muito maior sobre famílias mais pobres. Há que se considerar, contudo, que a maior parte dos beneficiários do PBF está no primeiro quartil.
Por região: Há variação regional substancial. O Nordeste apresenta o maior efeito (6,43 p.p.), seguido pelo Norte (4,17 p.p.), Centro-Oeste (2,40 p.p.), Sudeste (2,20 p.p.) e Sul (1,66 p.p.). As diferenças regionais são todas estatisticamente significativas, refletindo diferentes estruturas de mercado de trabalho e necessidades socioeconômicas regionais.
Por área urbana/rural: O efeito é significativamente maior em áreas rurais (7,22 p.p.) comparado às urbanas (2,52 p.p.), com diferença de 4,70 p.p. (p<0,001). Isso pode refletir menores oportunidades de trabalho formal em áreas rurais e maior dependência de transferências governamentais.
Interações significativas: O modelo inclui interações entre recebimento do Bolsa Família e área urbana/rural, bem como com quartis de idade, todas estatisticamente significativas, reforçando a heterogeneidade dos efeitos.
Interpretação
Os resultados da margem extensiva revelam que o Bolsa Família está associado a uma redução na participação na força de trabalho, mas com efeitos altamente heterogêneos. Os grupos mais vulneráveis socioeconomicamente (mulheres, pessoas com menor escolaridade, famílias mais pobres, residentes no Nordeste e áreas rurais) apresentam maior responsividade ao programa no que diz respeito ao impacto na participação na força de trabalho.
Esta heterogeneidade pode refletir:
Efeito renda: Para famílias muito pobres, o benefício pode proporcionar segurança suficiente para permitir que alguns membros se dediquem a outras atividades
Qualidade dos empregos: Em contextos de trabalhos precários e mal remunerados, o benefício pode tornar o “reservation wage” mais alto, o que teria impacto especialmente relevante em localidades onde os trabalhos ofertados possuem baixa remuneração
Custos de oportunidade: Para grupos com menores salários potenciais, o custo de oportunidade de não trabalhar é menor.
Condicionalidades: As condicionalidades educacionais podem incentivar que jovens permaneçam estudando ao invés de trabalhar.
Cuidado familiar: Especialmente para mulheres, o benefício pode facilitar a dedicação ao cuidado de crianças e idosos. Esse impacto é especialmente relevante considerando a tendência de envelhecimento populacional no Brasil.
Diagnósticos e Validação
A validação do modelo logístico para a margem extensiva passou por uma série de testes diagnósticos abrangentes para verificar a adequação da especificação e a robustez dos resultados.
Ajuste do modelo:
O modelo apresenta medidas de ajuste que parecem adequadas para um modelo logístico com dados individuais. O Pseudo-R² de McFadden foi de 0,196, o de Cox & Snell de 0,094, e o de Nagelkerke (Cragg-Uhler) de 0,238. Esses valores parecem satisfatórios pois se referem a um modelo de participação no mercado de trabalho, onde a variabilidade individual é naturalmente alta. O modelo também apresenta bom ajuste segundo o teste de Hosmer-Lemeshow (p-valor = 0,268 > 0,05), indicando que o modelo se ajusta adequadamente aos dados.
Heterocedasticidade:
Todos os testes de heterocedasticidade rejeitaram a hipótese nula de homocedasticidade. O teste de Breusch-Pagan (BP = 20.664, p<0,001), Goldfeld-Quandt (GQ = 6,753, p<0,001) e White (estatística = 241,22, p<0,001) indicaram presença de heterocedasticidade. Este resultado era esperado em dados de corte transversal com grande heterogeneidade individual, e foi adequadamente tratado através do uso de erros padrão robustos (HC1) nas estimações dos efeitos marginais.
Multicolinearidade:
Os Fatores de Inflação da Variância (VIF) mostram ausência de
multicolinearidade severa. A maioria das variáveis apresenta VIF abaixo
de 5, que é o limiar recomendado. As variáveis com VIF mais elevados são
recebeu_bf (4,18), grupo_idade (3,26) e
recebeu_bf:grupo_idade (4,75), mas todos permanecem abaixo
do limiar crítico de 10. Isso indica que a multicolinearidade não
compromete as estimações do modelo.
Significância individual e conjunta:
Os testes de Wald confirmaram a significância estatística individual
de todas as variáveis incluídas no modelo. Todas as variáveis principais
e interações são estatisticamente significativas (p<0,001), incluindo
as interações recebeu_bf:area_urbana e
recebeu_bf:grupo_idade. O teste de razão de verossimilhança
(LRT) rejeitou fortemente a hipótese nula de que o modelo completo é
igual ao modelo apenas com intercepto (LR = 18.688, p<0,001),
confirmando a significância conjunta do modelo.
Especificação funcional:
O teste RESET de Ramsey rejeitou a hipótese nula (RESET = 2.316,2, p<0,001), sugerindo possíveis problemas de especificação funcional. Contudo, para modelos logísticos, este resultado deve ser interpretado com cautela, pois a não-linearidade é inerente à função logística. A inclusão de termos de interação no modelo ajuda a capturar parte da heterogeneidade e não-linearidade nas relações. O teste de Hosmer-Lemeshow, mais específico para modelos logísticos, não rejeitou a hipótese de bom ajuste (p = 0,268), sugerindo que a especificação é adequada.
Linearidade:
Para modelos logísticos, a linearidade refere-se à relação linear entre as variáveis explicativas e o logit da probabilidade. A verificação através dos resíduos de Pearson mostrou média próxima de zero (-0,0002) e desvio padrão de 0,9998, indicando comportamento adequado. Os valores extremos (mín: -1,69, máx: 20,46) sugerem presença de outliers.
Normalidade dos resíduos:
Como esperado para modelos logísticos, o teste de Jarque-Bera rejeitou fortemente a normalidade dos resíduos de Pearson (X² = 19.899.021, p<0,001). Esta violação é esperada e não problemática em modelos logísticos, onde os resíduos não seguem distribuição normal por construção da função de verossimilhança.
Considerações sobre endogeneidade:
Embora testes formais de endogeneidade sejam limitados para modelos logísticos sem instrumentos apropriados, há pelo menos três fontes potenciais de endogeneidade que podem afetar os resultados: (1) variáveis omitidas (características não observadas que afetam tanto o recebimento do Bolsa Família quanto a participação no mercado de trabalho), (2) causalidade reversa (a situação de emprego pode influenciar a elegibilidade para o programa), e (3) erro de medida (possível subdeclaração ou erro na declaração do recebimento do benefício).
Robustez:
Para aumentar a robustez dos resultados, utilizaram-se erros padrão robustos à heterocedasticidade em todas as estimações de efeitos marginais. A consistência dos padrões encontrados entre diferentes subgrupos, a magnitude economicamente significativa dos efeitos estimados, e a adequação geral dos testes diagnósticos reforçam a confiabilidade dos achados principais. O modelo apresenta 188.910 observações e consegue explicar adequadamente a variação na participação da força de trabalho, mesmo reconhecendo as limitações metodológicas inerentes a estudos observacionais.
Margem Intensiva
Visão Geral
Os resultados da análise da margem intensiva mostram como o recebimento do Bolsa Família afeta as horas efetivas de trabalho entre pessoas já ocupadas.
Principais achados:
Efeito médio geral: O Bolsa Família reduz significativamente as horas trabalhadas entre os ocupados. O programa está associado a uma diminuição na intensidade do trabalho entre os beneficiários já inseridos no mercado de trabalho.
Heterogeneidade por características:
Por sexo: O efeito do Bolsa Família é maior para homens (-0,96 horas) comparado às mulheres (-0,52 horas). Esta diferença de 0,44 horas é estatisticamente significativa (p<0,001), indicando que homens beneficiários reduzem mais suas horas de trabalho. Há de se investigar se isso está relacionados à diferenças no perfil se inserção ocupacional de cada gênero, podendo os homens fazerem parte de posições na ocupação, por exemplo, com maior liberdade para ajustar a jornada.
Por idade: O efeito varia significativamente por quartil de idade, sendo progressivamente maior com o avanço da idade. É menor no primeiro quartil (-0,51 horas), aumenta no segundo (-0,63 horas) e terceiro quartis (-0,76 horas), e atinge seu máximo no quarto quartil (-1,24 horas). Isso sugere que trabalhadores mais velhos têm maior capacidade ou necessidade de reduzir suas horas quando recebem o benefício.
Por raça: Observa-se pouca variação racial nos efeitos. O impacto é ligeiramente menor entre brancos (-0,67 horas) comparado a pretos e pardos (-0,85 horas) e outros (-0,86 horas). As diferenças são estatisticamente significativas, mas relativamente pequenas em magnitude, indicando que a heterogeneidade racial é menos pronunciada na margem intensiva.
Por escolaridade: Há um gradiente muito claro por nível educacional. O efeito é mais elevado para pessoas com até 11 anos de estudo (-1,47 horas), substancialmente menor para aqueles com 12 a 15 anos (-0,41 horas), e não significativo para pessoas com 16 anos ou mais de estudo (-0,29 horas, p=0,495). Isso mostra que trabalhadores com menor escolaridade, provavelmente em empregos com maior flexibilidade horária, ajustam mais suas horas de trabalho. Trabalhadores com menor escolaridade também tem um custo de oportunidade menor para reduzir suas horas de trabalho.
Por condição no domicílio: Responsáveis/cônjuges apresentam o maior efeito (-0,85 horas), seguidos por filhos/enteados/netos (-0,52 horas) e pais/sogros/avós (-0,53 horas). Isso sugere que os principais provedores da família têm maior margem para reduzir suas horas quando recebem o benefício. Também pode estar relacionado ao perfil ocupacional.
Por renda domiciliar per capita: Existe um gradiente acentuado por renda. O efeito é muito forte no primeiro quartil (-1,53 horas), diminuindo progressivamente no segundo (-0,77 horas), terceiro (-0,41 horas) e quarto quartis (-0,38 horas, não significativo). Todas as diferenças entre quartis mais baixos são estatisticamente significativas, demonstrando que famílias mais pobres ajustam mais fortemente suas horas de trabalho. Há de se considerar, de qualquer forma, que a maior parte dos beneficiários do PBF está no primeiro quartil.
Por região: Há variação regional importante. O Nordeste e Norte apresentam os maiores efeitos (-1,22 e -1,20 horas, respectivamente), seguidos pelo Sul (-0,50 horas), Sudeste (-0,50 horas) e Centro-Oeste (-0,45 horas). As diferenças regionais são estatisticamente significativas, refletindo diferentes características dos mercados de trabalho regionais.
Por área urbana/rural: O efeito é expressivamente diferente entre áreas. Em áreas rurais, o efeito é pequeno e não significativo (-0,21 horas), enquanto em áreas urbanas é muito substancial (-2,84 horas). Esta diferença de 2,63 horas (p<0,001) pode refletir maior rigidez horária no trabalho rural e maior flexibilidade nos empregos urbanos.
Por posição na ocupação: Trabalhadores com carteira assinada aumentam suas horas (+1,27 horas), enquanto trabalhadores por conta própria (-2,69 horas) e sem carteira (-2,01 horas) reduzem significativamente. O efeito positivo para trabalhadores formais pode indicar maior estabilidade e possibilidade de horas extras (além de menor margem para redução de horas por rigidez na definição das jornadas de trabalho), enquanto os informais reduzem mais facilmente sua oferta de trabalho.
Interações significativas: O modelo inclui interações entre recebimento do Bolsa Família e área urbana/rural, bem como com posição na ocupação, todas estatisticamente significativas, confirmando a heterogeneidade dos efeitos.
Interpretação
Os resultados da margem intensiva revelam que o Bolsa Família está associado a uma redução nas horas trabalhadas entre ocupados, mas com efeitos altamente heterogêneos. Os grupos com menor escolaridade, famílias mais pobres, trabalhadores informais e residentes em áreas urbanas apresentam maior responsividade ao programa no que diz respeito à redução das horas trabalhadas.
Esta heterogeneidade pode refletir:
Flexibilidade ocupacional: Trabalhadores informais e por conta própria têm maior controle sobre suas horas, permitindo ajustes mais fáceis quando recebem transferências.
Efeito renda: Para famílias de baixa renda, o benefício proporciona segurança suficiente para reduzir a quantidade de horas trabalhadas.
Qualidade dos empregos: Trabalhadores em empregos precários podem usar o benefício para reduzir horas em atividades mal remuneradas.
Contexto urbano vs. rural: A maior rigidez do trabalho rural (especialmente agrícola) limita a capacidade de ajuste das horas, enquanto o ambiente urbano oferece mais flexibilidade.
Diagnósticos e Validação
A validação do modelo binomial negativa para a margem intensiva passou por uma série de testes diagnósticos abrangentes para verificar a adequação da especificação e a robustez dos resultados.
Ajuste do modelo:
O modelo apresenta medidas de ajuste razoáveis. O Pseudo-R² de McFadden foi de 0,076, o de Cox & Snell de 0,081, e o de Nagelkerke (Cragg-Uhler) de 0,121. O Pseudo-R² de Efron, específico para modelos de contagem, foi de 0,076. Embora esses valores sejam relativamente baixos, são consistentes com modelos de dados individuais em estudos de oferta de trabalho, onde a variabilidade individual é naturalmente alta.
Heterocedasticidade:
Todos os testes de heterocedasticidade rejeitaram a hipótese nula de homocedasticidade (p<0,001). O teste de Breusch-Pagan (BP = 4.166,8), Goldfeld-Quandt (GQ = 0,735) e White (estatística = 2.330,47) indicaram presença de heterocedasticidade. Este resultado era esperado em dados de corte transversal com grande heterogeneidade individual, e foi tratado através do uso de erros padrão robustos (HC1) nas estimações dos efeitos marginais.
Multicolinearidade:
Os Fatores de Inflação da Variância (VIF) mostram ausência de
multicolinearidade severa. A maioria das variáveis apresenta VIF abaixo
de 2, com apenas as interações recebeu_bf:pos_ocup (VIF =
8,55) apresentando multicolinearidade moderada, o que é aceitável e
esperado para termos de interação. Nenhuma variável apresentou VIF
superior a 10, indicando que a multicolinearidade não compromete as
estimações.
Significância individual e conjunta:
Os testes de Wald Tipo confirmaram a significância estatística individual das variáveis. Todas as variáveis principais e interações são estatisticamente significativas (p<0,05), exceto a variável raça (p=0,169). O teste de razão de verossimilhança (LRT) rejeitou fortemente a hipótese nula de que o modelo completo é igual ao modelo apenas com intercepto (LR = 11.254,85, p<0,001), confirmando a significância conjunta do modelo.
Especificação funcional:
O teste RESET de Ramsey rejeitou a hipótese nula (RESET = 95,058, p<0,001), sugerindo possíveis problemas de especificação funcional. Contudo, para modelos binomial negativa, este teste tem limitações interpretativas, pois a não-linearidade é inerente à função de ligação logarítmica. A inclusão de interações no modelo ajuda a capturar parte da não-linearidade nas relações.
Sobredispersão:
O teste específico para sobredispersão mostrou uma razão de dispersão de 0,994 (p=0,168), indicando ausência de sobredispersão significativa. O teste cujo modelo usava a distribuição de Poisson havia apresentado sobredispersão. Assim sendo, manteve-se o modelo binomial negativa por sua maior flexibilidade e robustez.
Normalidade dos resíduos:
Como esperado para modelos de contagem, o teste de Jarque-Bera rejeitou fortemente a normalidade dos resíduos de Pearson (X² = 82.063, p<0,001). Esta violação é esperada e não problemática em modelos binomial negativa, onde os resíduos não seguem distribuição normal por construção.
Linearidade:
A verificação através dos resíduos de Pearson mostrou média próxima de zero (0,000) e desvio padrão de 1,012, indicando comportamento adequado.
Considerações sobre endogeneidade:
Embora testes formais de endogeneidade sejam limitados para modelos logísticos sem instrumentos apropriados, há pelo menos três fontes potenciais de endogeneidade que podem afetar os resultados: (1) variáveis omitidas (características não observadas que afetam tanto o recebimento do Bolsa Família quanto a participação no mercado de trabalho), (2) causalidade reversa (a situação de emprego pode influenciar a elegibilidade para o programa), e (3) erro de medida (possível subdeclaração ou erro na declaração do recebimento do benefício).
Robustez:
Para aumentar a robustez dos resultados, utilizaram-se erros padrão robustos à heterocedasticidade em todas as estimações de efeitos marginais. A consistência dos padrões encontrados across diferentes subgrupos e a magnitude economicamente significativa dos efeitos estimados reforçam a confiabilidade dos achados principais, mesmo reconhecendo as limitações identificadas nos diagnósticos.
Limitações Metodológicas
Causalidade: A natureza observacional dos dados não permite estabelecer relações causais estritas. A relação estudada pode sofrer com “causalidade reversa”.
Seleção: O recebimento do Bolsa Família não é aleatório, podendo haver viés de seleção.
Variáveis omitidas: Características não observadas podem influenciar tanto o recebimento quanto os outcomes laborais.
Desenho amostral: Os cálculos foram realizados sem a calibragem da amostra pelos pesos amostrais disponíveis nos microdados anuais da PNADc.
Endogeneidade: É provável que o modelo sofra com problema de endogeneidade, pois há fatores que influenciam tanto a elegibilidade para o PBF quanto a participação no mercado de trabalho (sexo, idade, outros). Não foram desenvolvidas soluções com base em variáveis instrumentais, de forma que não foram elaboradas alternativas para lidar com eventuais problemas de endogeneidade.
Filtragem: a opção de excluir da amostra as pessoas de 14 anos ou mais fora da força de trabalho e fora da força de trabalho potencial influencia os resultados de maneira relevante. É uma opção que, por um lado, retira registros da amostra de pessoas que não poderiam entrar na força de trabalho de qualquer forma (pessoas com doenças), também retira indivíduos que eventualmente poderiam fazê-lo (donas de casa, estudantes em tempo integral)
Especificação: o resultado do teste RESET indica que há problema de especificação no modelo, podendo significar que há variáveis no modelo cuja relação com a variável dependente não é linear ou há variáveis faltando no modelo. A análise dos resíduos de Pearson também indica possível problema em relação a linearidade e eventual presença de outliers
Não declaração do benefício: há registros na amostra de pessoas que dizem não receber PBF mas que atendem a todos os critérios de elegibilidade. E os dados calculados através da PNADc não são coerentes com os dados dos registros administrativos do PBF (são menores). Por isso, provavelmente há na amostra pessoas que alegam não receber PBF mas que recebiam. Há métodos para tentar corrigir essa distorção e esforços futuros podem envolver refazer os cálculos aqui apresentados mediante amostra corrigida.
Anexo: Código R Completo
options(width = 80)
knitr::opts_chunk$set(
echo = TRUE,
warning = FALSE,
message = FALSE,
fig.align = "center",
tidy = TRUE, tidy.opts = list(width.cutoff = 80)
)
# Lista dos pacotes utilizados
packages <- c(
"PNADcIBGE",
"dplyr",
"tidyr",
"survey",
"jtools",
"MASS",
"DataExplorer",
"lares",
"performance",
"emmeans",
"sjstats",
"srvyr",
"broom",
"modelsummary",
"margins",
"marginaleffects",
"ggplot2",
"lmtest",
"sandwich",
"AER",
"car",
"knitr",
"kableExtra",
"tseries",
"skedastic",
"whitestrap",
"aod"
)
# Carregando ou instalando+carregando
for (pkg in packages) {
if (!requireNamespace(pkg, quietly = TRUE)) {
install.packages(pkg, dependencies = TRUE)
}
library(pkg, character.only = TRUE)
}
rm(pkg)
rm(packages)
"%ni%" <- Negate("%in%")
# Function to get PNADC data
get_pnadc_data <- function(year, interview, vars) {
get_pnadc(
year = year,
interview = interview,
vars = vars,
design = FALSE,
labels = FALSE
)
}
pnadc_vars <- c(
############### Filtros
"VD4003", ## Condição em relação à FT Potencial, Pessoa na FT Potencial = 1
############### Variáveis dependentes
## Participação no Mercado de Trabalho (Margem Extensiva)
"VD4001", ## Condição em relação à FT, Pessoa na FT = 1
"VD4002", ## Condição em relação a ocupação, Ocupados = 1
## Participação no Mercado de Trabalho (Margem Intensiva)
"VD4031", ## Horas habituais de trabalho (numérico)
"VD4035", ## Horas efetivas de trabalho (numérico)
## Participação no Mercado de Trabalho (remuneração)
"VD4016", ## Rendimento habitual do trabalho principal
"VD4017", ## Rendimento efetivo do trabalho principal
"VD4019", ## Rendimento habitual de todos os trabalhos
"VD4020", ## Rendimento efetivo de todos os trabalhos
############### Variáveis explicativas
"VD4047", ## Rendimento efetivo recebido de programas sociais
"V5002A", ## Se recebeu BF
"V5002A2", ## Valor efetivamente recebido via BF
############### Variáveis de controle
"UF",
"V1022", ## Urbano = 1, Rural = 2
"V1023", ## Capital = 1, Resto da RM = 2, Resto da RIDE = 3, Resto da UF = 4
"V2007", ## Sexo, Homem = 1
"V2009", ## Idade (numérico)
"V2010", ## Cor ou raça
"VD2002", ## Condição no domicílio
"VD2003", ## Número de pessoas no domicílio (numérico)
"VD3006", ## Faixa de anos de estudos
"VD4009", ## Posição na Ocupação
"VD4010", ## Grupamento de atividade do trabalho principal
"VD4011", ## Grupamento de ocupação do trabalho principal
"VD5005", ## Rendimento domiciliar per capita
"VD5006", ## Faixa de rendimento domiciliar per capita
############### Pesos
"V1032",
sprintf("V1032%03d", seq(1:200))
)
dados_pnadc <- get_pnadc_data(year = 2024, interview = 1, vars = pnadc_vars)
cat("Dimensões da base de dados:", dim(dados_pnadc), "\n")
analyze_and_plot_marginal_effects <- function(model, variable, by_vars, vcov_matrix, outcome_label, plot_title_prefix, comparison_type = "slopes") {
cat(paste0("\n=== EFEITO DO BOLSA FAMÍLIA POR ", toupper(paste(by_vars, collapse = " e ")), " ===\n"))
# Calcular efeitos marginais
if (comparison_type == "slopes") {
effects <- avg_slopes(
model,
variables = variable,
by = by_vars,
type = "response",
vcov = vcov_matrix
)
} else if (comparison_type == "comparisons") {
effects <- comparisons(
model,
variables = variable,
by = by_vars,
type = "response",
vcov = vcov_matrix
)
} else {
stop("Erro no cálculo dos efeitos marginais.")
}
# Printar efeitos marginais
cat("\n--- EFEITOS MARGINAIS ---\n")
print(effects)
cat("\n")
# Calcular predições
preds <- avg_predictions(model,
vcov = vcov_matrix,
by = c(variable, by_vars),
type = "response"
)
# Printar predições
cat("--- VALORES PREDITOS ---\n")
print(preds)
cat("\n")
# Verificar consistência entre efeitos marginais e predições
cat("--- VERIFICAÇÃO DE CONSISTÊNCIA ---\n")
if (comparison_type == "comparisons") {
# Para comparisons, verificar se as diferenças batem
preds_wide <- preds %>%
pivot_wider(names_from = all_of(variable),
values_from = estimate,
names_prefix = paste0(variable, "_"))
if(ncol(preds_wide) >= 4) { # Tem colunas para ambos os valores
col_0 <- paste0(variable, "_0")
col_1 <- paste0(variable, "_1")
if(col_0 %in% names(preds_wide) && col_1 %in% names(preds_wide)) {
preds_wide$diff_calculada <- preds_wide[[col_1]] - preds_wide[[col_0]]
# Merge com efeitos marginais para comparação
comparison_data <- merge(preds_wide, effects, by = by_vars, suffixes = c("_pred", "_effect"))
cat("Comparação entre diferenças calculadas e efeitos marginais:\n")
comparison_subset <- comparison_data %>%
select(all_of(by_vars), diff_calculada, estimate) %>%
mutate(diferenca_absoluta = abs(diff_calculada - estimate))
print(comparison_subset)
cat("\n")
}
}
}
# Gráfico de predições
plot_preds <- preds %>%
mutate(
Recebeu_BF = ifelse(.data[[variable]] == 1, "Recebe Bolsa Família", "Não Recebe")
) %>%
ggplot(aes(x = .data[[by_vars]], y = estimate,
color = Recebeu_BF, group = Recebeu_BF
)) +
geom_line(linewidth = 1.2) +
geom_point(size = 3) +
geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = Recebeu_BF),
alpha = 0.2, color = NA
) +
labs(
title = paste(outcome_label, "Preditas -", plot_title_prefix),
x = by_vars,
y = outcome_label,
color = NULL, fill = NULL
) +
scale_color_manual(values = c("#4DAF4A", "#984EA3")) +
scale_fill_manual(values = c("#4DAF4A", "#984EA3")) +
theme_minimal(base_size = 12) +
theme(legend.position = "top",
axis.text.x = element_text(angle = 45, hjust = 1))
print(plot_preds)
# Gráfico de efeitos marginais
plot_effects <- effects %>%
mutate(
Significativo = ifelse(p.value < 0.05, "Sim", "Não")
) %>%
ggplot(aes(x = .data[[by_vars]], y = estimate,
color = Significativo
)) +
geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
geom_point(size = 3) +
geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2) +
labs(
title = paste("Efeito do Bolsa Família -", plot_title_prefix),
subtitle = paste("[Recebeu vs Não Recebeu]"),
x = by_vars,
y = "Efeito Marginal",
color = "Estatisticamente\nSignificativo (p<0.05)"
) +
scale_color_manual(values = c("Sim" = "#E41A1C", "Não" = "#377EB8")) +
theme_minimal(base_size = 12) +
theme(legend.position = "top",
axis.text.x = element_text(angle = 45, hjust = 1))
print(plot_effects)
# Testes de hipóteses
num_levels <- nrow(effects)
if (num_levels > 1) {
hyp_combinations <- combn(1:num_levels, 2, simplify = FALSE)
hyp_tests <- c()
for (comb in hyp_combinations) {
hyp_tests <- c(hyp_tests, paste0("b", comb[1], " = b", comb[2]))
}
hyp_tests <- c(hyp_tests, paste0("b", 1:num_levels, " = 0"))
hypothesis_result <- effects %>%
hypotheses(
hyp_tests,
joint = FALSE
)
cat(paste0("Teste Wald (robusto) de diferenças entre ", tolower(by_vars), ":\n"))
print(hypothesis_result)
cat("\n")
}
return(list(
effects = effects,
predictions = preds,
plot_preds = plot_preds,
plot_effects = plot_effects
))
}
# Function to prepare data for extensive margin analysis
prepare_extensive_margin_data <- function(data) {
data |>
filter(VD4003 %ni% c("2")) |> # Excluindo pessoas fora da força de trabalho e fora da força de trabalho potencial. Poderia ser outro critério. A ideia é retirar quem não entraria no MT de qualquer jeito (aposentados, por exemplo)
mutate(
na_forca_trabalho = case_when(
VD4001 == "1" ~ "0",
VD4001 == "2" ~ "1",
TRUE ~ NA_character_
),
horas_ef = as.numeric(VD4035),
log_horas_ef = log(horas_ef + 0.01),
rem_ef_todos = as.numeric(VD4020),
log_rem_ef_todos = log(rem_ef_todos + 0.01),
recebeu_bf = case_when(
V5002A == "2" ~ "0",
V5002A == "1" ~ "1",
TRUE ~ NA_character_
),
sexo = as.factor(V2007),
idade = as.numeric(V2009),
log_idade = log(idade + 0.01),
raca = case_when(
V2010 == "1" ~ "B",
V2010 %in% c("2", "4") ~ "PP",
V2010 %in% c("3", "5", "9") ~ "O",
TRUE ~ NA_character_
) |> as.factor(),
anos_estudo = case_when(
VD3006 %in% c("1", "2", "3", "4") ~ "Até 11",
VD3006 %in% c("5") ~ "12 a 15",
VD3006 %in% c("6") ~ "16 ou mais",
TRUE ~ NA_character_
),
cond_dom = case_when(
VD2002 %in% c("01", "02") ~ "Responsável/Cônjuge",
VD2002 %in% c("03", "04", "05", "08", "09") ~ "Filho/Enteado/Genro/Nora/Neto/Bisneto",
VD2002 %in% c("06", "07", "11") ~ "Pai/Mãe/Padrasto/Madrasta/Sogro/Sogra/Avô/Avó",
TRUE ~ NA_character_
),
n_pessoas_dom = as.numeric(VD2003),
log_n_pessoas_dom = log(n_pessoas_dom + 0.01),
renda_dom_pc = as.numeric(VD5005),
log_renda_dom_pc = log(renda_dom_pc + 0.01),
UF = as.factor(UF),
regiao = case_when(
UF %in% c("11", "12", "13", "14", "15", "16", "17") ~ "N",
UF %in% c("21", "22", "23", "24", "25", "26", "27", "28", "29") ~ "NE",
UF %in% c("31", "32", "33", "35") ~ "SE",
UF %in% c("41", "42", "43") ~ "S",
UF %in% c("50", "51", "52", "53") ~ "CO",
TRUE ~ NA_character_
),
area_urbana = as.factor(V1022),
capital_regiao = as.factor(V1023),
pos_ocup = case_when(
VD4009 %in% c("01", "03", "05") ~ "Com carteira",
VD4009 %in% c("02", "04", "06") ~ "Sem carteira",
VD4009 %in% c("09") ~ "Conta própria",
VD4009 %in% c("07", "09", "10") ~ "Outros"
),
grup_ativ = as.factor(VD4010),
grup_ocup = as.factor(VD4011)
) |>
mutate(
na_forca_trabalho = factor(na_forca_trabalho, levels = c("0", "1")),
recebeu_bf = factor(recebeu_bf, levels = c("0", "1")),
anos_estudo = factor(anos_estudo, levels = c("Até 11", "12 a 15", "16 ou mais")),
raca = factor(raca, levels = c("B", "PP", "O")),
cond_dom = factor(cond_dom, levels = c("Responsável/Cônjuge", "Filho/Enteado/Genro/Nora/Neto/Bisneto", "Pai/Mãe/Padrasto/Madrasta/Sogro/Sogra/Avô/Avó")),
regiao = factor(regiao, levels = c("SE", "NE", "N", "CO", "S"))
) |>
mutate(
grupo_idade = cut(idade,
breaks = c(-Inf, 24, 64, Inf),
labels = c("Até 24", "25 a a 64", "65 ou mais"),
include.lowest = TRUE,
right = FALSE
),
grupo_renda_dom_pc = cut(renda_dom_pc,
breaks = c(-Inf, 218, 1412, 7060, Inf),
labels = c("Até R$ 218,00",
"R$ 218,01 a 1,00 SM",
"1,01 SM a 5,00 SM",
"Acima de 5,00 SM"),
include.lowest = TRUE,
right = FALSE
)
)
}
dados_pnadc_filtered_total <- prepare_extensive_margin_data(dados_pnadc)
vars_modelo_total <- c(
"na_forca_trabalho", "recebeu_bf", "sexo", "idade", "log_idade", "grupo_idade", "raca", "cond_dom",
"anos_estudo", "n_pessoas_dom", "log_n_pessoas_dom", "renda_dom_pc", "log_renda_dom_pc", "grupo_renda_dom_pc",
"regiao", "area_urbana", "capital_regiao"
)
dados_pnadc_filtered_total_semp <- dados_pnadc_filtered_total |>
drop_na(all_of(vars_modelo_total)) |>
dplyr::select(all_of(vars_modelo_total))
cat("Amostra para análise da margem extensiva:", nrow(dados_pnadc_filtered_total_semp), "observações\n")
# Function to fit the logistic regression model
fit_logit_model <- function(data) {
glm(na_forca_trabalho ~ recebeu_bf +
sexo + grupo_idade + raca + anos_estudo +
cond_dom + grupo_renda_dom_pc +
regiao + area_urbana,
family = binomial("logit"),
data = data
)
}
glm_base_total_final <- fit_logit_model(dados_pnadc_filtered_total_semp)
vcov_rob_total_3 <- sandwich::vcovHC(glm_base_total_final, type = "HC1")
# Sumário do modelo
summ(glm_base_total_final, exp = TRUE, confint = TRUE, vifs = TRUE, robust = "HC1")
# Função para diagnósticos completos do modelo logístico
perform_complete_diagnostics_logit <- function(model, data, vcov_model, model_name = "Modelo Logístico") {
cat("===============================================\n")
cat(paste("DIAGNÓSTICOS E VALIDAÇÃO DO", toupper(model_name), "\n"))
cat("===============================================\n\n")
# 1. ANÁLISE DO AJUSTE DO MODELO
cat("1. ANÁLISE DO AJUSTE DO MODELO\n")
cat("------------------------------\n")
# Sumário com R² do jtools
cat("Diferentes tipos de R²:\n")
tryCatch({
# Para modelos logísticos, calcular R² manualmente usando as informações do modelo
null_deviance <- model$null.deviance
residual_deviance <- model$deviance
n <- length(model$y)
# McFadden R²
mcfadden_r2 <- 1 - (residual_deviance / null_deviance)
cat("McFadden R²:", round(mcfadden_r2, 4), "\n")
# Cox & Snell R²
cox_snell_r2 <- 1 - exp((residual_deviance - null_deviance) / n)
cat("Cox & Snell R²:", round(cox_snell_r2, 4), "\n")
# Nagelkerke R²
nagelkerke_r2 <- cox_snell_r2 / (1 - exp(-null_deviance / n))
cat("Nagelkerke (Cragg-Uhler) R²:", round(nagelkerke_r2, 4), "\n")
# Tentar obter do summ() se disponível
model_summary <- summ(model, exp = FALSE, robust = "HC1")
if(is.list(model_summary) && "coeftable" %in% names(model_summary)) {
cat("Sumário do modelo obtido com sucesso\n")
}
cat("\n")
}, error = function(e) {
cat("Erro ao calcular R²:", e$message, "\n")
cat("Usando informações básicas do modelo:\n")
cat("Null Deviance:", round(model$null.deviance, 2), "\n")
cat("Residual Deviance:", round(model$deviance, 2), "\n")
cat("Diferença:", round(model$null.deviance - model$deviance, 2), "\n\n")
})
# Outros indicadores de ajuste e performance do modelo
cat("Outros indicadores de ajuste e performance do modelo:\n")
model_perf <- model_performance(model)
print(model_perf)
cat("\n")
# 2. HETEROCEDASTICIDADE
cat("2. HETEROCEDASTICIDADE\n")
cat("----------------------\n")
# Breusch-Pagan
cat("Teste de Breusch-Pagan:\n")
bp_test <- bptest(model)
print(bp_test)
cat("Interpretação: p-valor < 0.05 indica presença de heterocedasticidade\n\n")
# Goldfeld-Quandt (usando uma variável contínua como ordenação)
# Para modelos logísticos, valores ajustados
cat("Teste de Goldfeld-Quandt:\n")
tryCatch({
gq_test <- gqtest(model, order.by = fitted(model))
print(gq_test)
cat("Interpretação: p-valor < 0.05 indica presença de heterocedasticidade\n\n")
}, error = function(e) {
cat("Não foi possível executar o teste GQ para modelo logístico\n")
cat("Motivo: Este teste é mais adequado para modelos lineares\n\n")
})
# White Test (usando pacote skedastic)
cat("Teste de White:\n")
tryCatch({
# Teste de White
white_test <- skedastic::white(model)
print(white_test)
cat("Interpretação: p-valor < 0.05 indica presença de heterocedasticidade\n\n")
}, error = function(e) {
tryCatch({
# Fallback: tentar whitestrap
white_test <- whitestrap::white_test(model)
print(white_test)
cat("Interpretação: p-valor < 0.05 indica presença de heterocedasticidade\n\n")
}
)
})
# 3. MULTICOLINEARIDADE
cat("3. MULTICOLINEARIDADE\n")
cat("---------------------\n")
# VIF usando car
cat("Fatores de Inflação da Variância (VIF):\n")
vif_values <- car::vif(model)
print(vif_values)
cat("\nInterpretação:\n")
cat("- VIF < 5: Não há problema de multicolinearidade\n")
cat("- VIF entre 5-10: Multicolinearidade moderada\n")
cat("- VIF > 10: Multicolinearidade severa\n\n")
# Identificar variáveis problemáticas
if(any(vif_values > 10)) {
cat("ATENÇÃO: Variáveis com VIF > 10 (multicolinearidade severa):\n")
high_vif <- vif_values[vif_values > 10]
print(high_vif)
cat("\n")
}
# 4. SIGNIFICÂNCIA INDIVIDUAL
cat("4. SIGNIFICÂNCIA INDIVIDUAL\n")
cat("---------------------------\n")
# Teste de Wald Tipo II
cat("Teste de Wald (Tipo II):\n")
wald_test_2 <- car::Anova(model, type = "II", test.statistic = "Wald", white.adjust = "hc1")
print(wald_test_2)
cat("Interpretação: p-valor < 0.05 indica que a variável é estatisticamente significativa\n\n")
# Teste LRT Tipo II
cat("Teste de Razão de Verossimilhança (Tipo II):\n")
lr_test_2 <- car::Anova(model, type = "II", test.statistic = "LR", white.adjust = "hc1")
print(lr_test_2)
cat("Interpretação: p-valor < 0.05 indica que a variável é estatisticamente significativa\n\n")
# Teste de Wald Tipo III
cat("Teste de Wald (Tipo III):\n")
wald_test_3 <- car::Anova(model, type = "III", test.statistic = "Wald", white.adjust = "hc1")
print(wald_test_3)
cat("Interpretação: p-valor < 0.05 indica que a variável é estatisticamente significativa\n\n")
# Teste LRT Tipo III
cat("Teste de Razão de Verossimilhança (Tipo III):\n")
lr_test_3 <- car::Anova(model, type = "III", test.statistic = "LR", white.adjust = "hc1")
print(lr_test_3)
cat("Interpretação: p-valor < 0.05 indica que a variável é estatisticamente significativa\n\n")
# 5. SIGNIFICÂNCIA CONJUNTA
cat("5.4 Teste de Wald global\n")
cat("---------------------------------------\n")
tryCatch({
# Teste de Wald para todos os coeficientes exceto intercepto
n_coef <- length(coef(model))
if (n_coef > 1) {
terms_to_test <- 2:n_coef
wald_aod <- wald.test(b = coef(model), Sigma = vcov_model, Terms = terms_to_test)
print(wald_aod)
# Interpretação
p_val <- wald_aod$result$chi2[3]
cat("\nInterpretação:")
if (p_val < 0.001) {
cat(" Modelo altamente significativo (p < 0.001) ***\n")
} else if (p_val < 0.01) {
cat(" Modelo muito significativo (p < 0.01) **\n")
} else if (p_val < 0.05) {
cat(" Modelo significativo (p < 0.05) *\n")
} else {
cat(" Modelo não significativo (p >= 0.05)\n")
}
} else {
cat("Modelo tem apenas intercepto - teste não aplicável\n")
}
cat("\n")
}, error = function(e) {
cat("Erro no teste AOD/Wald:", e$message, "\n\n")
})
# 6. ESPECIFICAÇÃO DO MODELO
cat("6. ESPECIFICAÇÃO DO MODELO\n")
cat("--------------------------\n")
# RESET Test
cat("Teste RESET (Ramsey) - Termos quadráticos e cúbicos:\n")
tryCatch({
reset_test <- resettest(model, power = 2:3, type = "fitted")
print(reset_test)
cat("Interpretação: p-valor < 0.05 indica problemas de especificação funcional\n\n")
}, error = function(e) {
cat("Teste RESET não executado para modelo logístico\n")
cat("Motivo: Este teste é mais adequado para modelos lineares\n\n")
})
# Teste de Hosmer-Lemeshow (específico para modelos logísticos)
cat("Teste de Hosmer-Lemeshow (Goodness-of-Fit para modelos logísticos):\n")
hosmer_test <- performance_hosmer(model)
print(hosmer_test)
cat("Interpretação: p-valor > 0.05 indica bom ajuste do modelo\n\n")
# 7. LINEARIDADE
cat("7. LINEARIDADE (PARA MODELOS LOGÍSTICOS)\n")
cat("----------------------------------------\n")
cat("Para modelos logísticos, a linearidade se refere à relação linear entre\n")
cat("as variáveis explicativas e o logit da probabilidade.\n")
cat("Verificação através dos resíduos de Pearson:\n\n")
# Resíduos de Pearson
pearson_resid <- residuals(model, type = "pearson")
cat("Estatísticas dos Resíduos de Pearson:\n")
cat("Média:", round(mean(pearson_resid), 4), "\n")
cat("Desvio Padrão:", round(sd(pearson_resid), 4), "\n")
cat("Mín:", round(min(pearson_resid), 4), "Máx:", round(max(pearson_resid), 4), "\n\n")
# 8. NORMALIDADE DOS RESÍDUOS
cat("8. NORMALIDADE DOS RESÍDUOS\n")
cat("---------------------------\n")
cat("NOTA: Para modelos logísticos, os resíduos não seguem distribuição normal.\n")
cat("Os testes abaixo são informativos, mas não são requisitos do modelo.\n\n")
# Teste de Shapiro-Wilk (apenas para amostras pequenas)
if(length(pearson_resid) <= 5000) {
cat("Teste de Shapiro-Wilk (Resíduos de Pearson):\n")
shapiro_test <- shapiro.test(pearson_resid)
print(shapiro_test)
cat("\n")
} else {
cat("Teste de Shapiro-Wilk não executado (amostra muito grande)\n\n")
}
# Teste de Jarque-Bera
cat("Teste de Jarque-Bera (Resíduos de Pearson):\n")
jb_test <- jarque.bera.test(pearson_resid)
print(jb_test)
cat("Interpretação: Para modelos logísticos, rejeição da normalidade é esperada\n\n")
# 9. ENDOGENEIDADE (COM DEBUG)
cat("9. ENDOGENEIDADE\n")
cat("-----------------\n")
cat("Análise de correlação entre resíduos e variáveis explicativas\n")
cat("(Teste indireto para detecção de endogeneidade)\n\n")
model_resid <- residuals(model, type = "pearson")
test_vars <- c("recebeu_bf", "sexo", "grupo_idade", "raca", "anos_estudo",
"cond_dom", "grupo_renda_dom_pc", "regiao", "area_urbana")
test_vars <- test_vars[test_vars %in% names(data)]
correlation_results <- data.frame(
Variavel = character(),
Correlacao = numeric(),
P_valor = numeric(),
Significativo = character(),
stringsAsFactors = FALSE
)
cat("Correlações entre resíduos e variáveis explicativas:\n")
cat("---------------------------------------------------\n")
cat("Resíduos Pearson vs Correlação Pearson:\n")
cat("---------------------------------------------------\n")
for(i in 1:length(test_vars)) {
var <- test_vars[i]
tryCatch({
# Converter variável para numérica se necessário
if(is.factor(data[[var]])) {
var_numeric <- as.numeric(data[[var]])
} else if(is.character(data[[var]])) {
var_numeric <- as.numeric(as.factor(data[[var]]))
} else {
var_numeric <- as.numeric(data[[var]])
}
valid_indices <- !is.na(var_numeric) & !is.na(model_resid)
# Testar se há observações suficientes
if(sum(valid_indices) > 10) {
# Teste de correlação (usando Spearman por robustez)
cor_test <- cor.test(model_resid[valid_indices],
var_numeric[valid_indices],
method = "pearson")
# Determinar significância
significativo <- ifelse(cor_test$p.value < 0.05, "Sim", "Não")
# Adicionar aos resultados
correlation_results <- rbind(correlation_results,
data.frame(
Variavel = var,
Correlacao = round(cor_test$estimate, 4),
P_valor = round(cor_test$p.value, 6),
Significativo = significativo
))
# Mostrar resultado
if(var == "recebeu_bf") {
cat("*** ", var, " (VARIÁVEL PRINCIPAL): ",
round(cor_test$estimate, 4),
" (p = ", round(cor_test$p.value, 6), ") ",
significativo, " ***\n", sep="")
} else {
cat(var, ": ", round(cor_test$estimate, 4),
" (p = ", round(cor_test$p.value, 4), ") ",
significativo, "\n", sep="")
}
} else {
cat(var, ": Dados insuficientes para teste\n")
}
}, error = function(e) {
cat("ERRO na variável", var, ":", e$message, "\n")
})
}
cat("---------------------------------------------------\n")
cat("Resíduos Pearson vs Correlação Spearman:\n")
cat("---------------------------------------------------\n")
for(i in 1:length(test_vars)) {
var <- test_vars[i]
tryCatch({
# Converter variável para numérica se necessário
if(is.factor(data[[var]])) {
var_numeric <- as.numeric(data[[var]])
} else if(is.character(data[[var]])) {
var_numeric <- as.numeric(as.factor(data[[var]]))
} else {
var_numeric <- as.numeric(data[[var]])
}
valid_indices <- !is.na(var_numeric) & !is.na(model_resid)
# Testar se há observações suficientes
if(sum(valid_indices) > 10) {
# Teste de correlação (usando Spearman por robustez)
cor_test <- cor.test(model_resid[valid_indices],
var_numeric[valid_indices],
method = "spearman")
# Determinar significância
significativo <- ifelse(cor_test$p.value < 0.05, "Sim", "Não")
# Adicionar aos resultados
correlation_results <- rbind(correlation_results,
data.frame(
Variavel = var,
Correlacao = round(cor_test$estimate, 4),
P_valor = round(cor_test$p.value, 6),
Significativo = significativo
))
# Mostrar resultado
if(var == "recebeu_bf") {
cat("*** ", var, " (VARIÁVEL PRINCIPAL): ",
round(cor_test$estimate, 4),
" (p = ", round(cor_test$p.value, 6), ") ",
significativo, " ***\n", sep="")
} else {
cat(var, ": ", round(cor_test$estimate, 4),
" (p = ", round(cor_test$p.value, 4), ") ",
significativo, "\n", sep="")
}
} else {
cat(var, ": Dados insuficientes para teste\n")
}
}, error = function(e) {
cat("ERRO na variável", var, ":", e$message, "\n")
})
}
model_resid <- residuals(model, type = "deviance")
cat("---------------------------------------------------\n")
cat("Resíduos Deviance vs Correlação Pearson:\n")
cat("---------------------------------------------------\n")
for(i in 1:length(test_vars)) {
var <- test_vars[i]
tryCatch({
# Converter variável para numérica se necessário
if(is.factor(data[[var]])) {
var_numeric <- as.numeric(data[[var]])
} else if(is.character(data[[var]])) {
var_numeric <- as.numeric(as.factor(data[[var]]))
} else {
var_numeric <- as.numeric(data[[var]])
}
valid_indices <- !is.na(var_numeric) & !is.na(model_resid)
# Testar se há observações suficientes
if(sum(valid_indices) > 10) {
# Teste de correlação (usando Spearman por robustez)
cor_test <- cor.test(model_resid[valid_indices],
var_numeric[valid_indices],
method = "pearson")
# Determinar significância
significativo <- ifelse(cor_test$p.value < 0.05, "Sim", "Não")
# Adicionar aos resultados
correlation_results <- rbind(correlation_results,
data.frame(
Variavel = var,
Correlacao = round(cor_test$estimate, 4),
P_valor = round(cor_test$p.value, 6),
Significativo = significativo
))
# Mostrar resultado
if(var == "recebeu_bf") {
cat("*** ", var, " (VARIÁVEL PRINCIPAL): ",
round(cor_test$estimate, 4),
" (p = ", round(cor_test$p.value, 6), ") ",
significativo, " ***\n", sep="")
} else {
cat(var, ": ", round(cor_test$estimate, 4),
" (p = ", round(cor_test$p.value, 4), ") ",
significativo, "\n", sep="")
}
} else {
cat(var, ": Dados insuficientes para teste\n")
}
}, error = function(e) {
cat("ERRO na variável", var, ":", e$message, "\n")
})
}
cat("---------------------------------------------------\n")
cat("Resíduos Deviance vs Correlação Spearman:\n")
cat("---------------------------------------------------\n")
for(i in 1:length(test_vars)) {
var <- test_vars[i]
tryCatch({
# Converter variável para numérica se necessário
if(is.factor(data[[var]])) {
var_numeric <- as.numeric(data[[var]])
} else if(is.character(data[[var]])) {
var_numeric <- as.numeric(as.factor(data[[var]]))
} else {
var_numeric <- as.numeric(data[[var]])
}
valid_indices <- !is.na(var_numeric) & !is.na(model_resid)
# Testar se há observações suficientes
if(sum(valid_indices) > 10) {
# Teste de correlação (usando Spearman por robustez)
cor_test <- cor.test(model_resid[valid_indices],
var_numeric[valid_indices],
method = "spearman")
# Determinar significância
significativo <- ifelse(cor_test$p.value < 0.05, "Sim", "Não")
# Adicionar aos resultados
correlation_results <- rbind(correlation_results,
data.frame(
Variavel = var,
Correlacao = round(cor_test$estimate, 4),
P_valor = round(cor_test$p.value, 6),
Significativo = significativo
))
# Mostrar resultado
if(var == "recebeu_bf") {
cat("*** ", var, " (VARIÁVEL PRINCIPAL): ",
round(cor_test$estimate, 4),
" (p = ", round(cor_test$p.value, 6), ") ",
significativo, " ***\n", sep="")
} else {
cat(var, ": ", round(cor_test$estimate, 4),
" (p = ", round(cor_test$p.value, 4), ") ",
significativo, "\n", sep="")
}
} else {
cat(var, ": Dados insuficientes para teste\n")
}
}, error = function(e) {
cat("ERRO na variável", var, ":", e$message, "\n")
})
}
cat("===============================================\n")
cat("FIM DOS DIAGNÓSTICOS\n")
cat("===============================================\n\n")
}
# Executar diagnósticos para o modelo da margem extensiva
perform_complete_diagnostics_logit(glm_base_total_final, dados_pnadc_filtered_total_semp, vcov_rob_total_3, "Margem Extensiva")
cat("# EFEITOS MARGINAIS - MARGEM EXTENSIVA\n")
cat("======================================\n")
# Por sexo
result_sexo_ext <- analyze_and_plot_marginal_effects(glm_base_total_final, "recebeu_bf", "sexo", vcov_rob_total_3, "Probabilidade Predita", "por Sexo", comparison_type = "comparisons")
# Por idade
result_idade_ext <- analyze_and_plot_marginal_effects(glm_base_total_final, "recebeu_bf", "grupo_idade", vcov_rob_total_3, "Probabilidade Predita", "por Quartil de Idade", comparison_type = "comparisons")
# Por raça
result_raca_ext <- analyze_and_plot_marginal_effects(glm_base_total_final, "recebeu_bf", "raca", vcov_rob_total_3, "Probabilidade Predita", "por Raça", comparison_type = "comparisons")
# Por anos_estudo
result_anos_estudo_ext <- analyze_and_plot_marginal_effects(glm_base_total_final, "recebeu_bf", "anos_estudo", vcov_rob_total_3, "Probabilidade Predita", "por Faixa de Anos de Estudos", comparison_type = "comparisons")
# Por cond_dom
result_cond_dom_ext <- analyze_and_plot_marginal_effects(glm_base_total_final, "recebeu_bf", "cond_dom", vcov_rob_total_3, "Probabilidade Predita", "por Condição no Domicílio", comparison_type = "comparisons")
# Por renda
result_renda_ext <- analyze_and_plot_marginal_effects(glm_base_total_final, "recebeu_bf", "grupo_renda_dom_pc", vcov_rob_total_3, "Probabilidade Predita", "por Quartil de Renda per Capita", comparison_type = "comparisons")
# Por região
result_regiao_ext <- analyze_and_plot_marginal_effects(glm_base_total_final, "recebeu_bf", "regiao", vcov_rob_total_3, "Probabilidade Predita", "por Região", comparison_type = "comparisons")
# Por área urbana
result_area_ext <- analyze_and_plot_marginal_effects(glm_base_total_final, "recebeu_bf", "area_urbana", vcov_rob_total_3, "Probabilidade Predita", "por Área Urbana", comparison_type = "comparisons")
# Function to prepare data for intensive margin analysis
prepare_intensive_margin_data <- function(data) {
data |>
filter(VD4002 %in% "1") |>
mutate(
na_forca_trabalho = case_when(
VD4001 == "1" ~ "0",
VD4001 == "2" ~ "1",
TRUE ~ NA_character_
),
horas_ef = as.numeric(VD4035),
log_horas_ef = log(horas_ef + 0.01),
rem_ef_todos = as.numeric(VD4020),
log_rem_ef_todos = log(rem_ef_todos + 0.01),
recebeu_bf = case_when(
V5002A == "2" ~ "0",
V5002A == "1" ~ "1",
TRUE ~ NA_character_
),
sexo = as.factor(V2007),
idade = as.numeric(V2009),
log_idade = log(idade + 0.01),
raca = case_when(
V2010 == "1" ~ "B",
V2010 %in% c("2", "4") ~ "PP",
V2010 %in% c("3", "5", "9") ~ "O",
TRUE ~ NA_character_
) |> as.factor(),
anos_estudo = case_when(
VD3006 %in% c("1", "2", "3", "4") ~ "Até 11",
VD3006 %in% c("5") ~ "12 a 15",
VD3006 %in% c("6") ~ "16 ou mais",
TRUE ~ NA_character_
),
cond_dom = case_when(
VD2002 %in% c("01", "02") ~ "Responsável/Cônjuge",
VD2002 %in% c("03", "04", "05", "08", "09") ~ "Filho/Enteado/Genro/Nora/Neto/Bisneto",
VD2002 %in% c("06", "07", "11") ~ "Pai/Mãe/Padrasto/Madrasta/Sogro/Sogra/Avô/Avó",
TRUE ~ NA_character_
),
n_pessoas_dom = as.numeric(VD2003),
log_n_pessoas_dom = log(n_pessoas_dom + 0.01),
renda_dom_pc = as.numeric(VD5005),
log_renda_dom_pc = log(renda_dom_pc + 0.01),
UF = as.factor(UF),
regiao = case_when(
UF %in% c("11", "12", "13", "14", "15", "16", "17") ~ "N",
UF %in% c("21", "22", "23", "24", "25", "26", "27", "28", "29") ~ "NE",
UF %in% c("31", "32", "33", "35") ~ "SE",
UF %in% c("41", "42", "43") ~ "S",
UF %in% c("50", "51", "52", "53") ~ "CO",
TRUE ~ NA_character_
),
area_urbana = as.factor(V1022),
capital_regiao = as.factor(V1023),
pos_ocup = case_when(
VD4009 %in% c("01", "03", "05") ~ "Com carteira",
VD4009 %in% c("02", "04", "06") ~ "Sem carteira",
VD4009 %in% c("09") ~ "Conta própria",
VD4009 %in% c("07", "09", "10") ~ "Outros"
),
grup_ativ = as.factor(VD4010),
grup_ocup = as.factor(VD4011)
) |>
mutate(
na_forca_trabalho = factor(na_forca_trabalho, levels = c("0", "1")),
recebeu_bf = factor(recebeu_bf, levels = c("0", "1")),
anos_estudo = factor(anos_estudo, levels = c("Até 11", "12 a 15", "16 ou mais")),
raca = factor(raca, levels = c("B", "PP", "O")),
cond_dom = factor(cond_dom, levels = c("Responsável/Cônjuge", "Filho/Enteado/Genro/Nora/Neto/Bisneto", "Pai/Mãe/Padrasto/Madrasta/Sogro/Sogra/Avô/Avó")),
regiao = factor(regiao, levels = c("SE", "NE", "N", "CO", "S"))
) |>
mutate(
grupo_idade = cut(idade,
breaks = c(-Inf, 24, 64, Inf),
labels = c("Até 24", "25 a a 64", "65 ou mais"),
include.lowest = TRUE,
right = FALSE
),
grupo_renda_dom_pc = cut(renda_dom_pc,
breaks = c(-Inf, 218, 1412, 7060, Inf),
labels = c("Até R$ 218,00",
"R$ 218,01 a 1,00 SM",
"1,01 SM a 5,00 SM",
"Acima de 5,00 SM"),
include.lowest = TRUE,
right = FALSE
)
)
}
dados_pnadc_filtered_ocupados <- prepare_intensive_margin_data(dados_pnadc)
vars_modelo_ocupados <- c(
"horas_ef", "log_horas_ef", "recebeu_bf", "sexo", "idade", "log_idade", "grupo_idade", "raca",
"anos_estudo", "n_pessoas_dom", "log_n_pessoas_dom", "renda_dom_pc", "grupo_renda_dom_pc", "log_renda_dom_pc", "cond_dom",
"regiao", "area_urbana", "capital_regiao",
"pos_ocup", "grup_ativ", "grup_ocup", "rem_ef_todos", "log_rem_ef_todos"
)
dados_pnadc_filtered_ocupados_semp <- dados_pnadc_filtered_ocupados |>
drop_na(all_of(vars_modelo_ocupados)) |>
dplyr::select(all_of(vars_modelo_ocupados))
cat("Amostra para análise da margem intensiva:", nrow(dados_pnadc_filtered_ocupados_semp), "observações\n")
# Function to fit the negative binomial model
fit_nb_model <- function(data) {
glm.nb(horas_ef ~ recebeu_bf +
sexo + grupo_idade + anos_estudo + raca +
cond_dom +
grupo_renda_dom_pc +
regiao + area_urbana +
pos_ocup,
data = data
)
}
glm_base_ocupados_final <- fit_nb_model(dados_pnadc_filtered_ocupados_semp)
vcov_rob_ocupados_3 <- sandwich::vcovHC(glm_base_ocupados_final, type = "HC1")
# Sumário do modelo
summ(glm_base_ocupados_final, exp = TRUE, confint = TRUE, vif = TRUE, robust = "HC1")
# Função para diagnósticos completos do modelo binomial negativa
perform_complete_diagnostics_nb <- function(model, data, vcov_model, model_name = "Modelo Binomial Negativa") {
cat("===============================================\n")
cat(paste("DIAGNÓSTICOS E VALIDAÇÃO DO", toupper(model_name), "\n"))
cat("===============================================\n\n")
# 1. ANÁLISE DO AJUSTE DO MODELO
cat("1. ANÁLISE DO AJUSTE DO MODELO\n")
cat("------------------------------\n")
# Sumário com medidas do jtools
cat("Medidas de Ajuste (equivalentes ao R²):\n")
tryCatch({
# Para modelos binomial negativa, calcular pseudo-R² e outras medidas
null_deviance <- model$null.deviance
residual_deviance <- model$deviance
n <- length(model$y)
# McFadden Pseudo-R² (análogo ao R² para GLM)
mcfadden_r2 <- 1 - (residual_deviance / null_deviance)
cat("McFadden Pseudo-R²:", round(mcfadden_r2, 4), "\n")
# Cox & Snell Pseudo-R²
cox_snell_r2 <- 1 - exp((residual_deviance - null_deviance) / n)
cat("Cox & Snell Pseudo-R²:", round(cox_snell_r2, 4), "\n")
# Nagelkerke Pseudo-R² (Cragg-Uhler)
nagelkerke_r2 <- cox_snell_r2 / (1 - exp(-null_deviance / n))
cat("Nagelkerke Pseudo-R² (Cragg-Uhler):", round(nagelkerke_r2, 4), "\n")
# Efron's Pseudo-R² (específico para modelos de contagem)
y_obs <- model$y
y_pred <- fitted(model)
y_mean <- mean(y_obs)
efron_r2 <- 1 - sum((y_obs - y_pred)^2) / sum((y_obs - y_mean)^2)
cat("Efron's Pseudo-R²:", round(efron_r2, 4), "\n")
# Medidas adicionais de ajuste
cat("\nOutras medidas de ajuste:\n")
cat("AIC:", round(AIC(model), 2), "\n")
cat("BIC:", round(BIC(model), 2), "\n")
cat("Log-Likelihood:", round(as.numeric(logLik(model)), 2), "\n")
# Proporção de deviance explicada
prop_dev_explained <- (null_deviance - residual_deviance) / null_deviance
cat("Proporção de Deviance Explicada:", round(prop_dev_explained, 4), "\n")
cat("\n")
}, error = function(e) {
cat("Erro ao calcular medidas de ajuste:", e$message, "\n")
cat("Usando informações básicas do modelo:\n")
cat("AIC:", round(AIC(model), 2), "\n")
cat("BIC:", round(BIC(model), 2), "\n")
cat("Log-Likelihood:", round(as.numeric(logLik(model)), 2), "\n")
cat("Deviance:", round(deviance(model), 2), "\n\n")
})
# Outros indicadores de ajuste e performance do modelo
cat("Outros indicadores de ajuste e performance do modelo:\n")
model_perf <- model_performance(model)
print(model_perf)
cat("\n")
# 2. HETEROCEDASTICIDADE
cat("2. HETEROCEDASTICIDADE\n")
cat("----------------------\n")
# Breusch-Pagan
cat("Teste de Breusch-Pagan:\n")
bp_test <- bptest(model)
print(bp_test)
cat("Interpretação: p-valor < 0.05 indica presença de heterocedasticidade\n\n")
# Goldfeld-Quandt
cat("Teste de Goldfeld-Quandt:\n")
tryCatch({
gq_test <- gqtest(model, order.by = fitted(model))
print(gq_test)
cat("Interpretação: p-valor < 0.05 indica presença de heterocedasticidade\n\n")
}, error = function(e) {
cat("Teste GQ não executado para modelo binomial negativa\n")
cat("Motivo: Limitações computacionais ou características do modelo\n\n")
})
# White Test (usando pacote skedastic)
cat("Teste de White:\n")
tryCatch({
# Teste de White
white_test <- skedastic::white(model)
print(white_test)
cat("Interpretação: p-valor < 0.05 indica presença de heterocedasticidade\n\n")
}, error = function(e) {
tryCatch({
# Fallback: tentar whitestrap
white_test <- whitestrap::white_test(model)
print(white_test)
cat("Interpretação: p-valor < 0.05 indica presença de heterocedasticidade\n\n")
}
)
})
# 3. MULTICOLINEARIDADE
cat("3. MULTICOLINEARIDADE\n")
cat("---------------------\n")
# VIF usando car
cat("Fatores de Inflação da Variância (VIF):\n")
vif_values <- car::vif(model)
print(vif_values)
cat("\nInterpretação:\n")
cat("- VIF < 5: Não há problema de multicolinearidade\n")
cat("- VIF entre 5-10: Multicolinearidade moderada\n")
cat("- VIF > 10: Multicolinearidade severa\n\n")
# Identificar variáveis problemáticas
if(any(vif_values > 10)) {
cat("ATENÇÃO: Variáveis com VIF > 10 (multicolinearidade severa):\n")
high_vif <- vif_values[vif_values > 10]
print(high_vif)
cat("\n")
}
# 4. SIGNIFICÂNCIA INDIVIDUAL
cat("4. SIGNIFICÂNCIA INDIVIDUAL\n")
cat("---------------------------\n")
# Teste de Wald Tipo II
cat("Teste de Wald (Tipo II):\n")
wald_test_2 <- car::Anova(model, type = "II", test.statistic = "Wald", white.adjust = "hc1")
print(wald_test_2)
cat("Interpretação: p-valor < 0.05 indica que a variável é estatisticamente significativa\n\n")
# Teste LRT Tipo II
cat("Teste de Razão de Verossimilhança (Tipo II):\n")
lr_test_2 <- car::Anova(model, type = "II", test.statistic = "LR", white.adjust = "hc1")
print(lr_test_2)
cat("Interpretação: p-valor < 0.05 indica que a variável é estatisticamente significativa\n\n")
# Teste de Wald Tipo III
cat("Teste de Wald (Tipo III):\n")
wald_test_3 <- car::Anova(model, type = "III", test.statistic = "Wald", white.adjust = "hc1")
print(wald_test_3)
cat("Interpretação: p-valor < 0.05 indica que a variável é estatisticamente significativa\n\n")
# Teste LRT Tipo III
cat("Teste de Razão de Verossimilhança (Tipo III):\n")
lr_test_3 <- car::Anova(model, type = "III", test.statistic = "LR", white.adjust = "hc1")
print(lr_test_3)
cat("Interpretação: p-valor < 0.05 indica que a variável é estatisticamente significativa\n\n")
# 5. SIGNIFICÂNCIA CONJUNTA
cat("5.4 Teste de Wald global\n")
cat("---------------------------------------\n")
tryCatch({
# Teste de Wald para todos os coeficientes exceto intercepto
n_coef <- length(coef(model))
if (n_coef > 1) {
terms_to_test <- 2:n_coef
wald_aod <- wald.test(b = coef(model), Sigma = vcov_model, Terms = terms_to_test)
print(wald_aod)
# Interpretação
p_val <- wald_aod$result$chi2[3]
cat("\nInterpretação:")
if (p_val < 0.001) {
cat(" Modelo altamente significativo (p < 0.001) ***\n")
} else if (p_val < 0.01) {
cat(" Modelo muito significativo (p < 0.01) **\n")
} else if (p_val < 0.05) {
cat(" Modelo significativo (p < 0.05) *\n")
} else {
cat(" Modelo não significativo (p >= 0.05)\n")
}
} else {
cat("Modelo tem apenas intercepto - teste não aplicável\n")
}
cat("\n")
}, error = function(e) {
cat("Erro no teste AOD/Wald:", e$message, "\n\n")
})
# 6. ESPECIFICAÇÃO DO MODELO
cat("6. ESPECIFICAÇÃO DO MODELO\n")
cat("--------------------------\n")
# RESET Test
cat("Teste RESET (Ramsey):\n")
tryCatch({
reset_test <- resettest(model, power = 2:3, type = "fitted")
print(reset_test)
cat("Interpretação: p-valor < 0.05 indica problemas de especificação funcional\n\n")
}, error = function(e) {
cat("Teste RESET não executado para modelo binomial negativa\n")
cat("Motivo: Limitações do teste para modelos GLM não-lineares\n\n")
})
# 7. LINEARIDADE
cat("7. LINEARIDADE (PARA MODELOS DE CONTAGEM)\n")
cat("-----------------------------------------\n")
cat("Para modelos binomial negativa, a linearidade se refere à relação\n")
cat("entre as variáveis explicativas e o log da média condicional.\n")
cat("Verificação através dos resíduos de Pearson:\n\n")
# Resíduos de Pearson
pearson_resid <- residuals(model, type = "pearson")
cat("Estatísticas dos Resíduos de Pearson:\n")
cat("Média:", round(mean(pearson_resid), 4), "\n")
cat("Desvio Padrão:", round(sd(pearson_resid), 4), "\n")
cat("Mín:", round(min(pearson_resid), 4), "Máx:", round(max(pearson_resid), 4), "\n\n")
# 8. NORMALIDADE DOS RESÍDUOS
cat("8. NORMALIDADE DOS RESÍDUOS\n")
cat("---------------------------\n")
cat("NOTA: Para modelos binomial negativa, os resíduos não seguem distribuição normal.\n")
cat("Os testes abaixo são informativos, mas não são requisitos do modelo.\n\n")
# Teste de Shapiro-Wilk (apenas para amostras pequenas)
if(length(pearson_resid) <= 5000) {
cat("Teste de Shapiro-Wilk (Resíduos de Pearson):\n")
shapiro_test <- shapiro.test(pearson_resid)
print(shapiro_test)
cat("\n")
} else {
cat("Teste de Shapiro-Wilk não executado (amostra muito grande)\n\n")
}
# Teste de Jarque-Bera
cat("Teste de Jarque-Bera (Resíduos de Pearson):\n")
jb_test <- jarque.bera.test(pearson_resid)
print(jb_test)
cat("Interpretação: Para modelos de contagem, rejeição da normalidade é esperada\n\n")
# 9. SOBREDISPERSÃO
cat("9. SOBREDISPERSÃO\n")
cat("-----------------\n")
cat("Teste específico para modelos de contagem:\n")
# Teste de sobredispersão
overdispersion_test <- check_overdispersion(model)
print(overdispersion_test)
cat("\nInterpretação:\n")
cat("- p-valor < 0.05: Sobredispersão estatisticamente significativa\n\n")
# 10. ENDOGENEIDADE (COM DEBUG)
cat("10. ENDOGENEIDADE\n")
cat("-----------------\n")
cat("Análise de correlação entre resíduos e variáveis explicativas\n")
cat("(Teste indireto para detecção de endogeneidade)\n\n")
model_resid <- residuals(model, type = "pearson")
test_vars <- c("recebeu_bf", "sexo", "grupo_idade", "raca", "anos_estudo",
"cond_dom", "grupo_renda_dom_pc", "regiao", "area_urbana")
test_vars <- test_vars[test_vars %in% names(data)]
correlation_results <- data.frame(
Variavel = character(),
Correlacao = numeric(),
P_valor = numeric(),
Significativo = character(),
stringsAsFactors = FALSE
)
cat("Correlações entre resíduos e variáveis explicativas:\n")
cat("---------------------------------------------------\n")
cat("Resíduos Pearson vs Correlação Pearson:\n")
cat("---------------------------------------------------\n")
for(i in 1:length(test_vars)) {
var <- test_vars[i]
tryCatch({
# Converter variável para numérica se necessário
if(is.factor(data[[var]])) {
var_numeric <- as.numeric(data[[var]])
} else if(is.character(data[[var]])) {
var_numeric <- as.numeric(as.factor(data[[var]]))
} else {
var_numeric <- as.numeric(data[[var]])
}
valid_indices <- !is.na(var_numeric) & !is.na(model_resid)
# Testar se há observações suficientes
if(sum(valid_indices) > 10) {
# Teste de correlação (usando Spearman por robustez)
cor_test <- cor.test(model_resid[valid_indices],
var_numeric[valid_indices],
method = "pearson")
# Determinar significância
significativo <- ifelse(cor_test$p.value < 0.05, "Sim", "Não")
# Adicionar aos resultados
correlation_results <- rbind(correlation_results,
data.frame(
Variavel = var,
Correlacao = round(cor_test$estimate, 4),
P_valor = round(cor_test$p.value, 6),
Significativo = significativo
))
# Mostrar resultado
if(var == "recebeu_bf") {
cat("*** ", var, " (VARIÁVEL PRINCIPAL): ",
round(cor_test$estimate, 4),
" (p = ", round(cor_test$p.value, 6), ") ",
significativo, " ***\n", sep="")
} else {
cat(var, ": ", round(cor_test$estimate, 4),
" (p = ", round(cor_test$p.value, 4), ") ",
significativo, "\n", sep="")
}
} else {
cat(var, ": Dados insuficientes para teste\n")
}
}, error = function(e) {
cat("ERRO na variável", var, ":", e$message, "\n")
})
}
cat("---------------------------------------------------\n")
cat("Resíduos Pearson vs Correlação Spearman:\n")
cat("---------------------------------------------------\n")
for(i in 1:length(test_vars)) {
var <- test_vars[i]
tryCatch({
# Converter variável para numérica se necessário
if(is.factor(data[[var]])) {
var_numeric <- as.numeric(data[[var]])
} else if(is.character(data[[var]])) {
var_numeric <- as.numeric(as.factor(data[[var]]))
} else {
var_numeric <- as.numeric(data[[var]])
}
valid_indices <- !is.na(var_numeric) & !is.na(model_resid)
# Testar se há observações suficientes
if(sum(valid_indices) > 10) {
# Teste de correlação (usando Spearman por robustez)
cor_test <- cor.test(model_resid[valid_indices],
var_numeric[valid_indices],
method = "spearman")
# Determinar significância
significativo <- ifelse(cor_test$p.value < 0.05, "Sim", "Não")
# Adicionar aos resultados
correlation_results <- rbind(correlation_results,
data.frame(
Variavel = var,
Correlacao = round(cor_test$estimate, 4),
P_valor = round(cor_test$p.value, 6),
Significativo = significativo
))
# Mostrar resultado
if(var == "recebeu_bf") {
cat("*** ", var, " (VARIÁVEL PRINCIPAL): ",
round(cor_test$estimate, 4),
" (p = ", round(cor_test$p.value, 6), ") ",
significativo, " ***\n", sep="")
} else {
cat(var, ": ", round(cor_test$estimate, 4),
" (p = ", round(cor_test$p.value, 4), ") ",
significativo, "\n", sep="")
}
} else {
cat(var, ": Dados insuficientes para teste\n")
}
}, error = function(e) {
cat("ERRO na variável", var, ":", e$message, "\n")
})
}
model_resid <- residuals(model, type = "deviance")
cat("---------------------------------------------------\n")
cat("Resíduos Deviance vs Correlação Pearson:\n")
cat("---------------------------------------------------\n")
for(i in 1:length(test_vars)) {
var <- test_vars[i]
tryCatch({
# Converter variável para numérica se necessário
if(is.factor(data[[var]])) {
var_numeric <- as.numeric(data[[var]])
} else if(is.character(data[[var]])) {
var_numeric <- as.numeric(as.factor(data[[var]]))
} else {
var_numeric <- as.numeric(data[[var]])
}
valid_indices <- !is.na(var_numeric) & !is.na(model_resid)
# Testar se há observações suficientes
if(sum(valid_indices) > 10) {
# Teste de correlação (usando Spearman por robustez)
cor_test <- cor.test(model_resid[valid_indices],
var_numeric[valid_indices],
method = "pearson")
# Determinar significância
significativo <- ifelse(cor_test$p.value < 0.05, "Sim", "Não")
# Adicionar aos resultados
correlation_results <- rbind(correlation_results,
data.frame(
Variavel = var,
Correlacao = round(cor_test$estimate, 4),
P_valor = round(cor_test$p.value, 6),
Significativo = significativo
))
# Mostrar resultado
if(var == "recebeu_bf") {
cat("*** ", var, " (VARIÁVEL PRINCIPAL): ",
round(cor_test$estimate, 4),
" (p = ", round(cor_test$p.value, 6), ") ",
significativo, " ***\n", sep="")
} else {
cat(var, ": ", round(cor_test$estimate, 4),
" (p = ", round(cor_test$p.value, 4), ") ",
significativo, "\n", sep="")
}
} else {
cat(var, ": Dados insuficientes para teste\n")
}
}, error = function(e) {
cat("ERRO na variável", var, ":", e$message, "\n")
})
}
cat("---------------------------------------------------\n")
cat("Resíduos Deviance vs Correlação Spearman:\n")
cat("---------------------------------------------------\n")
for(i in 1:length(test_vars)) {
var <- test_vars[i]
tryCatch({
# Converter variável para numérica se necessário
if(is.factor(data[[var]])) {
var_numeric <- as.numeric(data[[var]])
} else if(is.character(data[[var]])) {
var_numeric <- as.numeric(as.factor(data[[var]]))
} else {
var_numeric <- as.numeric(data[[var]])
}
valid_indices <- !is.na(var_numeric) & !is.na(model_resid)
# Testar se há observações suficientes
if(sum(valid_indices) > 10) {
# Teste de correlação (usando Spearman por robustez)
cor_test <- cor.test(model_resid[valid_indices],
var_numeric[valid_indices],
method = "spearman")
# Determinar significância
significativo <- ifelse(cor_test$p.value < 0.05, "Sim", "Não")
# Adicionar aos resultados
correlation_results <- rbind(correlation_results,
data.frame(
Variavel = var,
Correlacao = round(cor_test$estimate, 4),
P_valor = round(cor_test$p.value, 6),
Significativo = significativo
))
# Mostrar resultado
if(var == "recebeu_bf") {
cat("*** ", var, " (VARIÁVEL PRINCIPAL): ",
round(cor_test$estimate, 4),
" (p = ", round(cor_test$p.value, 6), ") ",
significativo, " ***\n", sep="")
} else {
cat(var, ": ", round(cor_test$estimate, 4),
" (p = ", round(cor_test$p.value, 4), ") ",
significativo, "\n", sep="")
}
} else {
cat(var, ": Dados insuficientes para teste\n")
}
}, error = function(e) {
cat("ERRO na variável", var, ":", e$message, "\n")
})
}
cat("===============================================\n")
cat("FIM DOS DIAGNÓSTICOS\n")
cat("===============================================\n\n")
}
# Executar diagnósticos para o modelo da margem intensiva
perform_complete_diagnostics_nb(glm_base_ocupados_final, dados_pnadc_filtered_ocupados_semp, vcov_rob_ocupados_3, "Margem Intensiva")
vcov_rob_ocupados_3 <- sandwich::vcovHC(glm_base_ocupados_final, type = "HC1")
# Análise por diferentes características
cat("# EFEITOS MARGINAIS - MARGEM INTENSIVA\n")
cat("=====================================\n")
# Por sexo
result_sexo <- analyze_and_plot_marginal_effects(glm_base_ocupados_final, "recebeu_bf", "sexo", vcov_rob_ocupados_3, "Horas Preditas", "por sexo", comparison_type = "slopes")
# Por idade
result_idade <- analyze_and_plot_marginal_effects(glm_base_ocupados_final, "recebeu_bf", "grupo_idade", vcov_rob_ocupados_3, "Horas Preditas", "por Quartil de Idade", comparison_type = "slopes")
# Por raça
result_raca <- analyze_and_plot_marginal_effects(glm_base_ocupados_final, "recebeu_bf", "raca", vcov_rob_ocupados_3, "Horas Preditas", "por Raça", comparison_type = "slopes")
# Por anos_estudo
result_anos_estudo <- analyze_and_plot_marginal_effects(glm_base_ocupados_final, "recebeu_bf", "anos_estudo", vcov_rob_ocupados_3, "Horas Preditas", "por Faixa de Anos de Estudos", comparison_type = "slopes")
# Por cond_dom
result_cond_dom <- analyze_and_plot_marginal_effects(glm_base_ocupados_final, "recebeu_bf", "cond_dom", vcov_rob_ocupados_3, "Horas Preditas", "por Condição no Domicílio", comparison_type = "slopes")
# Por renda
result_renda <- analyze_and_plot_marginal_effects(glm_base_ocupados_final, "recebeu_bf", "grupo_renda_dom_pc", vcov_rob_ocupados_3, "Horas Preditas", "por Quartil de Renda per Capita", comparison_type = "slopes")
# Por região
result_regiao <- analyze_and_plot_marginal_effects(glm_base_ocupados_final, "recebeu_bf", "regiao", vcov_rob_ocupados_3, "Horas Preditas", "por Região", comparison_type = "slopes")
# Por área urbana
result_area <- analyze_and_plot_marginal_effects(glm_base_ocupados_final, "recebeu_bf", "area_urbana", vcov_rob_ocupados_3, "Horas Preditas", "por Área Urbana", comparison_type = "slopes")
# Por posição na ocupação
result_posicao <- analyze_and_plot_marginal_effects(glm_base_ocupados_final, "recebeu_bf", "pos_ocup", vcov_rob_ocupados_3, "Horas Preditas", "por Posição na Ocupação", comparison_type = "slopes")
# Todo o código utilizado na análise está disponível nos chunks acimaRelatório gerado em 2025-07-14 usando R Markdown.