Introdução

Este relatório analisa a heterogeneidade nos 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 o efeito via a margem extensiva, ou seja, na probabilidade de participação na força de trabalho. A teoria microeconômica tradicional indica que transferências governamentais implicam em diminuição da oferta de trabalho. A hipótese subjacente a este trabalho é de que este efeito existe, mas é heterogêneo entre os diversos grupos populacionais.

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", "huxtable")

# 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%")

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 245

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).

Outras Variáveis Explicativas:

  • Características individuais: sexo (V2007), idade (V2009 agrupada), raça (V2010 agrupada), escolaridade (VD3006 agrupada);
  • Características domiciliares: condição no domicílio (VD2002 agrupada), renda domiciliar per capita (VD5005, 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).
  • Variáveis de interação.

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 39, 40 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
  • Até 11 anos de estudos
  • Responsável/cônjuge no domicílio
  • Renda Domiciliar per capita de até 218
  • Região sudeste
  • Urbano

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. Ou seja, foram excluídos registros de pessoas que faziam parte da população em idade ativa, mas que não faziam parte da força de trabalho potencial. Essa é uma opção que impacta de maneira relevante os cálculos e seus impactos podem ser avaliados posteriormente.

  • Também foram excluídos da amostra qualquer registro individual com renda domiciliar per capita superior a 2 salários mínimos. Observou-se que parte expressiva dos registros de beneficiários do PBF na PNADc tinha renda domiciliar per capita efetiva superior ao critério de elegibilidade de R$ 218,00. A maior parte desses beneficiários tinha renda até 2SM.

# 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, 39, 64, Inf),
                    labels = c("Até 24", "25 a 39", "40 a 64", "65 ou mais"),
                    include.lowest = TRUE,
                    right = FALSE
  ),
  grupo_renda_dom_pc = cut(renda_dom_pc,
                           breaks = c(-Inf, 218, 1412, 2824, 7060, Inf),
                           labels = c("Até R$ 218,00", 
                                     "R$ 218,01 a 1,00 SM",
                                     "1,01 SM a 2,00 SM",
                                     "2,01 SM a 5,00 SM",
                                     "Acima de 5,00 SM"),
                           include.lowest = TRUE,
                           right = FALSE
  )
) |> 
    filter(grupo_renda_dom_pc %in% c("Até R$ 218,00", "R$ 218,01 a 1,00 SM", "1,01 SM a 2,00 SM"))
}

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: 147752 observações

Modelo Escolhido

Tipo de modelagem:

  • Como a variável dependente é binária (o indivíduo faz parte da força de trabalho ou não), optou-se por utilizar um Modelo Linear Generalizado do tipo binomial/Logit. O Logit foi escolhido em detrimento de outras opções binárias, como o Probit, por facilitar a interpretação. Além disso, o modelo estimado com logit alcançou Pseudo-R² e BIC superior ao estimado com probit.

Pesos amostrais:

  • Para facilitar as testagens de robustez 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. De todo modo, cabe mencionar que a saída dos modelos calculados 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 opção teórica. Foram selecionados fatores de ordem sociodemográficos, domiciliares e geográficos que poderiam ser relevantes na determinação da chance de participação do mercado de trabalho. A ideia foi fazer um modelo o mais detalhado possível, para evitar problemas com variáveis omitidas;

  • Foram adicionadas variáveis de interação entre cada um dos controles e a variável explicativa principal (recebeu_bf). Como o objeto central da análise é justamente a heterogeneidade destes efeitos, foram adicionadas o máximo de interações possíveis. Através de um teste de significância individual das variáveis, feitos via Wald e via LRT, tipo II ou III, percebeu-se que as interações de recebeu_bf com sexo, raça e região não eram significativas. Por isso, elas foram excluídas do modelo.

  • Em seguida, analisou-se os outliers do modelo para considerar novos aprimoramentos. Foram utilizados os resíduos padronizados (função rstandard(), para lidar melhor com a heterocedasticidade do modelo) do tipo “deviance” (baseados na função de máxima verossimilhança, mais adequada para um modelo do tipo logit). Assim, identificou-se que os outliers tinham predominância nos grupos: pessoas que não recebiam bolsa família; 40-64 anos; brancos; 12-15 anos de estudos; e de áreas urbanas;

  • A partir disso, foram testadas novas interações dois-a-dois e três-a-três das variáveis de controle. Dentre as que mais contribuíram para o poder explicativo do modelo, a de maior destaque foi anos_estudo:grupo_idade. Optamos por inserí-la pois consideramos que fazia sentido teoricamente: o efeito da escolaridade sobre a participação na FT varia conforme a idade (mais anos de estudos aumenta a chance de entrada na FT para jovens e de saída para idosos) - e virse versa.

A equação do modelo logístico estimado é:

\[\log\left(\frac{P(Y=1)}{1-P(Y=1)}\right) = \beta_0 + \beta_1 \text{recebeu_bf} + \beta_2 \text{sexo} + \beta_3 \text{grupo_idade} + \beta_4 \text{raca}\]

\[+ \beta_5 \text{anos_estudo} + \beta_6 \text{cond_dom} + \beta_7 \text{grupo_renda_dom_pc} + \beta_8 \text{regiao}\]

\[+ \beta_9 \text{area_urbana} + \beta_{10} (\text{recebeu_bf} \times \text{grupo_idade}) + \beta_{11} (\text{recebeu_bf} \times \text{anos_estudo})\]

\[+ \beta_{12} (\text{recebeu_bf} \times \text{cond_dom}) + \beta_{13} (\text{recebeu_bf} \times \text{grupo_renda_dom_pc})\]

