Tendências de Estratificação Horizontal

if(!require(pacman)) install.packages("pacman")
library(pacman)

pacman::p_load(dplyr, car, psych, nnet, AER, lmtest,
               gtsummary, reshape2, ggplot2, DescTools)

# Passo 2: Carregar o banco de dados
dados <- readxl::read_xlsx("data/Banco de Dados (2006-2023) - SIICUSP (1).xlsx")

# Passo 2.1: Transformação do banco de dados
# Modifiicar nome das colunas, e dos valores das células

## Renomeando as colunas
colnames(dados)[colnames(dados) %in%
                  c("SEXO", "COR", "RENDA MÉDIA FAMILIAR",
                    "ÁREA DE INTERESSE")] <- c("sexo", "cor", "rfm",
                                                       "areaint")
## Removendo valores indesejados
# Remover os valores "Todas" e "Nenhuma"
dados <- dados[!(dados$areaint %in% c(5, 6)),]

# Remover os valores 7 e 9 da coluna rfm
dados <- dados[!(dados$rfm %in% c(7, 9, 99)),]

# Remover "Outras respostas", "Outras/NR"
dados <- dados[!(dados$cor %in% c(6, 9)),]

## Renomeando as variáveis
dados$sexo <- ifelse(dados$sexo == 1, "Masculino", "Feminino")
dados$cor <- factor(dados$cor,
                    levels = c(1, 2, 3, 4, 5),
                    labels = c("Branca", "Preta", "Parda", "Amarela",
                               "Indígena"))

dados$areaint <- factor(dados$areaint,
                     levels = c(1, 2, 3, 4),
                     labels = c("Antropologia", "Ciência Política", "Sociologia",
                                "Ainda não sabe"))

## Formatando os dados como integer e factor

dados$sexo <- as.factor(dados$sexo)
dados$cor <- as.factor(dados$cor)
dados$areaint <- as.factor(dados$areaint)

dados$rfm <- as.integer(round(dados$rfm))

### Garantindo que o banco de dados não tenha valores N'As
dados <- dados[complete.cases(dados[, c("sexo", "cor", "areaint", "rfm")]), ]

## 3. Ausência de multicolinearidade
psych::pairs.panels(dados[, c(3, 4, 5, 6)])

m <- lm(as.numeric(areaint) ~ sexo + cor + rfm, data = dados)
car::vif(m)


## 4. Independência de alternativas irrelevantes (teste Hausman-McFadden)
# Críticas válidas aos testes que checam esse pressuposto: https://statisticalhorizons.com/iia
install.packages("mlogit")
library(mlogit)

modiia <- mlogit::mlogit(areaint ~ 1 | sexo + cor + rfm,
                         data = dados, shape = "wide",
                         reflevel = "Antropologia")

modiia2 <- mlogit::mlogit(areaint ~ 1 | sexo + cor + rfm,
                         data = dados, shape = "wide",
                         reflevel = "Antropologia",
                         alt.subset = c("Antropologia", "Ciência Política", "Ainda não sabe"))

modiia3 <- mlogit::mlogit(areaint ~ 1 | sexo + cor + rfm,
                          data = dados, shape = "wide",
                          reflevel = "Antropologia",
                          alt.subset = c("Antropologia", "Ciência Política", "Sociologia"))

modiia4 <- mlogit::mlogit(areaint ~ 1 | sexo + cor + rfm,
                          data = dados, shape = "wide",
                          reflevel = "Antropologia",
                          alt.subset = c("Antropologia", "Sociologia", "Ainda não sabe"))



mlogit::hmftest(modiia, modiia2)
mlogit::hmftest(modiia, modiia3)
mlogit::hmftest(modiia, modiia4)

# Passo 6. Construção do modelo e interpretação dos resultados

## Construção do modelo e do modelo nulo (usando o pacote nnet):
mod <- multinom(areaint ~ sexo + cor + rfm, data = dados, model = TRUE)
mod0 <- multinom(areaint ~ 1, data = dados, model = TRUE)

# Ajuste do modelo
anova(mod, mod0)
DescTools::PseudoR2(mod, which = "Nagelkerke")

# Overall effects
car::Anova(mod, type = "II", test = "Wald")

# Efeitos específicos
summary(mod)

## Obtenção dos valores de p - por Wald (pacote lmtest)
lmtest::coeftest(mod)

## Obtenção das razões de chance com IC 95% (usando log-likelihood)
exp(coef(mod))
exp(confint(mod))

## Tabela completa (pacote gtsummary)
gtsummary::tbl_regression(mod, exponentiate = FALSE)
gtsummary::tbl_regression(mod, exponentiate = TRUE)

# Filtrar os dados para remover a categoria "Ainda não sabe" de AreaInteresse
dados_prev_filtered <- dados_prev %>%
  filter(AreaInteresse != "Ainda não sabe") %>%
  mutate(AreaInteresse = fct_recode(AreaInteresse,
                                    "C. Política" = "Ciência Política"))

# Arredondar as probabilidades para números inteiros
dados_prev_filtered$Probabilidade <- round(dados_prev_filtered$Probabilidade * 100)

# Criando o gráfico
ggplot(dados_prev_filtered, aes(x = rfm, y = Probabilidade, color = sexo)) +
  geom_point(alpha = 0.6, size = 1.5) +  # Pontos com transparência
  geom_smooth(method = "loess", size = 0.8, linetype = "solid", se = TRUE) +  # Linha de suavização com intervalo de confiança
  labs(
    x = "Renda Familiar, em salários mínimos (R$)",
    y = "Probabilidade (%)",
    color = "Sexo"
  ) +
  scale_y_continuous(labels = scales::number_format(decimal.mark = ",", accuracy = 1)) +  # Mostrar como inteiros
  scale_x_continuous(labels = scales::number_format(decimal.mark = ",", accuracy = 1)) +  # Ajustar os rótulos do eixo x
  facet_grid(AreaInteresse ~ ., scales = "fixed") +  # Facetas com escalas fixas no eixo y
  theme_bw() +
  theme(
    text = element_text(size = 12),
    axis.title = element_text(size = 14),
    axis.text = element_text(size = 12),
    legend.title = element_text(size = 12),
    legend.text = element_text(size = 10),
    strip.text = element_text(size = 12),
    panel.grid.minor = element_blank()
  ) +
  guides(color = guide_legend(override.aes = list(fill = NA)))