Uma análise comparada entre Trabalhadores Formais e Trabalhadores por Conta-Própria com CNPJ

---

1 Objetivo Geral

• Apresentar um cenário comparativo a partir do salário-hora entre trabalhadores formais e trabalhadores por conta-própria que possuem CNPJ, a partir da pesquisa PNADc

2 Escolhas Metodológicas

• A base de dados utilizada foi a Pesquisa Nacional por Amostra de Domicílios Contínua - PNADc

• O recorte temporal foi o 3º trimestre/2024. Essa escolha ocorre tendo em vista as volatilidades do mercado de trabalho no início e fim de ano.

• A comparação de salário-hora não é direta. Para formalizados foram inseridos valores de 13º, férias, 1/3 de férias e FGTS. Para trabalhadores por conta-própria com PJ que recebem abaixo de R$ 81.000,00 foi retirada a contribuição de R$ 70,00 referente ao MEI. Para os demais PJs foram retirados 20% do valor do salário-mínimo. As fórmulas estão disponíveis na seção “Salário-Hora”

• Foram retirados outliers abaixo de 1% e acima de 99% a partir das variáveis de rendimento e de horas trabalhadas. Depois do cálculo de salário-hora foram retirados também outliers a partir dessa variável

• Foram retirados trabalhadores caracterizados como subocupados por insuficiência de horas.

• Tendo em vista que se trata de pesquisa amostral, quando foram analisados os valores totais (com peso pós-estratificado) qualquer valor total dos grupos (formal e conta-própria) que estivesse com coeficiente de variação acima de 25% foi retirado da tabela de resultados.