\[+ \beta_{14} (\text{recebeu_bf} \times \text{area_urbana}) + \beta_{15} (\text{anos_estudo} \times \text{grupo_idade}) + \varepsilon\]

onde \(Y = 1\) indica não participação na força de trabalho, \(P(Y=1)\) é a probabilidade de não participar da força de trabalho, os \(\beta\) são os coeficientes estimados para cada variável e suas interações, e \(\varepsilon\) é o termo de erro.

# 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 + recebeu_bf:grupo_idade +
        recebeu_bf:anos_estudo + recebeu_bf:cond_dom + recebeu_bf:grupo_renda_dom_pc +
        recebeu_bf:area_urbana + anos_estudo:grupo_idade, 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")

modelo_tidy <- tidy(glm_base_total_final, conf.int = TRUE, exponentiate = TRUE) %>%
    mutate(OR = round(estimate, 3), IC_inf = round(conf.low, 3), IC_sup = round(conf.high,
        3), IC_95 = paste0("(", IC_inf, " - ", IC_sup, ")"), p_valor = case_when(p.value <
        0.001 ~ "<0.001", p.value < 0.01 ~ "<0.01", p.value < 0.05 ~ "<0.05", TRUE ~
        as.character(round(p.value, 3))), Significancia = case_when(p.value < 0.001 ~
        "***", p.value < 0.01 ~ "**", p.value < 0.05 ~ "*", TRUE ~ "")) %>%
    select(term, OR, IC_95, p_valor, Significancia)

# Criar tabela
modelo_tidy %>%
    kbl(caption = "Saída do Modelo", col.names = c("Variável", "OR", "IC 95%",
        "p-valor", "Sig."), booktabs = TRUE, longtable = TRUE, format = "html") %>%
    kable_styling(latex_options = c("striped", "scale_down", "hold_position", "repeat_header"),
        font_size = 12) %>%
    footnote(general = "Odds Ratios com erros padrão robustos (HC1)", general_title = "Nota:",
        footnote_as_chunk = TRUE)
Saída do Modelo
Variável OR IC 95% p-valor Sig.
(Intercept) 0.282 (0.251 - 0.317) <0.001 ***
recebeu_bf1 1.314 (1.063 - 1.625) <0.05
sexo2 2.373 (2.267 - 2.484) <0.001 ***
grupo_idade25 a 39 0.375 (0.343 - 0.409) <0.001 ***
grupo_idade40 a 64 0.413 (0.379 - 0.45) <0.001 ***
grupo_idade65 ou mais 1.207 (1.064 - 1.367) <0.01 **
racaPP 1.077 (1.026 - 1.13) <0.01 **
racaO 1.630 (1.372 - 1.928) <0.001 ***
anos_estudo12 a 15 0.397 (0.367 - 0.429) <0.001 ***
anos_estudo16 ou mais 0.323 (0.23 - 0.443) <0.001 ***
cond_domFilho/Enteado/Genro/Nora/Neto/Bisneto 2.241 (2.101 - 2.391) <0.001 ***
cond_domPai/Mãe/Padrasto/Madrasta/Sogro/Sogra/Avô/Avó 0.614 (0.481 - 0.772) <0.001 ***
grupo_renda_dom_pcR$ 218,01 a 1,00 SM 0.232 (0.216 - 0.249) <0.001 ***
grupo_renda_dom_pc1,01 SM a 2,00 SM 0.093 (0.085 - 0.101) <0.001 ***
regiaoNE 2.138 (2.019 - 2.265) <0.001 ***
regiaoN 1.274 (1.184 - 1.369) <0.001 ***
regiaoCO 1.100 (1.004 - 1.204) <0.05
regiaoS 0.747 (0.682 - 0.816) <0.001 ***
area_urbana2 1.440 (1.37 - 1.513) <0.001 ***
recebeu_bf1:grupo_idade25 a 39 1.463 (1.248 - 1.717) <0.001 ***
recebeu_bf1:grupo_idade40 a 64 1.155 (0.981 - 1.362) 0.085
recebeu_bf1:grupo_idade65 ou mais 0.255 (0.138 - 0.441) <0.001 ***
recebeu_bf1:anos_estudo12 a 15 1.330 (1.198 - 1.477) <0.001 ***
recebeu_bf1:anos_estudo16 ou mais 1.359 (1.002 - 1.822) <0.05
recebeu_bf1:cond_domFilho/Enteado/Genro/Nora/Neto/Bisneto 0.681 (0.579 - 0.8) <0.001 ***
recebeu_bf1:cond_domPai/Mãe/Padrasto/Madrasta/Sogro/Sogra/Avô/Avó 0.949 (0.533 - 1.62) 0.853
recebeu_bf1:grupo_renda_dom_pcR$ 218,01 a 1,00 SM 1.275 (1.098 - 1.48) <0.01 **
recebeu_bf1:grupo_renda_dom_pc1,01 SM a 2,00 SM 1.001 (0.772 - 1.29) 0.992
recebeu_bf1:area_urbana2 1.205 (1.098 - 1.321) <0.001 ***
grupo_idade25 a 39:anos_estudo12 a 15 1.554 (1.389 - 1.74) <0.001 ***
grupo_idade40 a 64:anos_estudo12 a 15 1.625 (1.454 - 1.815) <0.001 ***
grupo_idade65 ou mais:anos_estudo12 a 15 2.858 (2.235 - 3.629) <0.001 ***
grupo_idade25 a 39:anos_estudo16 ou mais 1.419 (1 - 2.053) 0.056
grupo_idade40 a 64:anos_estudo16 ou mais 1.245 (0.871 - 1.811) 0.24
grupo_idade65 ou mais:anos_estudo16 ou mais 3.778 (2.208 - 6.387) <0.001 ***
Nota: Odds Ratios com erros padrão robustos (HC1)

Testagens sobre o Modelo

  • Análise do Ajuste do Modelo: O modelo apresenta qualidade de ajuste razoável, com Pseudo-R² de McFadden de 0,18 e de Nagelkerke de 0,23. O AIC e BIC sugerem um modelo razoável. O PCP (Percentage Correctly Predicted) de 86,7% demonstram capacidade preditiva adequada.

  • Heterocedasticidade: Todos os testes aplicados (Breusch-Pagan e White) rejeitam a hipótese nula de homocedasticidade com p-valores < 0,001, confirmando a presença de heterocedasticidade nos resíduos. Este problema foi tratado através do uso de erros padrão robustos (HC1) em todas as estimações e testes subsequentes, garantindo a validade das inferências estatísticas.

  • Multicolinearidade: A análise dos fatores de inflação da variância (GVIFs) revela problemas relevantes de multicolinearidade. Os GVIFs ajustados para variáveis categóricas mostram valores aceitáveis (< 2,2) para a maioria das variáveis, mas a multicolinearidade entre educação, idade e suas interações deve ser considerada na interpretação dos resultados.

  • Significância Individual: Os testes de Razão de Verossimilhança (tipo III) confirma que todas as variáveis principais e suas interações são estatisticamente significativas (p < 0,001). A variável principal recebeu_bf mantém significância mesmo controlando por todas as demais variáveis e interações, o que implicou em alta multicolinearidade.

  • Significância Global: O teste de Razão de Verossimilhança global (p < 0,001) garante a significância geral do modelo.