• A PNADc utiliza a Classificação de Ocupações para as Pesquisas Domiciliares – COD (https://ftp.ibge.gov.br/Censos/Censo_Demografico_2010/metodologia/anexos/anexo_7_ocupacao_cod.pdf) que tem a seguinte estrutura:

---

3 Desenho da Base

3.1 Variáveis de Interesse

• Ano
• Trimestre
• UF
• V4010 - Código da ocupação (cargo ou função)
  • É utilizada a Classificação de Ocupações para as Pesquisas Domiciliares – COD
• VD4002 - Ocupados/Desocupados
• VD4009 - Posição na ocupação e categoria do emprego
  • 1 - Empregado no setor privado com carteira de trabalho assinada
  • 9 - Conta-própria #VD4010 - Grupamentos de atividade principal do empreendimento do trabalho
• VD4016 - Rendimento mensal habitual do trabalho principal
• VD4010 - Grupamentos de atividade principal do empreendimento do trabalho principal
  • Ex.: Agricultura, Indústria, Construção etc
• V4039 - “Quantas horas … trabalhava normalmente, por semana, nesse trabalho principal?”
• V4019 - “Esse negócio/empresa era registrado no Cadastro Nacional da Pessoa Jurídica - CNPJ?”
• VD4004A - Subocupação por insuficiência de horas habitualmente trabalhadas em todos os trabalhos
• VD3004 - Escolaridade
• VD3005 - Anos de Estudo
• V2007 - Sexo
• V2010 - Cor/Raça
• V2009 - Idade
• VD4011 - Grupamento Ocupacional
• VD4012 - Contribuição para instituto de previdência em qualquer trabalho da semana de referência para pessoas de 14 anos ou mais de idade

3.2 Pré-Processamento

#Temas para Rmarkdown - https://rpubs.com/ranydc/rmarkdown_themes

library(tidyverse)
library(readxl)
library(DT)
library(PNADcIBGE)
library(survey)
library(srvyr)
library(scales)
library(jtools) #função summ - parte inferencial

Sys.setlocale("LC_ALL", "pt_BR.UTF-8") #problemas com encoding

Será feita uma consulta somente à base de dados de 2024 - 3º Trimestre, tendo em vista a limitação de capacidade de processamento.

O código abaixo fará o download dos dados e salvará em um objeto do tipo .rds para uma pasta local.

##########RODAR SOMENTE A PRIMEIRA VEZ###############

# variaveis_selecionadas <- c("Ano", "Trimestre", "UF","V4010",
#                             "VD4002","VD4009","VD4010","VD4016",
#                             "V4039","V4019", "VD4004A", "VD3004", "VD3005",
#                             "V2007", "V2010", "V2009", "VD4012", "VD4011")
# 
# pnadc <- get_pnadc(year=2024, quarter=3,
#                    vars=variaveis_selecionadas, design = T, labels = F)
# class(pnadc)

# #TRANSFORMAR PARA TIBBLE PARA USAR A GRAMÁTICA TIDYVERSE
# pnad_srvyr <- as_survey(pnadc)
# class(pnad_srvyr)

# #PRIMEIRO FILTRO - SOMENTE OCUPADOS
# ocupados <- pnad_srvyr |>
#   filter(VD4002 == 1) |>
#   mutate(V1028 = as.numeric(V1028))
# class(ocupados)

# rm(pnad_srvyr)
# rm(pnadc)

# #SALVAR PARA USAR OFFLINE
# saveRDS(ocupados, file = "c:/Users/raphael.rsl/Desktop/R/pjclt/ocupados3.rds")

3.3 Preparação da Base de Dados

Aqui será realizada a junção entre a PNADc pura e a tabela COD das profissões.

ocupados <- readRDS("c:/Users/raphael.rsl/Desktop/R/pjclt/ocupados3.rds")
#class(ocupados)

#LER A TABELA DE OCUPAÇÕES
cod <- read_xls("c:/Users/raphael.rsl/Desktop/R/pjclt/Estrutura_Ocupacao_COD.xls") |>
  rename(grande_grupo = `Grande Grupo`,
         subgrupo_principal = `Subgrupo principal`,
         subgrupo = Subgrupo,
         grupo_base = `Grupo de base`,
         denominacao = Denominação) |>
  mutate( #NÃO TRAZER NAs, ou seja, se não for NA traz "denominacao"
    denom_subgrupo_principal = if_else(!is.na(subgrupo_principal), denominacao, NA_character_),
    denom_subgrupo = if_else(!is.na(subgrupo), denominacao, NA_character_),
    denom_grupo_base = if_else(!is.na(grupo_base), denominacao, NA_character_)) |>
  fill(grande_grupo, subgrupo_principal, denom_subgrupo_principal, 
       subgrupo, denom_subgrupo, grupo_base, denom_grupo_base, .direction = "down")
       #propaga os valores para baixo

#RETIRANDO NAs e possíveis duplicados
cod_map <- cod |>
  filter(!is.na(grupo_base)) |>
  select(grupo_base, denom_grupo_base, subgrupo, denom_subgrupo,
         subgrupo_principal, denom_subgrupo_principal, grande_grupo) |>
  distinct()

#aplicando a base COD na base geral
ocupados2 <- ocupados |> 
  mutate(V4010 = as.character(V4010),
         subgrupo_cod = substr(V4010, 1, 3),
         # busca grupo_base e subgrupo e faz o match
         grupo_base = cod_map$denom_grupo_base[match(V4010, cod_map$grupo_base)],
         subgrupo = cod_map$denom_subgrupo[match(subgrupo_cod, cod_map$subgrupo)])

# Retirados indivíduos classificados como subocupados
ocupados2 <- ocupados2 |> 
  filter(is.na(VD4004A) | VD4004A != 1)
#É necessário manter NA por conta da estrutura do tbl_svy

#outliers - retira 1% superior e 1% inferior
ocupados2 <- ocupados2 |> 
  filter(VD4016 >= quantile(VD4016, 0.01, na.rm = TRUE),
         VD4016 <= quantile(VD4016, 0.99, na.rm = TRUE),
         V4039 >= quantile(V4039, 0.01, na.rm = TRUE),
         V4039 <= quantile(V4039, 0.99, na.rm = TRUE)) 

Gráficos Exploratórios

#VERIFICAR SE É O CASO DE RETIRAR MAIS OUTLIERS A PARTIR DOS BOXPLOTS

# Boxplot para Rendimento
ggplot(ocupados2, aes(x = "", y = VD4016)) +
  geom_boxplot(outlier.color = "red") +
  labs(title = "Boxplot de Rendimento", y = "Rendimento") +
  theme_minimal()

# Boxplot para Jornada
ggplot(ocupados2, aes(x = "", y = V4039)) +
  geom_boxplot(outlier.color = "red") +
  labs(title = "Boxplot de Jornada", y = "Horas") +
  theme_minimal()  

3.4 Salário-Hora

Criou-se uma variável auxiliar: salario_hora com a seguinte composição

Para formais, consideramos: - 12 salários mensais habituais; - 13º salário (1 salário adicional); - 1/3 adicional de férias; - 8% de FGTS; - Apenas 11 meses de trabalho efetivo (pois 1 mês é de férias).

A fórmula do salário-hora é:

\[ \text{Salário-hora}_{formal} = \frac{(12 + 1 + \tfrac{1}{3}) \cdot \text{Salário Mensal Habitual} \cdot (1 + \text{FGTS})} {\text{Horas Semanais Habitual} \cdot \tfrac{52}{12} \cdot 11} \]

Para trabalhadores autônomos registrados como MEI, consideramos:
- Rendimento mensal habitual declarado;
- Desconto fixo da guia DAS (imposto mensal do MEI);
- Jornada habitual semanal multiplicada por \(52/12\) (para obter as horas mensais).

A fórmula do salário-hora é:

\[ \text{Salário-hora}_{MEI} = \frac{\text{Salário Mensal Habitual} - \text{DAS}} {\text{Horas Semanais Habitual} \cdot \tfrac{52}{12}} \]

Autônomos acima do limite MEI

Para trabalhadores por conta própria com rendimento acima do limite MEI, consideramos:

• Contribuição proporcional: aplica-se 11% diretamente sobre o rendimento mensal, considerado o valor declarado como pro-labore.

\[ \text{Salário-hora}_{CP}^{prop} = \frac{\text{Salário Mensal Habitual} \cdot (1 - 11\%)} {\text{Horas Semanais Habitual} \cdot \tfrac{52}{12}} \]

Além disso, foi criada uma variável auxiliar chamada categoria_emprego para facilitar a filtragem de formal, cnpj e outros.

# Parâmetros
sal_min <- 1412          # salário-mínimo em 2024
mei_mensal <- 81000 / 12  # R$ 6750 por mês
mei_das <- 71.60          # DAS-MEI 2025
fgts <- 0.08         # 8% FGTS
semana_mes <- 52/12  # ~4.333
opcao_inss <- "NOT_fixed_min_wage_11pct" # regra usada para contribuição direta no salário minimo - não usada
#ou 20% segundo teste foi para 20%
inss_valor <- 0.11

# --- converter campos para numeric / proteger divisões ---
ocupados2 <- ocupados2  |> 
  mutate(VD4016 = as.numeric(as.character(VD4016)),  
         V4039  = as.numeric(as.character(V4039)),   
         VD4009 = as.character(VD4009))

# --- cálculo atualizado ---
ocupados3 <- ocupados2 |>
  mutate(
    # horas efetivamente trabalhadas no mês
    horas_mes = if_else(!is.na(V4039) & V4039 > 0, V4039 * semana_mes, NA_real_),

    # salário-hora
    salario_hora = case_when(
      # Formal com férias, 13º e FGTS (11 meses de trabalho no ano)
      VD4009 == "1" & !is.na(horas_mes) ~ {
        sal_anual <- VD4016 * (12 + 1 + 1/3)      # 12 salários + 13º + 1/3 férias
        horas_ano <- horas_mes * 11               # só 11 meses de trabalho, tendo em vista férias
        (sal_anual * (1 + fgts)) / horas_ano},

      # Autônomo / MEI até limite mensal = desconta DAS fixo
      VD4009 != "1" & !is.na(horas_mes) & VD4016 <= mei_mensal ~
        pmax((VD4016 - mei_das) / horas_mes, 0), #pmax garante que não vem valor negativo

      # Autônomo / CNPJ acima do limite MEI = aplica regra do INSS
      VD4009 != "1" & !is.na(horas_mes) & VD4016 > mei_mensal ~ {
        if (opcao_inss == "fixed_min_wage_11pct") {
          inss_val <- sal_min * inss_valor
          pmax((VD4016 - inss_val) / horas_mes, 0)} 
        else {
          pmax((VD4016 * (1 - inss_valor)) / horas_mes, 0)}},
      TRUE ~ NA_real_  ))

# --- categorias de emprego ---
ocupados3$variables <- ocupados3$variables |>
  mutate(categoria_emprego = case_when(
    VD4009 == "01" ~ "formal",
    VD4009 == "09" & V4019 == "1" ~ "cnpj",
    TRUE ~ "outros"))

rm(ocupados)

ggplot(ocupados3, aes(x = "", y = salario_hora)) +
  geom_boxplot(outlier.color = "red") +
  labs(title = "Boxplot de Salário-Hora", y = "Horas") +
  theme_minimal()  

#retirar outliers de salario_hora
ocupados3 <- ocupados3 |>
  filter(salario_hora >= quantile(salario_hora, 0.01, na.rm = TRUE),
         salario_hora <= quantile(salario_hora, 0.99, na.rm = TRUE))

4 Análise Descritiva

4.1 Iniciais - Por Categoria

4.1.1 Salário-Hora

resumo_sal <- ocupados3 |>
  group_by(categoria_emprego) |>
  summarise(n = survey_total(na.rm = TRUE),
            mediana = survey_median(salario_hora, na.rm = TRUE),
            .groups = "drop") |>
  mutate(n_label = paste0("n = ", scales::comma(n)),
         mediana_label = paste0("Mediana = R$ ", round(mediana, 2)))

ggplot(ocupados3, aes(x = categoria_emprego, y = salario_hora, fill = categoria_emprego)) +
  geom_boxplot(alpha = 0.7, outlier.colour = "red") +
  scale_y_continuous(labels = scales::dollar_format(prefix = "R$ ")) +
  coord_cartesian(ylim = c(0, 50)) +
  labs(title = "Distribuição do Salário-Hora por Categoria de Emprego \n (até R$ 50)",
       x = "Categoria de Emprego",
       y = "Salário-Hora") +
  theme_minimal(base_size = 14) +
  theme(legend.position = "none") +
  geom_text(data = resumo_sal,
            aes(x = categoria_emprego, y = 5, label = n_label), # coloca o texto no topo
            inherit.aes = FALSE,
            vjust = -0.5,
            size = 4) +
  geom_text(data = resumo_sal,
            aes(x = categoria_emprego, y = mediana, label = round(mediana, 1)),
            inherit.aes = FALSE,
            vjust = -0.5,
            fontface = "bold",
            size = 4,
            color = "blue")

4.1.2 Escolaridade

ocupados3 <- ocupados3 |> 
  mutate(escolaridade = factor(VD3004,
                               levels = 1:7,
                               labels = c("Sem instrução",
                                          "Fundamental incompleto",
                                          "Fundamental completo",
                                          "Médio incompleto",
                                          "Médio completo",
                                          "Superior incompleto","Superior completo")))


resumo_escolaridade <- ocupados3 |>
  filter(!is.na(escolaridade)) |>
  group_by(categoria_emprego, escolaridade) |>
  summarise(n = survey_total(na.rm = TRUE)) |>
  group_by(categoria_emprego) |>
  mutate(freq_relativa = n / sum(n) * 100)

ggplot(resumo_escolaridade, aes(x = categoria_emprego, y = freq_relativa, fill = escolaridade)) +
  geom_bar(stat = "identity", position = "stack") +
  scale_y_continuous(labels = scales::percent_format(scale = 1)) +
  labs(title = "Distribuição Percentual da Escolaridade \n por Categoria de Emprego",
       x = "Categoria de Emprego",
       y = "Frequência Relativa (%)",
       fill = "Faixa de Escolaridade") +
  theme_minimal()

4.1.3 E a questão regional?

ufs_labels <- c(
  "11" = "RO", "12" = "AC", "13" = "AM", "14" = "RR", "15" = "PA",
  "16" = "AP", "17" = "TO", "21" = "MA", "22" = "PI", "23" = "CE",
  "24" = "RN", "25" = "PB", "26" = "PE", "27" = "AL", "28" = "SE",
  "29" = "BA", "31" = "MG", "32" = "ES", "33" = "RJ", "35" = "SP",
  "41" = "PR", "42" = "SC", "43" = "RS", "50" = "MT", "51" = "MS",
  "52" = "GO", "53" = "DF")

ocupados3 <- ocupados3 |>
  mutate(uf_sigla = ufs_labels[as.character(UF)])

# Filtro apenas para duas categorias para não ficar poluído
dados_plot_uf <- ocupados3 |>
  filter(categoria_emprego %in% c("formal", "cnpj")) |>
  mutate(UF_label = factor(ufs_labels[as.character(UF)], levels = ufs_labels))

# Selecionar só UFs que tenham os dois grupos: formal e cnpj
ufs_completas <- dados_plot_uf |>
  group_by(UF_label, categoria_emprego) |>
  summarise(contagem = n(), .groups = 'drop') |>
  pivot_wider(names_from = categoria_emprego, values_from = contagem, values_fill = 0) |>
  filter(formal > 0 & cnpj > 0) |>
  pull(UF_label)

# Filtrar pelos UFs completos e criar fator para categoria de emprego
dados_plot_uf <- dados_plot_uf |>
  filter(UF_label %in% ufs_completas) |>
  mutate(categoria_emprego = factor(categoria_emprego, levels = c("formal", "cnpj")))

# Calcular medianas por UF e categoria de emprego para mostrar no gráfico
medianas_uf <- dados_plot_uf |>
  group_by(UF_label, categoria_emprego) |>
  summarise(mediana = survey_median(salario_hora, na.rm = TRUE), .groups = 'drop')

# Boxplot com facet_wrap por UF
ggplot(dados_plot_uf, aes(x = categoria_emprego, y = salario_hora, fill = categoria_emprego)) +
  geom_boxplot(alpha = 0.7, outlier.colour = "red") +
  geom_text(data = medianas_uf, 
            aes(x = categoria_emprego, 
                y = mediana, 
                label = paste0("R$ ", round(mediana, 2))),
            vjust = -0.5, color = "black", fontface = "bold", size = 4) +
  scale_y_continuous(labels = scales::dollar_format(prefix = "R$ ")) +
  coord_cartesian(ylim = c(0, 75)) +
  labs(title = "Distribuição do Salário-Hora por Categoria de Emprego e UF\n(visualização até R$ 75)",
       x = "Categoria de Emprego",
       y = "Salário-Hora") +
  facet_wrap(~ UF_label, scales = "free_x", ncol = 5) +
  theme_minimal(base_size = 14) +
  theme(legend.position = "none")

4.1.4 Contribuição Previdenciária

resumo_contribuinte <- ocupados3 |>
  filter(categoria_emprego %in% c("formal", "cnpj")) |>
  mutate(contribuinte_cat = case_when(VD4012 == 1 ~ "Contribuinte",
                                      VD4012 == 2 ~ "Não contribuinte",
                                      TRUE ~ "Não aplicável")) |>
  group_by(categoria_emprego, contribuinte_cat) |>
  summarise(n = survey_total(na.rm = TRUE)) |>
  group_by(categoria_emprego) |>
  mutate(prop = n / sum(n)) |>
  ungroup()

ggplot(resumo_contribuinte, aes(x = categoria_emprego, y = prop, fill = contribuinte_cat)) +
  geom_col(position = "fill") +
  scale_y_continuous(labels = scales::percent_format()) +
  labs(title = "Proporção de Contribuintes por Categoria de Emprego",
       x = "Categoria de Emprego",
       y = "Proporção (%)",
       fill = "Categoria de Contribuinte") +
  theme_minimal()

4.2 Grupo-Base - Profissões

Para a análise das profissões foram criados três objetos.

tab_total que sumariza o total por profissão (grupo_base).

tab_formal que sumariza por trabalhadores formais

tab_cnpj que sumariza por trabalhadores por conta-própria que possuem CNPJ

Foram criadas duas variáveis auxiliares para análise: a dif_perc_cnpj_formal que consiste em verificar a relação percentual da média entre o salário-hora do CNPJ e o salário-hora do formal e a *dif_perc_medianas_cnpj_formal** para verificar a mesma relação, mas a partir da mediana.

Uma análise adequada a partir das profissões deve considerar as duas variáveis em conjunto.

tab_total <- ocupados3 |> 
  group_by(grupo_base) |> 
  summarise(total = survey_total(na.rm = TRUE),
            salario_hora = survey_mean(salario_hora, na.rm = TRUE))

tab_total <- tab_total |> 
  mutate(total_cv = 100 * total_se / total)

tab_formal <- ocupados3 |> 
  filter(categoria_emprego == "formal") |>
  group_by(grupo_base) |> 
  summarise(total_formal = survey_total(),
            salario_hora_formal = survey_mean(salario_hora, na.rm = TRUE, vartype = "se"))

tab_formal <- tab_formal |> 
  mutate(total_formal_cv = 100 * total_formal_se / total_formal)

tab_cnpj <- ocupados3 |> 
  filter(categoria_emprego == "cnpj") |>
  group_by(grupo_base) |> 
  summarise(total_cnpj = survey_total(),
            salario_hora_cnpj = survey_mean(salario_hora, na.rm = TRUE, vartype = "se"))

tab_cnpj <- tab_cnpj |> 
  mutate(total_cnpj_cv = 100 * total_cnpj_se / total_cnpj)

tab_geral <- tab_total |> 
  full_join(tab_formal, by = "grupo_base") |> 
  full_join(tab_cnpj, by = "grupo_base")

#Variável Auxiliar - Diferença entre CNPJ e Formal
#Se CNPJ > Formal = positivo
tab_geral <- tab_geral |> 
  mutate(dif_perc_cnpj_formal = 100 * (salario_hora_cnpj - salario_hora_formal) / salario_hora_formal)

# Indicadores de precisão - Coeficiente de variação acima de 25%
tab_geral <- tab_geral |> 
  mutate(low_precision_formal = total_formal_cv > 25, #indicadores de baixa precisão - coeficiente de variação acima de 25%
    low_precision_cnpj = total_cnpj_cv > 25,
    total_formal_exib = if_else(low_precision_formal, "NS", format(round(total_formal, 0), big.mark = ",")),
    total_cnpj_exib = if_else(low_precision_cnpj, "NS", format(round(total_cnpj, 0), big.mark = ",")),
    salario_hora_formal_exib = if_else(low_precision_formal, "NS", format(round(salario_hora_formal, 2))),
    salario_hora_cnpj_exib = if_else(low_precision_cnpj, "NS", format(round(salario_hora_cnpj, 2))),
    dif_perc_cnpj_formal = if_else(
      low_precision_formal | low_precision_cnpj,
      NA_real_,
      100 * (salario_hora_cnpj - salario_hora_formal) / salario_hora_formal))

Medianas

medianas_formal <- ocupados3 |> 
  filter(categoria_emprego == "formal") |>
  group_by(grupo_base) |> 
  summarise(mediana_salario_hora_formal = survey_median(salario_hora, na.rm = TRUE, vartype = "se"),
            mediana_salario_hora_formal_se = attr(mediana_salario_hora_formal, "se"))

medianas_cnpj <- ocupados3 |> 
  filter(categoria_emprego == "cnpj") |>
  group_by(grupo_base) |> 
  summarise(mediana_salario_hora_cnpj = survey_median(salario_hora, na.rm = TRUE, vartype = "se"),
            mediana_salario_hora_cnpj_se = attr(mediana_salario_hora_cnpj, "se"))

# Indicadores de precisão para medianas
tab_geral <- tab_geral |> 
  left_join(medianas_formal, by = "grupo_base") |>
  left_join(medianas_cnpj, by = "grupo_base") |>
  mutate(
    # CV das medianas
    mediana_formal_cv = 100 * mediana_salario_hora_formal_se / mediana_salario_hora_formal,
    mediana_cnpj_cv = 100 * mediana_salario_hora_cnpj_se / mediana_salario_hora_cnpj,
    
    # Indicadores de baixa precisão
    low_precision_median_formal = mediana_formal_cv > 25,
    low_precision_median_cnpj = mediana_cnpj_cv > 25,
    
    # Exibição condicional para medianas
    mediana_salario_hora_formal_exib = if_else(low_precision_median_formal, "NS", format(round(mediana_salario_hora_formal, 2))),
    mediana_salario_hora_cnpj_exib = if_else(low_precision_median_cnpj, "NS", format(round(mediana_salario_hora_cnpj, 2))),
    
    # Diferença % medianas condicionada à precisão
    dif_perc_medianas_cnpj_formal = if_else(
      low_precision_median_formal | low_precision_median_cnpj,
      NA_real_,
      100 * (mediana_salario_hora_cnpj - mediana_salario_hora_formal) / mediana_salario_hora_formal))

4.3 Tabela de Resultados - Grupo-Base de Profissões

A tabela abaixo apresenta os resultados a partir do dado mais granular: as profissões.

É importante fazer uma leitura que considere os números ponderados, as médias e também as medianas.

tab_geral_sem_se <- tab_geral |> 
  select(
    grupo_base,
    total,
    total_formal,
    total_cnpj,
    salario_hora,
    salario_hora_formal,
    salario_hora_cnpj,
    dif_perc_cnpj_formal,
    mediana_salario_hora_formal_exib,
    mediana_salario_hora_cnpj_exib,
    dif_perc_medianas_cnpj_formal) |> 
  arrange(desc(total_cnpj))

colnames(tab_geral_sem_se) <- c(
  "Profissão",
  "Total da Profissão",
  "Total Formal",
  "Total CNPJ",
  "Salário-hora - Profissão (Média)",
  "Salário-hora Formal (Média)",
  "Salário-hora CNPJ (Média)",
  "Dif. % Média CNPJ vs. Formal",
  "Mediana Formal",
  "Mediana CNPJ",
  "Dif. % Mediana CNPJ vs. Formal")

# Atualizar o datatable para exibir e formatar as novas colunas
datatable(tab_geral_sem_se,
  options = list(
    pageLength = 20,
    columnDefs = list(
      list(targets = 0, createdCell = JS("function(td){$(td).css({'border-right': '4px double #333'});}")),
      list(targets = 3, createdCell = JS("function(td){$(td).css({'border-right': '4px double #333'});}")),
      list(targets = 7, createdCell = JS("function(td){$(td).css({'border-right': '4px double #333'});}")),
      list(targets = 10, createdCell = JS("function(td){$(td).css({'border-right': '4px double #333'});}"))
      #list(targets = 10, createdCell = JS("function(td){$(td).css({'border-right': '4px double #333'});}"))
    )
  ),
  rownames = FALSE,
  class = 'cell-border stripe',
  caption = htmltools::tags$caption(
    style = 'caption-side: top; text-align: left; font-size: 16px; font-weight: bold;',
    "Tabela Comparativa por Profissões com Médias e Medianas",
    htmltools::tags$div(
      style = 'font-size: 12px; font-weight: normal; color: gray; margin-top: 4px;',
      "Nota: Valores com coeficiente de variação acima de 25% foram retirados --- A ordenação-padrão está em ordem decrescente para o Total de CNPJ" 
    ))
) |> 
  formatCurrency(columns = c('Salário-hora - Profissão (Média)', 'Salário-hora Formal (Média)', 'Salário-hora CNPJ (Média)', 
                            'Dif. % Média CNPJ vs. Formal', 'Mediana Formal', 
                            'Mediana CNPJ', 'Dif. % Mediana CNPJ vs. Formal'), currency = "", digits = 2, interval = 3, mark = "."
  ) |> 
  formatRound(
    columns = c('Total da Profissão', 'Total Formal', 'Total CNPJ'),
    digits = 0, mark = "."
  ) |>
  formatStyle(
    columns = 1:ncol(tab_geral_sem_se),
    fontFamily = "Calibri, sans-serif")

4.3.1 Conclusões Parciais

Uma análise a partir das profissões nos permite avaliar distintas questões:

Em números absolutos, as 5 profissões com mais trabalhadores são:

• Escriturário Geral (ou seja, trabalhadores de apoio administrativo geral)
• Balconistas e vendedores de lojas
• Trabalhadores dos serviços domésticos em geral
  
• Comerciantes de lojas
• Trabalhadores de limpeza de interior de edifícios, escritórios, hotéis e outros estabelecimentos

Em números absolutos, as 5 profissões com mais trabalhadores formalizados são:

• Balconistas e vendedores de lojas
• Escriturários gerais
  
• Trabalhadores de limpeza de interior de edifícios, escritórios, hotéis e outros estabelecimentos
  
• Condutores de caminhões pesados
  
• Guardas de segurança

Em números absolutos, as 5 profissões com mais trabalhadores por cont-própria que possuem CNPJ são:

• Comerciantes de lojas
• Cabeleireiros
• Especialistas em tratamento de beleza e afins
• Condutores de automóveis, taxis e caminhonetes
  
• Pedreiros

Uma análise da tabela com uma desagregação muito grande deve considerar os valores de salário-hora a partir do número ponderado de profissionais.

Assim, as cinco maiores médias de salário-hora de formalizados, que considere somente acima de 20 mil profissionais, são as seguintes profissões:

• Médicos especialistas - R$ 54.14
  
• Engenheiros mecânicos - R$ 45.17
  
• Dirigentes de serviços de tecnologia da informação e comunicações - R$ 45.11
  
• Engenheiros industriais e de produção - R$ 42.21
• Engenheiros não classificados anteriormente - R$ 38.74

Por sua vez, as cinco maiores médias de salário-hora de conta-própria com CNPJ, que considere somente acima de 20 mil profissionais, são as seguintes profissões:

• Médicos especialistas - R$ 50.89
• Analistas de sistemas - R$ 43.11
• Advogados e juristas - R$ 34.99
• Engenheiros civis - R$ 34.46
• Analista de Gestão e Administração - R$ 33.51

Dado o objetivo do estudo, entretanto, é relevante direcionar o olhar para a diferença entre profissionais de ocupação semelhante. Nesse sentido, é importante sempre olhar a diferença a partir de uma análise coordenada de média e mediana. A partir do qual pode-se concluir:

As maiores diferenças em que é o PJ recebe mais que o formalizado encontram-se nas seguintes profissões

• Agentes imobiliários - onde um PJ (~89 mil) recebe uma média de R$ 30,50 e o CLT (~29 mil) recebe uma média de R$ 16,09. A análise da mediana também indica uma diferença relevante com formal a R$ 11.69 e PJ a R$ 24,72. É importante destacar, entretanto que a maioria dos profissionais (cerca de 63%) não estão em nenhuma dessas duas classificações.

Análise semelhante deve ser feita para Especialistas em tratamento de beleza e afins que apresenta diferença significativa entre as duas categorias (82% a mais para PJ na média e 56% a mais na mediana), mas tem-se cerca de 75% de todos os profissionais fora dessas duas categorias.

As demais diferenças que podem ser destacadas estão nas seguintes profissões:

• Trabalhadores de agências de viagem
  
• Condutores de automóveis, taxis e caminhonetes
  
• Representantes comerciais
• Lavadores de veículos
• Condutores de caminhões pesados

4.4 Grupamentos de Atividade

vd4010_labels <- c(
  "01" = "Agricultura, pecuária, produção florestal, pesca e aquicultura",
  "02" = "Indústria geral",
  "03" = "Construção",
  "04" = "Comércio, reparação de veículos automotores e motocicletas",
  "05" = "Transporte, armazenagem e correio",
  "06" = "Alojamento e alimentação",
  "07" = "Informação, comunicação e atividades financeiras, imobiliárias, profissionais e administrativas",
  "08" = "Administração pública, defesa e seguridade social",
  "09" = "Educação, saúde humana e serviços sociais",
  "10" = "Outros Serviços",
  "11" = "Serviços domésticos",
  "12" = "Atividades mal definidas")

ocupados3 <- ocupados3 |>
  mutate(VD4010_label = vd4010_labels[as.character(VD4010)])

resumo_salario <- ocupados3 |>
  group_by(VD4010_label, categoria_emprego) |>
  summarise(salario_medio = mean(salario_hora, na.rm = TRUE),
            n = survey_total(na.rm = TRUE)) |>
  ungroup() |>
  mutate(n_label = paste0("n = ", scales::comma(round(as.numeric(n)))))

# Gráfico de barras com facet
ggplot(resumo_salario, aes(x = categoria_emprego, y = salario_medio, fill = categoria_emprego)) +
  geom_col(position = "dodge") +
  geom_text(aes(label = paste0("n=", scales::comma(round(n))), y = salario_medio + max(salario_medio, na.rm = TRUE) * 0.03),
            position = position_dodge(width = 0.9), size = 4, color = "black") +
  facet_wrap(~ VD4010_label, scales = "free_y", ncol = 3,
             labeller = labeller(VD4010_label = label_wrap_gen(width = 40))) +
  scale_y_continuous(labels = scales::dollar_format(prefix = "R$ ")) +
  labs(title = "Média do Salário-Hora por Categoria de Emprego e Atividade\n(n ponderado por barra)",
       x = "Categoria de Emprego",
       y = "Salário-Hora Médio") +
  theme_minimal(base_size = 13) +
  theme(axis.text.x = element_text(angle = 0, hjust = 0.5), legend.position = "none")

max_salario <- 50 # ponto próximo da borda direita do gráfico, tendo em vista outliers

#Boxplot
ggplot(ocupados3, aes(x = salario_hora, y = categoria_emprego, fill = categoria_emprego)) +
  geom_boxplot(alpha = 0.7, outlier.colour = "red") +
  facet_wrap(~ VD4010_label, scales = "free_y", ncol = 3,
             labeller = labeller(VD4010_label = label_wrap_gen(width = 40))) +
  scale_x_continuous(labels = dollar_format(prefix = "R$ ")) +
  coord_cartesian(xlim = c(0, 60)) +
  labs(
    title = "Distribuição do Salário-Hora por Categoria de Emprego e Atividade\n(Limite de visualização: R$ 60)",
    x = "Salário-Hora",
    y = "Categoria de Emprego") +
  theme_minimal(base_size = 13) +
  theme(legend.position = "none") +
  geom_text(
    data = resumo_salario,
    mapping = aes(x = max_salario, y = categoria_emprego, label = n_label),
    inherit.aes = FALSE,
    hjust = 0,
    size = 3)

4.5 Grupamentos Ocupacionais

# Dicionário dos rótulos das ocupações - VD4011
vd4011_labels <- c(
  "01" = "Diretores e gerentes",
  "02" = "Profissionais das ciências e intelectuais",
  "03" = "Técnicos e profissionais de nível médio",
  "04" = "Trabalhadores de apoio administrativo",
  "05" = "Trabalhadores dos serviços, vendedores dos comércios e mercados",
  "06" = "Trabalhadores qualificados da agropecuária, florestais, da caça e da pesca",
  "07" = "Trabalhadores qualificados, operários e artesões da construção, das artes mecânicas e outros ofícios",
  "08" = "Operadores de instalações e máquinas e montadores",
  "09" = "Ocupações elementares",
  "10" = "Membros das forças armadas, policiais e bombeiros militares",
  "11" = "Ocupações maldefinidas")

ocupados3 <- ocupados3 |>
  mutate(VD4011_label = vd4011_labels[as.character(VD4011)] |> factor(levels = vd4011_labels))

resumo_salario2 <- ocupados3 |>
  group_by(VD4011_label, categoria_emprego) |>
  summarise(salario_medio = mean(salario_hora, na.rm = TRUE),
            n = survey_total(na.rm = TRUE)) |>
  ungroup() |>
  mutate(n_label = paste0("n = ", scales::comma(round(as.numeric(n)))))

# Gráfico de barras com facet
ggplot(resumo_salario2, aes(x = categoria_emprego, y = salario_medio, fill = categoria_emprego)) +
  geom_col(position = "dodge") +
  geom_text(aes(label = paste0("n=", scales::comma(round(n))), y = salario_medio + max(salario_medio, na.rm = TRUE) * 0.03),
            position = position_dodge(width = 0.9), size = 4, color = "black") +
  facet_wrap(~ VD4011_label, scales = "free_y", ncol = 3,
             labeller = labeller(VD4011_label = label_wrap_gen(width = 40))) +
  scale_y_continuous(labels = scales::dollar_format(prefix = "R$ ")) +
  labs(title = "Média do Salário-Hora por Categoria de Emprego e Grupo Ocupacional\n(n ponderado por barra)",
       x = "Grupo Ocupacional",
       y = "Salário-Hora Médio") +
  theme_minimal(base_size = 13) +
  theme(axis.text.x = element_text(angle = 0, hjust = 0.5), legend.position = "none")

# Plot boxplot horizontal com n como texto
ggplot(ocupados3, aes(x = salario_hora, y = categoria_emprego, fill = categoria_emprego)) +
  geom_boxplot(alpha = 0.7, outlier.colour = "red") +
  facet_wrap(~ VD4011_label, scales = "free_y", ncol = 3,
             labeller = labeller(VD4011_label = label_wrap_gen(width = 40))) +
  scale_x_continuous(labels = dollar_format(prefix = "R$ ")) +
  coord_cartesian(xlim = c(0, 60)) +
  labs(title = "Distribuição do Salário-Hora por Categoria de Emprego e Grupamento Ocupacional\n(Limite de visualização: R$ 100)",
       x = "Salário-Hora",
       y = "Categoria de Emprego") +
  theme_minimal(base_size = 13) +
  theme(legend.position = "none") +
  geom_text(data = resumo_salario2,
            mapping = aes(x = max_salario, y = categoria_emprego, label = n_label),
            inherit.aes = FALSE,
            hjust = 0,
            size = 3)

5 Análise Inferencial

5.1 Mincer

Existe uma diferença salarial entre a condição “Conta-Própria com CNPJ” e “Formal - CLT”.

A pergunta que se coloca, entretanto, é: o quanto essa diferença salarial deve-se à condição PJ/Formal?

Para isso, vamos fazer um modelo de regressão a partir das equações de Mincer.

5.1.1 Recodificação de Variáveis

No código abaixo foi criado o log(salario_hora), colocados os rótulos para sexo e raça/cor.

Além disso, será criada a variável idade_centralizada que consiste em

\[ \text{Idade_Centralizada}_{i} = Idade_{i} - \overline{Idade} \] Essa variável permite diminuir a multicolinearidade entre idade e idade^2

Para uma aproximação da variavel experiência será feito o seguinte cálculo

\[ \text{Experiência} = Idade - Anos de Escolaridade - 5 \]

ocupados3 <- ocupados3 |>
  mutate(
    # log do salário-hora (apenas se válido)
    ln_salario_hora = if_else(salario_hora > 0, log(salario_hora), NA_real_),
    
    # Sexo
    sexo = factor(V2007,
      levels = c(1, 2),
      labels = c("Homem", "Mulher")),
    
    # Raça
    raca = factor(V2010,
      levels = c(1, 2, 3, 4, 5, 9),
      labels = c("Branca", "Preta", "Amarela", "Parda", "Indígena", "Ignorado")),
    
    # Idade
    idade = as.numeric(V2009),

    # Anos de escolaridade (converter de character p/ numeric)
    anos_escolaridade = as.numeric(VD3005)) |>
  mutate(
    # Idade centralizada pela média
    idade_centralizada = idade - mean(idade, na.rm = TRUE),

    # Experiência potencial
    experiencia = idade - anos_escolaridade - 5,
    
    # Experiência centralizada pela média
    experiencia_c = experiencia - mean(experiencia, na.rm = TRUE))

5.1.2 Regressões

A equação de Mincer pode ser expressa como:

\[ \ln(w_i) = \beta_0 + \beta_1 Escolaridade_i + \beta_2 Experiência_i + \beta_3 Experiência_i^2 + \varepsilon_i \] O modelo estimado é dado por:

\[ \begin{aligned} \ln(\text{SalárioHora}_i) =\ & \beta_0 + \beta_1 \text{Categoria de Emprego}_i + \beta_2 \overline{Experiência}_i \\ &+ \beta_3 \overline{Experiência}_i^2 + \beta_4 \text{Raça}_i \\ &+ \beta_5 \text{Sexo}_i + \beta_6 \text{SetorAtividade}_i + \varepsilon_i \end{aligned} \]

Foram realizados dois testes:

1 - O primeiro que utiliza idade_centralizada + idade_centralizada^2 2 - O segundo que utiliza o proxy experiencia_centralizada

Em ambos os modelos o poder explicativo ficou em 25%. Assim, para manter fidelidade com o modelo de Mincer optou-se por seguir

#Não substituir UF e VD4010 pelos labels porque o tempo de processamento aumenta demais
library(jtools)
modelo_mincer <- summ(
  svyglm(ln_salario_hora ~ categoria_emprego + anos_escolaridade + experiencia_c + I(experiencia_c^2) + raca + sexo + VD4010,
  design = ocupados3,
  family = gaussian()),
  digits = 2, confint = FALSE)

coefs <- modelo_mincer$coeftable |> 
  as.data.frame() |> 
  rownames_to_column("term")

# filtrar para p <= 0.05 e sem o intercepto
sig_coefs <- coefs |>
  filter(p <= 0.05, term != "(Intercept)") |>
  pull(term)

plot_summs(modelo_mincer, coefs = sig_coefs)

#Referências
#CNPJ
#Escolaridade - Sem Instrução
#Raça - Branca
#Sexo - Masculino
#VD4010 1 - Agricultura, pecuária, produção florestal, pesca e aquicultura 

#Percentual

library(broom)
resultados <- broom::tidy(modelo_mincer) |>
  mutate(efeito_pct = (exp(estimate) - 1) * 100,   # transformar log em %
         signif = ifelse(p.value < 0.05, "Significativo", "Não significativo"))

#apenas significativos
resultados_sig <- resultados |>
  filter(signif == "Significativo") |>
  select(term, estimate, efeito_pct, std.error, statistic, p.value)

# renomear colunas
resultados_sig <- resultados_sig |>
  rename(Variável = term,
         Coef_Log = estimate,
         `Efeito (%)` = efeito_pct,
         `Erro Padrão` = std.error,
         `t valor` = statistic,
         `p-valor` = p.value)

library(kableExtra)

resultados_sig |>
  kable(digits = 2, caption = "Efeitos percentuais significativos no salário-hora (modelo Mincer)") |>
  kable_styling(full_width = FALSE, position = "center")

Efeitos percentuais significativos no salário-hora (modelo Mincer)

Variável	Coef_Log	Efeito (%)	Erro Padrão	t valor	p-valor
(Intercept)	1.84	532.77	0.02	109.93	0.00
categoria_empregoformal	-0.19	-17.37	0.01	-18.79	0.00
categoria_empregooutros	-0.26	-22.95	0.01	-24.87	0.00
anos_escolaridade	0.08	8.22	0.00	89.07	0.00
experiencia_c	0.01	1.16	0.00	68.89	0.00
I(experiencia_c^2)	0.00	-0.02	0.00	-23.57	0.00
racaPreta	-0.18	-16.33	0.01	-28.67	0.00
racaParda	-0.18	-16.66	0.00	-44.07	0.00
racaIndígena	-0.20	-18.26	0.03	-6.80	0.00
sexoMulher	-0.19	-16.95	0.00	-48.13	0.00
VD401002	0.21	22.77	0.01	19.07	0.00
VD401003	0.17	19.07	0.01	16.11	0.00
VD401004	0.08	8.09	0.01	7.02	0.00
VD401005	0.20	21.61	0.01	16.56	0.00
VD401006	0.04	3.84	0.01	3.11	0.00
VD401007	0.31	36.41	0.01	28.22	0.00
VD401008	0.54	72.18	0.01	44.64	0.00
VD401009	0.40	49.58	0.01	35.17	0.00
VD401010	0.13	14.20	0.01	10.46	0.00
VD401011	0.02	2.51	0.01	2.05	0.04

Colinearidade

Muita colinearidade para UF

# usa os dados que estão dentro do design survey
data_svy <- ocupados3$variables

# monta a matrix de regressors como o modelo que você citou
mm <- model.matrix(~ categoria_emprego + anos_escolaridade + experiencia_c + I(experiencia_c^2) + sexo + raca + VD4010,
                   data = data_svy)

# remove a coluna do intercepto
preds <- as.data.frame(mm[, -1, drop = FALSE])

# obtém os pesos de amostra do design (tipo "sampling")
wt <- weights(ocupados3, type = "sampling")

# função para R2 ponderado (usaremos em lm com weights)
weighted_R2 <- function(y, yhat, w) {
  # soma dos quadrados dos resíduos ponderados
  sse <- sum(w * (y - yhat)^2, na.rm = TRUE)
  # soma total ponderada
  ybar_w <- stats::weighted.mean(y, w, na.rm = TRUE)
  sst <- sum(w * (y - ybar_w)^2, na.rm = TRUE)
  1 - sse / sst
}

# calcula VIF para cada coluna da matrix de regressors
vif_table <- data.frame(variable = character(), R2 = numeric(), VIF = numeric(), stringsAsFactors = FALSE)

for (j in seq_len(ncol(preds))) {
  y <- preds[[j]]
  Xother <- preds[, -j, drop = FALSE]

  # monta data frame temporário
  df_tmp <- cbind(y = y, Xother)

  # regressão ponderada com lm (resposta y contra os demais)
  fit <- lm(y ~ ., data = df_tmp, weights = wt)

  # prediz e calcula R2 ponderado
  yhat <- predict(fit, newdata = df_tmp)
  R2w <- weighted_R2(y, yhat, wt)

  vif_j <- 1 / (1 - R2w)
  vif_table <- rbind(vif_table, data.frame(variable = colnames(preds)[j], R2 = R2w, VIF = vif_j))
}

# ordena e mostra
vif_table[order(-vif_table$VIF), ]

O modelo explica cerca de 32% da diferença salarial.

Formais ganham, em média, cerca de 19% menos que o conta-própria com CNPJ, mantendo escolaridade, idade, raça, sexo, e setor econômico fixos.

Mas qual variável conta mais?

5.2 Comparar modelos

Testando função para ver o incremento

# Função compacta para calcular pseudo-R² para svyglm (Gaussian)
pseudo_r2 <- function(model) {
  y <- model$y
  yhat <- model$fitted.values
  ss_res <- sum((y - yhat)^2)
  ss_tot <- sum((y - mean(y))^2)
  1 - ss_res / ss_tot}

# Ajuste incremental, reutilizando o objeto de design e limpando explicitamente
r2s <- numeric(0)

formulas <- list(
  ln_salario_hora ~ categoria_emprego,
  ln_salario_hora ~ categoria_emprego + sexo,
  ln_salario_hora ~ categoria_emprego + sexo + raca,
  ln_salario_hora ~ categoria_emprego + sexo + raca + VD4010,
  ln_salario_hora ~ categoria_emprego + sexo + raca + VD4010 + anos_escolaridade,
  ln_salario_hora ~ categoria_emprego + sexo + raca + VD4010 + anos_escolaridade + experiencia_c,
  ln_salario_hora ~ categoria_emprego + sexo + raca + VD4010 + anos_escolaridade + experiencia_c + I(experiencia_c^2)
)

for (f in formulas) {
  modelo <- svyglm(f, design = ocupados3, family = gaussian())
  r2s <- c(r2s, pseudo_r2(modelo))
  rm(modelo) # remove objeto grande
  gc() # força o garbage collector a liberar memória
}

names(r2s) <- c("categoria_emprego", "+sexo", "+raça", "+grupo_ativ", "+escolaridade", "+experiencia", "+experiencia2")
print(r2s)

## categoria_emprego             +sexo             +raça       +grupo_ativ 
##        0.01351587        0.01244125        0.05583003        0.20438819 
##     +escolaridade      +experiencia     +experiencia2 
##        0.28018363        0.31889390        0.32297245

5.2.1 Conclusões Parciais

Ainda que exista uma diferença entre as categorias de emprego percebe-se que o fator que mais influencia o salário-hora é a escolaridade e o grupo de atividade.

A categoria emprego sozinha explica cerca de 1% da variação salarial.

5.3 Propensity Score Matching

Propensity Score Matching cria um grupo de controle artificial para comparar com o grupo de tratamento, com o objetivo de eliminar vieses.

A escolha por esse teste permite responder à pergunta:

“Dados dois trabalhadores com perfis semelhantes (mesma escolaridade, experiência, sexo e raça), como fica a diferença salário-hora entre aquele formalizado e conta-própria com CNPJ?”

library(MatchIt)

# Variável de tratamento
ocupados_psm <- ocupados3$variables |>
  filter(categoria_emprego %in% c("formal", "cnpj")) |>
  mutate(formal = if_else(categoria_emprego == "formal", 1, 0),
         formal = factor(formal, levels = c(0, 1), labels = c("Conta própria", "Formal")))

ocupados_psm <- as.data.frame(ocupados_psm)

# Matching com logit 
# Foi retirado VD4010 pois demorava cerca de 40 minutos.
system.time({
psm_model <- matchit(formal ~ sexo + raca + anos_escolaridade + experiencia_c + I(experiencia_c^2),
                     data = ocupados_psm,
                     method = "nearest",
                     distance = "logit")
})

##   usuário   sistema decorrido 
##     47.97      2.63     51.64

matched_data <- match.data(psm_model)

library(cobalt)
love.plot(psm_model, binary = "std")  # gráfico de diferenças padronizadas

5.3.1 Teste t

t_test <- t.test(salario_hora ~ formal, data = matched_data, weights = weights)
t_test

## 
##  Welch Two Sample t-test
## 
## data:  salario_hora by formal
## t = 74.234, df = 12181, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group Conta própria and group Formal is not equal to 0
## 95 percent confidence interval:
##  10.77835 11.36300
## sample estimates:
## mean in group Conta própria        mean in group Formal 
##                   20.512122                    9.441448

# Visualizar balanceamento (standardized mean differences)
#plot(psm_model, type = "jitter", interactive = FALSE)
#plot(psm_model, type = "qq")
#plot(psm_model, type = "hist")

matched_data$final_weights <- matched_data$weights * matched_data$V1028

options(survey.lonely.psu = "adjust")

design_psm <- svydesign(
  ids = ~1,  # ignora cluster
  weights = ~final_weights,
  data = matched_data
)

# Teste t ponderado
svy_test <- svyttest(salario_hora ~ formal, design = design_psm)
svy_test

## 
##  Design-based t-test
## 
## data:  salario_hora ~ formal
## t = -53.95, df = 21342, p-value < 2.2e-16
## alternative hypothesis: true difference in mean is not equal to 0
## 95 percent confidence interval:
##  -11.76878 -10.94361
## sample estimates:
## difference in mean 
##          -11.35619

O test t indica R$ 20,51 para conta própria e R$ 9,44 para o formalizado (cerca de R$ 11,07 a menos) O test t para objeto do tipo survey indica uma diferença de cerca de -R$ 11,35 para o formalizado

5.4 Quantílica com PSM

Como já foi possível visualizar, especialmente na parte de análise descritiva, as variáveis que envolvem salários tem boxplots com o terceiro quartil mais alongado e o quarto quartil com muitos outliers.

Assim, a pergunta seguinte é:

A diferença salário-hora entre formais e PJs se mantêm entre todos os níveis salariais?

Para responder essa pergunta será feita uma análise quantílica. Ou seja, o salário-hora será dividido em dez partes e será observado o comportamento diferencial entre cada decil.

library(quantreg)

taus <- seq(0.0, 1.0, by = 0.1)
rq_model_all <- rq(salario_hora ~ formal, 
                   data = matched_data, 
                   tau = taus, 
                   weights = weights)

system.time({
rq_summary <- summary(rq_model_all, se = "boot", R = 50)
})

##   usuário   sistema decorrido 
##    108.61      0.22    109.98

5.4.1 Resultados

#Sys.setlocale("LC_ALL", "pt_BR.UTF-8")

results <- do.call(rbind, lapply(rq_summary, function(s) {
  data.frame(
    tau = s$tau,
    coef = s$coefficients["formalFormal", 1],
    se   = s$coefficients["formalFormal", 2]
  )}))

#erros-padrão
results <- results |>
  mutate(lower = coef - 1.96 * se,
         upper = coef + 1.96 * se,
         type = "Quantílico")

#resultados dos dois tests t - adicionados
mean_effects <- data.frame(
  tau = c(0.5, 0.5),
  coef = c(t_test$estimate[2] - t_test$estimate[1], svy_test$estimate),     # Formal - Conta própria
  lower = c(-t_test$conf.int[2], svy_test$conf.int[1]),                     # Inverta upper (vai para lower)
  upper = c(-t_test$conf.int[1], svy_test$conf.int[2]),                     # Inverta lower (vai para upper)
  type = c("T-test", "Svy T-test"))

plot_data <- bind_rows(results, mean_effects)

g1 <- ggplot() +
  # efeitos quantílicos
  geom_line(data = results, aes(x = tau, y = coef, color = type), size = 1) +
  geom_point(data = results, aes(x = tau, y = coef, color = type), size = 2) +
  geom_ribbon(data = results, aes(x = tau, ymin = lower, ymax = upper, fill = type), alpha = 0.15) +

  # rótulos nos efeitos quantílicos
  geom_text(data = results, 
            aes(x = tau, y = coef, label = round(coef, 1), color = type), 
            vjust = -1, size = 3.5, show.legend = FALSE) +

  # médias destacadas
  geom_point(data = mean_effects, aes(x = tau, y = coef, shape = type, color = type), size = 4) +
  geom_errorbar(data = mean_effects, aes(x = tau, ymin = lower, ymax = upper, color = type), 
                width = 0.05, size = 1) +

  # rótulos nas médias
  geom_text(data = mean_effects, 
            aes(x = tau, y = coef, label = round(coef, 1), color = type), 
            vjust = -1.2, size = 4, fontface = "bold", show.legend = FALSE) +

  # linha de referência
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +

  # eixo de 0% a 100%
  scale_x_continuous(breaks = seq(0, 1, 0.1),  # decil a decil
                     labels = scales::percent) +

  labs(x = "Quantil da distribuição salarial",
       y = "Diferença (Formal - Conta própria)",
       title = "Diferença Salário-Hora (Formal vs PJ) por Decil - Usando PSM \n Efeitos médios (t-test / survey t-test) vs quantílicos") +
  theme_minimal(base_size = 14)

g1

5.4.2 Composição dos Quantis

library(Hmisc)

# Calcular os cortes de salário-hora ponderados
quantile_cuts <- Hmisc::wtd.quantile(matched_data$salario_hora,
                                     weights = matched_data$final_weights,
                                     probs = taus)

quantile_cuts

##        0%       10%       20%       30%       40%       50%       60%       70% 
##  1.284923  6.873846  7.491608  8.016084  9.065035 10.493846 12.592615 15.358741 
##       80%       90%      100% 
## 20.603497 30.218462 82.153846

# Classificar cada indivíduo no seu intervalo de quantil
matched_data <- matched_data |>
  mutate(quantil = cut(salario_hora,
                       breaks = quantile_cuts,
                       labels = c("0-10%", "10-20%", "20-30%", "30-40%", "40-50%", 
                                  "50-60%", "60-70%", "70-80%", "80-90%", "90-100%"),
                       right = TRUE,
                       include.lowest = TRUE))

dist_quantis <- matched_data |>
  group_by(quantil, formal) |>
  summarise(n_ponderado = sum(final_weights, na.rm = TRUE), .groups = "drop") |>
  filter(!is.na(quantil)) |>
  mutate(prop_total = n_ponderado / sum(n_ponderado))

g2 <- ggplot(dist_quantis, aes(x = quantil, y = prop_total, fill = formal)) +
  geom_col(position = position_dodge(width = 0.8), width = 0.7) +
  geom_text(aes(label = scales::percent(prop_total, accuracy = 0.1)),
            position = position_dodge(width = 0.8),
            vjust = -0.3, size = 3.5) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  labs(title = "Proporção de Formal vs Conta Própria \n por Quantil (com PSM)",
       x = "Quantis do salário-hora", y = "Proporção do total geral",
       fill = "Categoria de emprego") +
  theme_minimal(base_size = 14)

g2

5.4.3 Média salarial

medias_decis <- matched_data |>
  group_by(quantil, formal) |>
  summarise(media_salario = weighted.mean(salario_hora, w = final_weights, na.rm = TRUE),
            .groups = "drop")

# se os quantis estão como fator, converter para numérico para alinhar no eixo
medias_decis <- medias_decis |>
  mutate(quantil_num = as.numeric(gsub("(%|\\-).*", "", quantil)) / 100) #ponto médio

medias_decis <- medias_decis |> filter(!is.na(quantil))

g3 <- ggplot(medias_decis, aes(x = quantil, y = media_salario, fill = formal)) +
  geom_col(position = "dodge", alpha = 0.8) +
  geom_text(aes(label = round(media_salario, 1)),
            position = position_dodge(width = 0.9),
            vjust = -0.1, size = 3.5) +
  labs(
    title = "Salário-Hora Médio por Decil da Distribuição",
    subtitle = "Comparação entre trabalhadores formais e conta própria \n com PSM",
    x = "Quantil da distribuição salarial",
    y = "Salário-hora médio (R$)",
    fill = "Categoria") +
  scale_y_continuous(labels = scales::dollar_format(prefix = "R$ ", big.mark = ".", decimal.mark = ",")) +
  theme_minimal(base_size = 14)

g3

5.5 Quantílica sem PSM

5.5.1 Resultados

t_test <- t.test(salario_hora ~ formal, data = ocupados_psm, weights = V1028)
t_test

## 
##  Welch Two Sample t-test
## 
## data:  salario_hora by formal
## t = 44.486, df = 12230, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group Conta própria and group Formal is not equal to 0
## 95 percent confidence interval:
##  6.341628 6.926236
## sample estimates:
## mean in group Conta própria        mean in group Formal 
##                    20.51212                    13.87819

rq_model_all_SEMPSM <- rq(salario_hora ~ formal, 
                   data = ocupados_psm, 
                   tau = taus, 
                   weights = V1028)

#boot com R = 50 para diminuir o tempo de processamento. Diminuiu de 1h20 para 30m
system.time({
rq_summary_SEMPSM <- summary(rq_model_all_SEMPSM, se = "boot", R=50)
})

##   usuário   sistema decorrido 
##   1131.08      1.33   1135.28

results_SEMPSM <- do.call(rbind, lapply(rq_summary_SEMPSM, function(s) {
  data.frame(
    tau = s$tau,
    coef = s$coefficients["formalFormal", 1],
    se   = s$coefficients["formalFormal", 2])}))

results_SEMPSM <- results_SEMPSM |>
  mutate(lower = coef - 1.96 * se,
         upper = coef + 1.96 * se,
         type = "Quantílico")

g4 <- ggplot() +
  # efeitos quantílicos
  geom_line(data = results_SEMPSM, aes(x = tau, y = coef, color = type), size = 1) +
  geom_point(data = results_SEMPSM, aes(x = tau, y = coef, color = type), size = 2) +
  geom_ribbon(data = results_SEMPSM, aes(x = tau, ymin = lower, ymax = upper, fill = type), 
              alpha = 0.15) +

  # rótulos nos efeitos quantílicos
  geom_text(data = results_SEMPSM, 
            aes(x = tau, y = coef, label = round(coef, 1), color = type), 
            vjust = -0.3, size = 3.5, show.legend = FALSE) +

  # linha de referência
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +

  # eixo de 0% a 100%
  scale_x_continuous(breaks = seq(0, 1, 0.1),  # decil a decil
                     labels = scales::percent) +
  
  labs(x = "Quantil da distribuição salarial",
       y = "Diferença (Formal - Conta própria)",
       title = "Diferença Salário-Hora (Formal vs PJ) por Decil - Sem PSM") +
  theme_minimal(base_size = 14)

library(patchwork)
g1 / g4

5.5.2 Composição dos Quantis

quantile_cuts_puro <- Hmisc::wtd.quantile(ocupados_psm$salario_hora,
                                          weights = ocupados_psm$V1028,
                                          probs = taus)

ocupados_psm <- ocupados_psm |>
  mutate(quantil = cut(salario_hora,
                       breaks = quantile_cuts_puro,
                       labels = c("0-10%", "10-20%", "20-30%", "30-40%", "40-50%", 
                                  "50-60%", "60-70%", "70-80%", "80-90%", "90-100%"),
                       include.lowest = TRUE))

dist_quantis_puro <- ocupados_psm |>
  group_by(quantil, formal) |>
  summarise(n_ponderado = sum(V1028, na.rm = TRUE), .groups = "drop") |>
  filter(!is.na(quantil)) |>
  mutate(prop_total = n_ponderado / sum(n_ponderado))

g5 <- ggplot(dist_quantis_puro, aes(x = quantil, y = prop_total, fill = formal)) +
  geom_col(position = position_dodge(width = 0.8), width = 0.7) +
  geom_text(aes(label = scales::percent(prop_total, accuracy = 0.1)),
            position = position_dodge(width = 0.8),
            vjust = -0.3, size = 3.5) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  labs(title = "Proporção de Formal vs Conta Própria \n por Quantil (sem PSM)",
       x = "Quantis do salário-hora", y = "Proporção do total geral",
       fill = "Categoria de emprego") +
  theme_minimal(base_size = 14)


## Comparação lado a lado
g2 / g5

table(matched_data$formal)

## 
## Conta própria        Formal 
##         10672         10672

table(ocupados_psm$formal)

## 
## Conta própria        Formal 
##         10672         67340

prop.table(table(matched_data$formal))

## 
## Conta própria        Formal 
##           0.5           0.5

prop.table(table(ocupados_psm$formal))

## 
## Conta própria        Formal 
##     0.1367995     0.8632005

5.5.3 Média Salarial

medias_sempsm <- ocupados_psm |>
  group_by(quantil, formal) |>
  summarise(media_salario = weighted.mean(salario_hora, w = V1028, na.rm = TRUE),
            .groups = "drop") |>
  filter(!is.na(quantil)) |>
  mutate(tipo = "Sem PSM")


g6 <- ggplot(medias_sempsm, aes(x = quantil, y = media_salario, fill = formal)) +
  geom_col(position = "dodge", alpha = 0.8) +
  geom_text(aes(label = round(media_salario, 1)),
            position = position_dodge(width = 0.9),
            vjust = -0.1, size = 3.5) +
  labs(
    title = "Salário-Hora Médio por Decil da Distribuição",
    subtitle = "Comparação entre trabalhadores formais e conta própria \n sem PSM",
    x = "Quantil da distribuição salarial",
    y = "Salário-hora médio (R$)",
    fill = "Categoria") +
  scale_y_continuous(labels = scales::dollar_format(prefix = "R$ ", big.mark = ".", decimal.mark = ",")) +
  theme_minimal(base_size = 14)


g3 / g6

5.5.4 Conclusões Parciais

A formalidade estabiliza (menos dispersão), mas no topo implica em salários-hora consistentemente mais baixos.

O efeito não é constante: é pequeno na base e grande no topo.

Isso indica que a formalidade atua como uma “âncora” salarial: protege contra rendas muito baixas, mas também limita a chance de rendas muito altas.

No percentil 10, a diferença é positiva para formais, a medida que o salário-hora sobe a diferença se inverte, mas é gradual.

Isso confirma que a formalidade cresce ao longo da distribuição: quanto mais alto o decil, maior a perda relativa dos formais.

Já os pontos de t-test e survey t-test ficam próximos de -13 R$/h, refletindo a diferença média encontrada antes, mas demonstrando como a análise quantílica é relevante.

O resultado médio mascara a heterogeneidade: ser formal implica, em média, salário-hora menor, mas a magnitude depende fortemente da posição na distribuição.