# 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")

    # 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, type = "predictor")
    print(vif_values)
    cat("\nInterpretação VIF:\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")
    cat("\nInterpretação GVIF (mais adequado para variáveis categóricas):\n")
    cat("- GVIF < 2.2: Multicolinearidade aceitável\n")
    cat("- GVIF entre 2.2-3.2: Multicolinearidade moderada\n")
    cat("- GVIF > 3.2: Multicolinearidade severa\n\n")

    # 4. SIGNIFICÂNCIA INDIVIDUAL
    cat("4. SIGNIFICÂNCIA INDIVIDUAL\n")
    cat("---------------------------\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.1Teste de Razão de Verossimilhança (LRT) global\n")
    cat("---------------------------------------\n")
    tryCatch({
        # Teste LRT global: modelo completo vs. modelo só com intercepto
        model_null <- glm(na_forca_trabalho ~ 1, family = binomial, data = model$data)

        lrt_test <- anova(model_null, model, test = "LRT")
        print(lrt_test)

        # Interpretação
        p_val_lrt <- lrt_test$`Pr(>Chi)`[2]
        cat("\nInterpretação LRT:")
        if (p_val_lrt < 0.001) {
            cat(" Modelo altamente significativo (p < 0.001) ***\n")
        } else if (p_val_lrt < 0.01) {
            cat(" Modelo muito significativo (p < 0.01) **\n")
        } else if (p_val_lrt < 0.05) {
            cat(" Modelo significativo (p < 0.05) *\n")
        } else {
            cat(" Modelo não significativo (p >= 0.05)\n")
        }
        cat("\n")
    }, error = function(e) {
        cat("Erro no teste LRT:", e$message, "\n\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.1847 
## Cox & Snell R²: 0.1015 
## Nagelkerke (Cragg-Uhler) R²: 0.2308 
## 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
## ------------------------------------------------------------------------
## 69915.660 | 69915.677 | 70262.275 |     0.139 | 0.259 | 1.000 |    0.236
## 
## AIC       | Score_log | Score_spherical |   PCP
## -----------------------------------------------
## 69915.660 |      -Inf |       1.226e-05 | 0.867
## 
## 2. HETEROCEDASTICIDADE
## ----------------------
## Teste de Breusch-Pagan:
## 
##  studentized Breusch-Pagan test
## 
## data:  model
## BP = 16299, df = 34, p-value < 2.2e-16
## 
## 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: 170.35
## 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                     22.295381  1        4.721798
## sexo                            1.290222  1        1.135880
## grupo_idade                    12.226134  3        1.517801
## raca                            1.113980  2        1.027352
## anos_estudo                    39.617119  2        2.508827
## cond_dom                        2.999644  2        1.316035
## grupo_renda_dom_pc              1.711732  2        1.143823
## regiao                          1.203771  4        1.023453
## area_urbana                     1.525441  1        1.235087
## recebeu_bf:grupo_idade         19.453301  3        1.639956
## recebeu_bf:anos_estudo          2.825257  2        1.296476
## recebeu_bf:cond_dom             1.636783  2        1.131092
## recebeu_bf:grupo_renda_dom_pc  11.242746  2        1.831125
## recebeu_bf:area_urbana          2.416817  1        1.554611
## grupo_idade:anos_estudo       118.055447  6        1.488243
## 
## Interpretação VIF:
## - VIF < 5: Não há problema de multicolinearidade
## - VIF entre 5-10: Multicolinearidade moderada
## - VIF > 10: Multicolinearidade severa
## 
## 
## Interpretação GVIF (mais adequado para variáveis categóricas):
## - GVIF < 2.2: Multicolinearidade aceitável
## - GVIF entre 2.2-3.2: Multicolinearidade moderada
## - GVIF > 3.2: Multicolinearidade severa
## 
## 4. SIGNIFICÂNCIA INDIVIDUAL
## ---------------------------
## 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                        6.34  1  0.0118235 *  
## sexo                           1400.02  1  < 2.2e-16 ***
## grupo_idade                     834.39  3  < 2.2e-16 ***
## raca                             33.06  2  6.631e-08 ***
## anos_estudo                     565.52  2  < 2.2e-16 ***
## cond_dom                        619.45  2  < 2.2e-16 ***
## grupo_renda_dom_pc             2640.74  2  < 2.2e-16 ***
## regiao                         1267.50  4  < 2.2e-16 ***
## area_urbana                     203.61  1  < 2.2e-16 ***
## recebeu_bf:grupo_idade           73.55  3  7.396e-16 ***
## recebeu_bf:anos_estudo           29.92  2  3.186e-07 ***
## recebeu_bf:cond_dom              22.13  2  1.564e-05 ***
## recebeu_bf:grupo_renda_dom_pc    14.11  2  0.0008648 ***
## recebeu_bf:area_urbana           15.59  1  7.859e-05 ***
## grupo_idade:anos_estudo         142.75  6  < 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.1Teste de Razão de Verossimilhança (LRT) global
## ---------------------------------------
## Analysis of Deviance Table
## 
## Model 1: na_forca_trabalho ~ 1
## Model 2: na_forca_trabalho ~ recebeu_bf + sexo + grupo_idade + raca + 
##     anos_estudo + cond_dom + grupo_renda_dom_pc + regiao + area_urbana + 
##     recebeu_bf:grupo_idade + recebeu_bf:anos_estudo + recebeu_bf:cond_dom + 
##     recebeu_bf:grupo_renda_dom_pc + recebeu_bf:area_urbana + 
##     anos_estudo:grupo_idade
##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
## 1    147751      85666                          
## 2    147717      69846 34    15820 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Interpretação LRT: Modelo altamente significativo (p < 0.001) ***
## 
## ===============================================
## FIM DOS DIAGNÓSTICOS
## ===============================================

Resultados: APEs

  • Como as modelagens realizadas possuem muitas variáveis categóricas de controle, a interpretação do coeficiente de recebeu_bf e de suas interações como saem no modelo seriam inapropriada: a comparação seria sempre em relação ao grupo de referência. Além disso, o número de dummys é muito alto. Por isso, usamos o pacote marginaleffects para calcular o efeito parcial médio (APE) geral de recebeu_bf e seus efeitos agrupados pelas variáveis de controle. Isso foi feito usando a função comparisons(), indicada pelo pacote como melhor opção para variáveis dependentes binárias. A significância estatística dos APEs também foi testada e seus resultados são reportados;

Efeito Parcial Médio Geral

  • Efeito parcial médio: O recebimento do Bolsa Família está associado a um aumento de 4,81 pontos percentuais na probabilidade de não participar da força de trabalho.

  • Interpretação do resultado: Este achado argumenta no sentido da validade da ideia de que a ampliação da renda não laboral implica em redução oferta de trabalho.

  • Interpretação da significância estatística: O resultado é altamente significativo (p < 0,001) com intervalo de confiança entre 4,19% e 5,42%.

# Função para calcular efeito médio geral (versão final)
analyze_average_marginal_effect <- function(model, variable, vcov_matrix, outcome_label) {

    # Calcular efeito marginal médio geral
    avg_effect <- avg_comparisons(model, variables = variable, type = "response",
        vcov = vcov_matrix)

    # Gráfico do efeito marginal médio
    plot_effect <- avg_effect %>%
        mutate(Significativo = ifelse(p.value < 0.05, "Sim", "Não")) %>%
        ggplot(aes(x = "Efeito do BF", y = estimate, color = Significativo)) + geom_hline(yintercept = 0,
        linetype = "dashed", color = "gray50") + geom_point(size = 5) + geom_errorbar(aes(ymin = conf.low,
        ymax = conf.high), width = 0.1, size = 1) + geom_text(aes(label = paste0(scales::percent(estimate,
        accuracy = 0.01), "\n(p = ", round(p.value, 3), ")")), vjust = -1.2, size = 3.5,
        fontface = "bold") + labs(title = "Efeito Médio do Bolsa Família", subtitle = paste("Mudança na probabilidade de",
        outcome_label, "[Recebeu vs Não Recebeu]"), x = "", y = "Efeito Marginal Médio",
        color = "Significativo (p<0.05)") + scale_color_manual(values = c(Sim = "#E41A1C",
        Não = "#377EB8")) + scale_y_continuous(labels = scales::percent_format()) +
        theme_minimal(base_size = 12) + theme(legend.position = "top", plot.subtitle = element_text(size = 10,
        color = "gray60"))

    return(list(avg_effect = avg_effect, plot_effects = plot_effect))
}

# ============================================ EXECUÇÃO DA ANÁLISE
# ============================================

# Calcular efeito médio geral
resultado_medio <- analyze_average_marginal_effect(model = glm_base_total_final,
    variable = "recebeu_bf", vcov_matrix = vcov_rob_total_3, outcome_label = "Participação na Força de Trabalho")

# Mostrar resultados
cat("=== EFEITO MÉDIO GERAL DO BOLSA FAMÍLIA ===\n\n")
## === EFEITO MÉDIO GERAL DO BOLSA FAMÍLIA ===
cat("--- EFEITO MARGINAL MÉDIO ---\n")
## --- EFEITO MARGINAL MÉDIO ---
print(resultado_medio$avg_effect)
## 
##  Estimate Std. Error    z Pr(>|z|)     S  2.5 % 97.5 %
##    0.0481    0.00313 15.3   <0.001 174.1 0.0419 0.0542
## 
## Term: recebeu_bf
## Type: response
## Comparison: 1 - 0
cat("\n")
# Extrair valores para resumo
efeito <- resultado_medio$avg_effect$estimate[1]
p_valor <- resultado_medio$avg_effect$p.value[1]
ic_low <- resultado_medio$avg_effect$conf.low[1]
ic_high <- resultado_medio$avg_effect$conf.high[1]
erro_padrao <- resultado_medio$avg_effect$std.error[1]
estatistica_z <- resultado_medio$avg_effect$statistic[1]

cat("--- RESUMO DO EFEITO ---\n")
## --- RESUMO DO EFEITO ---
cat(sprintf("Efeito marginal médio: %.4f (%.2f p.p.)\n", efeito, efeito * 100))
## Efeito marginal médio: 0.0481 (4.81 p.p.)
cat(sprintf("Erro padrão: %.4f\n", erro_padrao))
## Erro padrão: 0.0031
cat(sprintf("Estatística Z: %.3f\n", estatistica_z))
## Estatística Z: 15.344
cat(sprintf("P-valor: %.6f\n", p_valor))
## P-valor: 0.000000
cat(sprintf("Intervalo de confiança (95%%): [%.4f, %.4f]\n", ic_low, ic_high))
## Intervalo de confiança (95%): [0.0419, 0.0542]
cat(sprintf("Significativo a 5%%: %s\n", ifelse(p_valor < 0.05, "SIM", "NÃO")))
## Significativo a 5%: SIM
cat("\n")
# Mostrar gráficos
resultado_medio$plot_effects

Efeitos Parciais Médios por grupo

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")

    # 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, plot_effects = plot_effects))
}

Sexo

  • Efeito parcial médio: Para homens, o programa está associado a um aumento de 3,52 pontos percentuais na probabilidade de não participar da força de trabalho. Para mulheres, este efeito é substancialmente maior, chegando a 6,41 pontos percentuais.

  • Interpretação do resultado: A heterogeneidade por gênero revela que as mulheres são mais responsivas ao efeito renda do Bolsa Família do que os homens. Isso pode refletir diferenças nos papéis de gênero tradicionais, onde mulheres frequentemente assumem responsabilidades domésticas e de cuidado que tornam a saída da força de trabalho mais viável quando há renda alternativa disponível.

  • Interpretação da significância estatística: Ambos os efeitos são estatisticamente significativos (p < 0,001), e o teste de diferença entre os grupos confirma que o efeito é significativamente maior para mulheres (diferença de 2,89 pontos percentuais, p < 0,001).

cat("# EFEITOS MARGINAIS - MARGEM EXTENSIVA\n")
## # EFEITOS MARGINAIS - MARGEM EXTENSIVA
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")
## 
## === EFEITO DO BOLSA FAMÍLIA POR SEXO ===
## 
## --- EFEITOS MARGINAIS ---
## 
##  sexo Estimate Std. Error    z Pr(>|z|)     S  2.5 % 97.5 %
##     1   0.0352    0.00277 12.7   <0.001 120.6 0.0298 0.0407
##     2   0.0641    0.00379 16.9   <0.001 210.8 0.0567 0.0716
## 
## Term: recebeu_bf
## Type: response
## Comparison: 1 - 0

## Teste Wald (robusto) de diferenças entre sexo:
## 
##  Hypothesis Estimate Std. Error     z Pr(>|z|)     S   2.5 %  97.5 %
##       b1=b2  -0.0289    0.00184 -15.7   <0.001 181.8 -0.0325 -0.0253
##       b1=0    0.0352    0.00277  12.7   <0.001 120.6  0.0298  0.0407
##       b2=0    0.0641    0.00379  16.9   <0.001 210.8  0.0567  0.0716

Idade

  • Efeito parcial médio: Para jovens até 24 anos, o programa está associado a um aumento de 4,35 pontos percentuais na probabilidade de não participar da força de trabalho. Para adultos de 25 a 39 anos, o efeito é maior (6,69 pontos percentuais), assim como para aqueles de 40 a 64 anos (4,51 pontos percentuais). Para idosos de 65 anos ou mais, o efeito é negativo (-5,33 pontos percentuais), indicando que o programa está associado a maior participação laboral neste grupo.

  • Interpretação do resultado: A heterogeneidade etária sugere diferentes mecanismos de resposta ao programa. Adultos em idade produtiva plena (25-39 anos) mostram maior sensibilidade ao efeito renda, possivelmente porque têm maior flexibilidade para escolher entre trabalho e outras atividades. O efeito negativo entre idosos muito provavelmente reflete diferenças de seletividade: idosos beneficiários são predominantemente de baixa renda e continuam trabalhando por necessidade econômica, enquanto idosos não-beneficiários incluem maior proporção de pessoas com renda suficiente (aposentadorias, pensões) que permite a saída da força de trabalho. O programa, neste caso, funciona como complemento insuficiente para permitir a inatividade.

  • Interpretação da significância estatística: Todos os efeitos são estatisticamente significativos (p < 0,001), e os testes de diferença entre grupos confirmam heterogeneidade significativa, especialmente entre idosos e os demais grupos.

# 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 Faixa Etária",
    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.0435    0.01169  3.72   <0.001  12.3  0.0206  0.0664
##   25 a 39      0.0669    0.00373 17.95   <0.001 236.9  0.0596  0.0742
##   40 a 64      0.0451    0.00289 15.59   <0.001 179.6  0.0394  0.0507
##   65 ou mais  -0.0533    0.01397 -3.82   <0.001  12.9 -0.0807 -0.0260
## 
## Term: recebeu_bf
## Type: response
## Comparison: 1 - 0

## Teste Wald (robusto) de diferenças entre grupo_idade:
## 
##  Hypothesis Estimate Std. Error      z Pr(>|z|)     S   2.5 %   97.5 %
##       b1=b2 -0.02341    0.01127 -2.076   0.0379   4.7 -0.0455 -0.00131
##       b1=b3 -0.00156    0.01170 -0.133   0.8942   0.2 -0.0245  0.02137
##       b1=b4  0.09685    0.01804  5.370   <0.001  23.6  0.0615  0.13220
##       b2=b3  0.02185    0.00416  5.258   <0.001  22.7  0.0137  0.03000
##       b2=b4  0.12025    0.01430  8.408   <0.001  54.4  0.0922  0.14828
##       b3=b4  0.09840    0.01411  6.972   <0.001  38.2  0.0707  0.12606
##       b1=0   0.04351    0.01169  3.722   <0.001  12.3  0.0206  0.06643
##       b2=0   0.06692    0.00373 17.951   <0.001 236.9  0.0596  0.07422
##       b3=0   0.04507    0.00289 15.589   <0.001 179.6  0.0394  0.05073
##       b4=0  -0.05334    0.01397 -3.818   <0.001  12.9 -0.0807 -0.02595

Raça

  • Efeito parcial médio: Para brancos, o programa está associado a um aumento de 3,51 pontos percentuais na probabilidade de não participar da força de trabalho. Para pretos e pardos, o efeito é maior (5,45 pontos percentuais), e para outros grupos raciais, o efeito é ainda mais pronunciado (6,85 pontos percentuais).

  • Interpretação do resultado: A heterogeneidade racial sugere que grupos historicamente marginalizados são mais responsivos ao efeito renda do Bolsa Família. Isso pode refletir diferenças nas oportunidades de trabalho disponíveis, qualidade dos empregos acessíveis, e níveis de renda prévia. Pretos, pardos e outros grupos raciais podem enfrentar maior discriminação no mercado de trabalho e ter acesso predominante a empregos de menor qualidade e remuneração, tornando a transferência de renda uma alternativa mais atrativa.

  • Interpretação da significância estatística: Todos os efeitos são estatisticamente significativos (p < 0,001), e os testes de diferença entre grupos confirmam que o efeito é significativamente maior para pretos/pardos e outros grupos em relação aos brancos.

# 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.0351    0.00277 12.7   <0.001 119.9 0.0297 0.0405
##    PP   0.0545    0.00338 16.1   <0.001 192.0 0.0478 0.0611
##    O    0.0685    0.00526 13.0   <0.001 126.5 0.0582 0.0788
## 
## Term: recebeu_bf
## Type: response
## Comparison: 1 - 0

## Teste Wald (robusto) de diferenças entre raca:
## 
##  Hypothesis Estimate Std. Error      z Pr(>|z|)     S   2.5 %  97.5 %
##       b1=b2  -0.0194    0.00124 -15.59   <0.001 179.6 -0.0218 -0.0169
##       b1=b3  -0.0334    0.00390  -8.56   <0.001  56.3 -0.0410 -0.0257
##       b2=b3  -0.0140    0.00374  -3.75   <0.001  12.5 -0.0214 -0.0067
##       b1=0    0.0351    0.00277  12.67   <0.001 119.9  0.0297  0.0405
##       b2=0    0.0545    0.00338  16.13   <0.001 192.0  0.0478  0.0611
##       b3=0    0.0685    0.00526  13.03   <0.001 126.5  0.0582  0.0788

Anos de estudo

  • Efeito parcial médio: Para pessoas com até 11 anos de estudo, o programa está associado a um aumento de 4,76 pontos percentuais na probabilidade de não participar da força de trabalho. Para aqueles com 12 a 15 anos de estudo, o efeito é ligeiramente maior (5,22 pontos percentuais). Para pessoas com 16 anos ou mais de estudo, o efeito é menor (3,23 pontos percentuais).

  • Interpretação do resultado: A heterogeneidade educacional revela que pessoas com menor escolaridade são mais responsivas ao efeito renda do Bolsa Família. Isso pode refletir que indivíduos com menor educação formal têm acesso predominante a empregos de baixa remuneração, condições precárias e menor estabilidade, tornando a transferência de renda uma alternativa mais atrativa. Para pessoas com ensino superior, que geralmente têm acesso a melhores oportunidades de trabalho e maior potencial de renda, o efeito do programa é reduzido.

  • Interpretação da significância estatística: Todos os efeitos são estatisticamente significativos (p < 0,001). Os testes de diferença entre grupos mostram que apenas a diferença entre os grupos de 12-15 anos e 16+ anos é significativa (p = 0,0157), enquanto as demais comparações não apresentam diferenças estatisticamente significativas.

# 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.0476    0.00396 12.01   <0.001 108.0 0.0398 0.0553
##   12 a 15      0.0522    0.00386 13.51   <0.001 135.8 0.0446 0.0598
##   16 ou mais   0.0323    0.00781  4.14   <0.001  14.8 0.0170 0.0476
## 
## Term: recebeu_bf
## Type: response
## Comparison: 1 - 0

## Teste Wald (robusto) de diferenças entre anos_estudo:
## 
##  Hypothesis Estimate Std. Error     z Pr(>|z|)     S    2.5 %  97.5 %
##       b1=b2 -0.00465    0.00438 -1.06   0.2881   1.8 -0.01323 0.00393
##       b1=b3  0.01524    0.00835  1.83   0.0677   3.9 -0.00111 0.03160
##       b2=b3  0.01989    0.00823  2.42   0.0157   6.0  0.00376 0.03603
##       b1=0   0.04755    0.00396 12.01   <0.001 108.0  0.03979 0.05531
##       b2=0   0.05220    0.00386 13.51   <0.001 135.8  0.04463 0.05977
##       b3=0   0.03231    0.00781  4.14   <0.001  14.8  0.01700 0.04762

Condição no domicílio

  • Efeito parcial médio: Para responsáveis/cônjuges, o programa está associado a um aumento de 4,98 pontos percentuais na probabilidade de não participar da força de trabalho. Para filhos/enteados/genros/noras/netos/bisnetos, o efeito é ligeiramente menor (4,47 pontos percentuais). Para pais/mães/padrastos/madrastas/sogros/sogras/avôs/avós, o efeito é substancialmente menor e não significativo (1,92 pontos percentuais, p = 0,16).

  • Interpretação do resultado: A heterogeneidade por condição domiciliar sugere que responsáveis pelo domicílio são mais sensíveis ao efeito renda do programa, assim como filhos/netos. Para outros parentes (pais/sogros/avós), o efeito é estatisticamente não significativo, sugerindo que sua participação laboral é menos influenciada pelo recebimento do benefício.

  • Interpretação da significância estatística: Os efeitos para responsáveis/cônjuges e filhos/dependentes são estatisticamente significativos (p < 0,001), enquanto o efeito para outros parentes não é significativo (p = 0,16). Os testes de diferença entre grupos confirmam heterogeneidade significativa.

# 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.0498    0.00247 20.17
##  Filho/Enteado/Genro/Nora/Neto/Bisneto           0.0447    0.00950  4.71
##  Pai/Mãe/Padrasto/Madrasta/Sogro/Sogra/Avô/Avó   0.0192    0.01369  1.40
##  Pr(>|z|)     S   2.5 % 97.5 %
##    <0.001 298.1  0.0450 0.0546
##    <0.001  18.6  0.0261 0.0634
##      0.16   2.6 -0.0076 0.0461
## 
## Term: recebeu_bf
## Type: response
## Comparison: 1 - 0

## Teste Wald (robusto) de diferenças entre cond_dom:
## 
##  Hypothesis Estimate Std. Error      z Pr(>|z|)     S    2.5 % 97.5 %
##       b1=b2  0.00506    0.00953  0.531   0.5954   0.7 -0.01362 0.0237
##       b1=b3  0.03057    0.01374  2.224   0.0262   5.3  0.00363 0.0575
##       b2=b3  0.02550    0.01658  1.539   0.1239   3.0 -0.00698 0.0580
##       b1=0   0.04981    0.00247 20.168   <0.001 298.1  0.04497 0.0546
##       b2=0   0.04474    0.00950  4.711   <0.001  18.6  0.02613 0.0634
##       b3=0   0.01924    0.01369  1.405   0.1600   2.6 -0.00760 0.0461

Renda Dom. P.C.

  • Efeito parcial médio: Para famílias com renda até R$ 218,00, o programa está associado a um aumento de 10,67 pontos percentuais na probabilidade de não participar da força de trabalho. Para famílias com renda entre R$ 218,01 e 1,00 salário mínimo, o efeito é menor (6,89 pontos percentuais). Para famílias com renda entre 1,01 e 2,00 salários mínimos, o efeito é substancialmente reduzido (1,12 pontos percentuais).

  • Interpretação do resultado: A heterogeneidade por renda revela um padrão claro: quanto menor a renda domiciliar, maior o efeito do Bolsa Família na redução da participação laboral. Para famílias em extrema pobreza (até R$ 218,00), o programa representa uma fonte de renda proporcionalmente muito significativa, criando forte incentivo para reduzir a oferta de trabalho. À medida que a renda domiciliar aumenta, o valor do benefício torna-se proporcionalmente menor em relação à renda total, reduzindo o efeito sobre as decisões laborais.

  • Interpretação da significância estatística: Todos os efeitos são estatisticamente significativos, sendo altamente significativos para as duas primeiras faixas (p < 0,001) e significativo para a terceira faixa (p = 0,0037). Os testes de diferença entre grupos confirmam heterogeneidade significativa entre todas as comparações.

# 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 Faixa de Renda Domiciliar 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.1067    0.01548  6.89   <0.001  37.4 0.07631 0.1370
##  R$ 218,01 a 1,00 SM   0.0689    0.00403 17.09   <0.001 215.1 0.06099 0.0768
##  1,01 SM a 2,00 SM     0.0112    0.00387  2.90   0.0037   8.1 0.00365 0.0188
## 
## Term: recebeu_bf
## Type: response
## Comparison: 1 - 0

## 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.0378    0.01510  2.50   0.0124   6.3 0.00817 0.0674
##       b1=b3   0.0954    0.01571  6.07   <0.001  29.6 0.06462 0.1262
##       b2=b3   0.0576    0.00514 11.21   <0.001  94.4 0.04756 0.0677
##       b1=0    0.1067    0.01548  6.89   <0.001  37.4 0.07631 0.1370
##       b2=0    0.0689    0.00403 17.09   <0.001 215.1 0.06099 0.0768
##       b3=0    0.0112    0.00387  2.90   0.0037   8.1 0.00365 0.0188

Região

  • Efeito parcial médio: No Sudeste, o programa está associado a um aumento de 2,94 pontos percentuais na probabilidade de não participar da força de trabalho. No Nordeste, o efeito é substancialmente maior (7,78 pontos percentuais). No Norte, o efeito é de 4,93 pontos percentuais. No Centro-Oeste, é de 3,08 pontos percentuais. No Sul, o efeito é o menor (2,19 pontos percentuais).

  • Interpretação do resultado: A heterogeneidade regional revela que o Nordeste apresenta a maior responsividade ao programa, seguido pelo Norte. Isso pode refletir diferenças nas estruturas econômicas regionais, oportunidades de trabalho disponíveis, e níveis de desenvolvimento. Regiões com mercados de trabalho menos desenvolvidos, menor diversificação econômica e predominância de empregos de baixa remuneração (como Nordeste e Norte) podem tornar o Bolsa Família uma alternativa mais atrativa. Regiões mais desenvolvidas economicamente (Sudeste e Sul) apresentam menores efeitos, possivelmente devido a melhores oportunidades de trabalho e maior competitividade salarial.

  • Interpretação da significância estatística: Todos os efeitos regionais são estatisticamente significativos (p < 0,001). Os testes de diferença entre grupos confirmam heterogeneidade significativa, especialmente entre o Nordeste e as demais regiões.

# 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.0294    0.00252 11.7   <0.001 102.0 0.0244 0.0343
##      NE   0.0778    0.00448 17.4   <0.001 221.6 0.0690 0.0866
##      N    0.0493    0.00334 14.8   <0.001 161.4 0.0427 0.0558
##      CO   0.0308    0.00284 10.8   <0.001  88.6 0.0253 0.0364
##      S    0.0219    0.00210 10.4   <0.001  82.4 0.0178 0.0261
## 
## Term: recebeu_bf
## Type: response
## Comparison: 1 - 0

## Teste Wald (robusto) de diferenças entre regiao:
## 
##  Hypothesis Estimate Std. Error      z Pr(>|z|)     S    2.5 %    97.5 %
##       b1=b2 -0.04844    0.00259 -18.71   <0.001 257.0 -0.05351 -0.043361
##       b1=b3 -0.01990    0.00173 -11.48   <0.001  99.0 -0.02330 -0.016504
##       b1=b4 -0.00147    0.00120  -1.22    0.222   2.2 -0.00382  0.000887
##       b1=b5  0.00743    0.00113   6.58   <0.001  34.3  0.00521  0.009641
##       b2=b3  0.02854    0.00175  16.34   <0.001 196.9  0.02511  0.031960
##       b2=b4  0.04697    0.00269  17.45   <0.001 224.0  0.04169  0.052244
##       b2=b5  0.05586    0.00304  18.37   <0.001 247.9  0.04990  0.061826
##       b3=b4  0.01843    0.00193   9.53   <0.001  69.0  0.01464  0.022223
##       b3=b5  0.02733    0.00208  13.12   <0.001 128.1  0.02324  0.031411
##       b4=b5  0.00890    0.00145   6.13   <0.001  30.1  0.00605  0.011740
##       b1=0   0.02937    0.00252  11.67   <0.001 102.0  0.02444  0.034305
##       b2=0   0.07781    0.00448  17.35   <0.001 221.6  0.06902  0.086596
##       b3=0   0.04927    0.00334  14.76   <0.001 161.4  0.04273  0.055812
##       b4=0   0.03084    0.00284  10.84   <0.001  88.6  0.02527  0.036413
##       b5=0   0.02194    0.00210  10.44   <0.001  82.4  0.01782  0.026060

Urbano ou Rural

  • Efeito parcial médio: Em áreas urbanas, o programa está associado a um aumento de 3,59 pontos percentuais na probabilidade de não participar da força de trabalho. Em áreas rurais, o efeito é substancialmente maior (7,91 pontos percentuais).

  • Interpretação do resultado: A heterogeneidade entre áreas urbanas e rurais mostra que populações rurais são mais responsivas ao efeito renda do Bolsa Família. Isso pode refletir diferenças nas estruturas de mercado de trabalho: áreas rurais frequentemente oferecem menos oportunidades de emprego formal, menor diversificação ocupacional e predominância de atividades agrícolas sazonais ou de baixa remuneração. Em contraste, áreas urbanas geralmente apresentam mercados de trabalho mais diversificados, maior competição salarial e mais oportunidades de emprego formal, tornando o trabalho uma alternativa mais atrativa mesmo com o recebimento do benefício.

  • Interpretação da significância estatística: Ambos os efeitos são estatisticamente significativos (p < 0,001). O teste de diferença entre grupos confirma que o efeito é significativamente maior em áreas rurais (diferença de 4,32 pontos percentuais, p < 0,001).

# 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.0359    0.00311 11.5   <0.001  99.7 0.0298 0.0420
##            2   0.0791    0.00521 15.2   <0.001 170.4 0.0689 0.0893
## 
## Term: recebeu_bf
## Type: response
## Comparison: 1 - 0

## Teste Wald (robusto) de diferenças entre area_urbana:
## 
##  Hypothesis Estimate Std. Error     z Pr(>|z|)     S   2.5 %  97.5 %
##       b1=b2  -0.0432    0.00487 -8.88   <0.001  60.4 -0.0528 -0.0337
##       b1=0    0.0359    0.00311 11.52   <0.001  99.7  0.0298  0.0420
##       b2=0    0.0791    0.00521 15.17   <0.001 170.4  0.0689  0.0893

Limitações Metodológicas

  1. Seleção não aleatória: O recebimento do Bolsa Família não é aleatório, podendo haver viés de seleção.

  2. Simultaneidade: A participação ou não na força de trabalho também pode afetar o recebimento de Bolsa Família.

  3. 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.

  4. Desenho amostral: Os cálculos foram realizados sem a calibragem da amostra pelos pesos amostrais disponíveis nos microdados anuais da PNADc.

  5. 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).