Carregamento de Pacotes e Dados

library(haven)
library(dplyr)
library(ggplot2)
library(broom)
library(knitr)
library(purrr)
library(tidyr)
library(gridExtra)
library(corrplot)
library(car)
# Carregamento da base (assumindo que já está filtrada para Q1 do ESCS)
dados <- read_sav("../2-base-de-analise/PISA2022_LATAM_resiliencia.sav")

# Seleção de variáveis necessárias
variaveis <- c("CNT", "resil_math", "resil_read", "resil_scie", "ICTREG", "ICTEFFIC", "ST004D01T", "ESCS")
variaveis_existentes <- variaveis[variaveis %in% names(dados)]
dados_trabalho <- dados %>% select(all_of(variaveis_existentes))

# Checagem opcional: garantir Q1 de ESCS, caso a base não esteja filtrada
if("ESCS" %in% names(dados_trabalho)){
  # Mantém apenas o primeiro quartil dentro da própria amostra (se já estiver filtrada, nada muda)
  q1 <- quantile(dados_trabalho$ESCS, probs = 0.25, na.rm = TRUE)
  n_before <- nrow(dados_trabalho)
  dados_trabalho <- dados_trabalho %>% filter(!is.na(ESCS) & ESCS <= q1)
  cat("Observações mantidas após garantir Q1 do ESCS:", nrow(dados_trabalho), "(de", n_before, ")
")
}
## Observações mantidas após garantir Q1 do ESCS: 4788 (de 19143 )
# Limpeza essencial
dados_limpos <- dados_trabalho %>%
  filter(!is.na(CNT), !is.na(ICTREG), !is.na(ICTEFFIC),
         !is.na(resil_math), !is.na(resil_read), !is.na(resil_scie)) %>%
  mutate(
    genero = if("ST004D01T" %in% names(.)) ifelse(ST004D01T == 1, "Feminino", "Masculino") else NA_character_,
    ICTREG_std = as.numeric(scale(ICTREG)),
    ICTEFFIC_std = as.numeric(scale(ICTEFFIC))
  ) %>%
  mutate(genero = if("ST004D01T" %in% names(.)) factor(genero) else NULL)

str(dados_limpos)
## tibble [742 × 11] (S3: tbl_df/tbl/data.frame)
##  $ CNT         : chr+lbl [1:742] ARG, ARG, ARG, ARG, ARG, ARG, ARG, ARG, ARG, ARG, ARG,...
##    ..@ label        : chr "Country code 3-character"
##    ..@ format.spss  : chr "A3"
##    ..@ display_width: int 3
##    ..@ labels       : Named chr [1:81] "MDA" "THA" "BRA" "FRA" ...
##    .. ..- attr(*, "names")= chr [1:81] "Republic of Moldova" "Thailand" "Brazil" "France" ...
##  $ resil_math  : num [1:742] 0 0 0 0 0 0 0 0 0 0 ...
##   ..- attr(*, "format.spss")= chr "F8.0"
##  $ resil_read  : num [1:742] 0 0 0 0 0 0 0 0 0 0 ...
##   ..- attr(*, "format.spss")= chr "F8.0"
##  $ resil_scie  : num [1:742] 0 0 0 1 0 0 0 0 0 0 ...
##   ..- attr(*, "format.spss")= chr "F8.0"
##  $ ICTREG      : dbl+lbl [1:742]  0.1612,  0.6260,  3.2020,  0.3759,  1.1799,  0.1169, ...
##    ..@ label      : chr "Views of regulated ICT use in school (WLE)"
##    ..@ format.spss: chr "F7.4"
##    ..@ labels     : Named num [1:4] 95 97 98 99
##    .. ..- attr(*, "names")= chr [1:4] "Valid Skip" "Not Applicable" "Invalid" "No Response"
##  $ ICTEFFIC    : dbl+lbl [1:742]  0.0490, -1.1925, -0.5172, -1.2300,  0.2302,  0.6282, ...
##    ..@ label      : chr "Self-efficacy in digital competencies (WLE)"
##    ..@ format.spss: chr "F7.4"
##    ..@ labels     : Named num [1:4] 95 97 98 99
##    .. ..- attr(*, "names")= chr [1:4] "Valid Skip" "Not Applicable" "Invalid" "No Response"
##  $ ST004D01T   : dbl+lbl [1:742] 2, 2, 2, 1, 2, 2, 2, 1, 1, 1, 2, 1, 1, 1, 1, 2, 2, 2, ...
##    ..@ label      : chr "Student (Standardized) Gender"
##    ..@ format.spss: chr "F1.0"
##    ..@ labels     : Named num [1:6] 1 2 5 7 8 9
##    .. ..- attr(*, "names")= chr [1:6] "Female" "Male" "Valid Skip" "Not Applicable" ...
##  $ ESCS        : dbl+lbl [1:742] -2.90, -2.98, -2.89, -3.09, -3.12, -3.22, -3.16, -3.20...
##    ..@ label      : chr "Index of economic, social and cultural status"
##    ..@ format.spss: chr "F7.4"
##    ..@ labels     : Named num [1:4] 95 97 98 99
##    .. ..- attr(*, "names")= chr [1:4] "Valid Skip" "Not Applicable" "Invalid" "No Response"
##  $ genero      : Factor w/ 2 levels "Feminino","Masculino": 2 2 2 1 2 2 2 1 1 1 ...
##  $ ICTREG_std  : num [1:742] -0.0266 0.3504 2.4396 0.1475 0.7996 ...
##  $ ICTEFFIC_std: num [1:742] 0.707 -0.446 0.181 -0.481 0.875 ...
cat("Amostra final:", nrow(dados_limpos), "observações; Países:", length(unique(dados_limpos$CNT)), "
")
## Amostra final: 742 observações; Países: 6

Exploração Rápida

# Estatísticas de resiliência por país no Q1
desc_pais <- dados_limpos %>%
  group_by(CNT) %>%
  summarise(
    n = n(),
    resil_math_prop = mean(resil_math),
    resil_read_prop = mean(resil_read),
    resil_scie_prop = mean(resil_scie),
    ICTREG_media = mean(ICTREG),
    ICTEFFIC_media = mean(ICTEFFIC),
    .groups = 'drop'
  )

kable(desc_pais, digits = 3, caption = "Descritivas por País (Q1 ESCS)")
Descritivas por País (Q1 ESCS)
CNT n resil_math_prop resil_read_prop resil_scie_prop ICTREG_media ICTEFFIC_media
ARG 147 0.041 0.061 0.061 0.191 -0.899
BRA 282 0.089 0.121 0.096 0.175 -0.652
CHL 24 0.000 0.042 0.042 -0.169 -0.516
DOM 82 0.134 0.085 0.073 0.019 -0.530
PAN 138 0.051 0.051 0.036 0.511 -0.890
URY 69 0.072 0.058 0.087 -0.020 -0.489 
# Correlação entre preditores (sem CNT)
vars_num <- dados_limpos %>% select(where(is.numeric)) %>% select(-starts_with("resil_")) %>% na.omit()
if(ncol(vars_num) > 1){
  corrplot(cor(vars_num), method = "color", type = "upper", order = "hclust", tl.cex = 0.8)
}

Estratégia de Modelagem

Serão estimados, para cada disciplina, três modelos: - Modelo 1: ICTREG_std + ICTEFFIC_std - Modelo 2: Modelo 1 + gênero - Modelo 3: Modelo 2 + efeitos fixos de país (CNT)

fit_series <- function(df, y){
  # Constrói fórmulas
  f1 <- as.formula(paste(y, "~ ICTREG_std + ICTEFFIC_std"))
  f2 <- as.formula(paste(y, "~ ICTREG_std + ICTEFFIC_std + genero"))
  f3 <- as.formula(paste(y, "~ ICTREG_std + ICTEFFIC_std + genero + CNT"))

  # Ajuste com tryCatch para robustez
  safe_glm <- function(form){
    tryCatch(glm(form, data = df, family = binomial(link = "logit")), error = function(e) NULL)
  }

  m1 <- safe_glm(f1)
  m2 <- if("genero" %in% names(df)) safe_glm(f2) else NULL
  m3 <- if(!is.null(m2)) safe_glm(f3) else if("genero" %in% names(df)) NULL else safe_glm(as.formula(paste(y, "~ ICTREG_std + ICTEFFIC_std + CNT")))

  list(m1 = m1, m2 = m2, m3 = m3)
}

extract_table <- function(model, nome){
  if(is.null(model)) return(NULL)
  out <- tidy(model)
  # filtra apenas termos de interesse
  out <- out %>% filter(grepl("ICTREG_std|ICTEFFIC_std|genero", term)) %>%
    mutate(
      modelo = nome,
      OR = exp(estimate),
      OR_lower = exp(estimate - 1.96*std.error),
      OR_upper = exp(estimate + 1.96*std.error),
      sig = ifelse(p.value < 0.05, "Sim", "Não")
    ) %>%
    select(modelo, term, estimate, std.error, OR, OR_lower, OR_upper, statistic, p.value, sig)
  out
}

Matemática

mods_math <- fit_series(dados_limpos, "resil_math")
res_math <- bind_rows(
  extract_table(mods_math$m1, "Modelo 1: ICT"),
  extract_table(mods_math$m2, "Modelo 2: ICT + Gênero"),
  extract_table(mods_math$m3, "Modelo 3: ICT + Gênero + País")
)

kable(res_math, digits = 4, caption = "Regressão Logística Múltipla - Matemática (Q1 ESCS)")
Regressão Logística Múltipla - Matemática (Q1 ESCS)
modelo term estimate std.error OR OR_lower OR_upper statistic p.value sig
Modelo 1: ICT ICTREG_std 0.2132 0.1575 1.2377 0.9090 1.6851 1.3542 0.1757 Não
Modelo 1: ICT ICTEFFIC_std 0.3645 0.1311 1.4398 1.1135 1.8616 2.7798 0.0054 Sim
Modelo 2: ICT + Gênero ICTREG_std 0.2067 0.1578 1.2296 0.9024 1.6753 1.3094 0.1904 Não
Modelo 2: ICT + Gênero ICTEFFIC_std 0.3645 0.1311 1.4399 1.1135 1.8618 2.7801 0.0054 Sim
Modelo 2: ICT + Gênero generoMasculino 0.1220 0.2871 1.1298 0.6435 1.9834 0.4249 0.6709 Não
Modelo 3: ICT + Gênero + País ICTREG_std 0.2406 0.1648 1.2720 0.9209 1.7569 1.4600 0.1443 Não
Modelo 3: ICT + Gênero + País ICTEFFIC_std 0.3231 0.1331 1.3815 1.0643 1.7931 2.4283 0.0152 Sim
Modelo 3: ICT + Gênero + País generoMasculino 0.1288 0.2903 1.1375 0.6439 2.0092 0.4437 0.6572 Não 
comp_math <- data.frame(
  Modelo = c("M1", "M2", "M3"),
  AIC = c(if(!is.null(mods_math$m1)) AIC(mods_math$m1) else NA,
          if(!is.null(mods_math$m2)) AIC(mods_math$m2) else NA,
          if(!is.null(mods_math$m3)) AIC(mods_math$m3) else NA),
  BIC = c(if(!is.null(mods_math$m1)) BIC(mods_math$m1) else NA,
          if(!is.null(mods_math$m2)) BIC(mods_math$m2) else NA,
          if(!is.null(mods_math$m3)) BIC(mods_math$m3) else NA)
)

kable(comp_math, digits = 2, caption = "Comparação de Modelos - Matemática")
Comparação de Modelos - Matemática
Modelo AIC BIC
M1 382.90 396.73
M2 384.72 403.16
M3 383.63 425.11

Leitura

mods_read <- fit_series(dados_limpos, "resil_read")
res_read <- bind_rows(
  extract_table(mods_read$m1, "Modelo 1: ICT"),
  extract_table(mods_read$m2, "Modelo 2: ICT + Gênero"),
  extract_table(mods_read$m3, "Modelo 3: ICT + Gênero + País")
)

kable(res_read, digits = 4, caption = "Regressão Logística Múltipla - Leitura (Q1 ESCS)")
Regressão Logística Múltipla - Leitura (Q1 ESCS)
modelo term estimate std.error OR OR_lower OR_upper statistic p.value sig
Modelo 1: ICT ICTREG_std 0.0361 0.1396 1.0368 0.7886 1.3631 0.2589 0.7957 Não
Modelo 1: ICT ICTEFFIC_std 0.4104 0.1219 1.5074 1.1870 1.9142 3.3658 0.0008 Sim
Modelo 2: ICT + Gênero ICTREG_std 0.0417 0.1405 1.0425 0.7916 1.3730 0.2965 0.7669 Não
Modelo 2: ICT + Gênero ICTEFFIC_std 0.4107 0.1220 1.5078 1.1872 1.9149 3.3673 0.0008 Sim
Modelo 2: ICT + Gênero generoMasculino -0.1175 0.2738 0.8891 0.5199 1.5205 -0.4292 0.6677 Não
Modelo 3: ICT + Gênero + País ICTREG_std 0.0695 0.1467 1.0719 0.8040 1.4291 0.4735 0.6358 Não
Modelo 3: ICT + Gênero + País ICTEFFIC_std 0.3936 0.1238 1.4824 1.1629 1.8896 3.1790 0.0015 Sim
Modelo 3: ICT + Gênero + País generoMasculino -0.1185 0.2764 0.8882 0.5167 1.5268 -0.4289 0.6680 Não 
comp_read <- data.frame(
  Modelo = c("M1", "M2", "M3"),
  AIC = c(if(!is.null(mods_read$m1)) AIC(mods_read$m1) else NA,
          if(!is.null(mods_read$m2)) AIC(mods_read$m2) else NA,
          if(!is.null(mods_read$m3)) AIC(mods_read$m3) else NA),
  BIC = c(if(!is.null(mods_read$m1)) BIC(mods_read$m1) else NA,
          if(!is.null(mods_read$m2)) BIC(mods_read$m2) else NA,
          if(!is.null(mods_read$m3)) BIC(mods_read$m3) else NA)
)

kable(comp_read, digits = 2, caption = "Comparação de Modelos - Leitura")
Comparação de Modelos - Leitura
Modelo AIC BIC
M1 421.37 435.20
M2 423.19 441.63
M3 424.99 466.47

Ciências

mods_scie <- fit_series(dados_limpos, "resil_scie")
res_scie <- bind_rows(
  extract_table(mods_scie$m1, "Modelo 1: ICT"),
  extract_table(mods_scie$m2, "Modelo 2: ICT + Gênero"),
  extract_table(mods_scie$m3, "Modelo 3: ICT + Gênero + País")
)

kable(res_scie, digits = 4, caption = "Regressão Logística Múltipla - Ciências (Q1 ESCS)")
Regressão Logística Múltipla - Ciências (Q1 ESCS)
modelo term estimate std.error OR OR_lower OR_upper statistic p.value sig
Modelo 1: ICT ICTREG_std 0.0728 0.1510 1.0755 0.8000 1.4458 0.4819 0.6299 Não
Modelo 1: ICT ICTEFFIC_std 0.4523 0.1282 1.5720 1.2227 2.0210 3.5282 0.0004 Sim
Modelo 2: ICT + Gênero ICTREG_std 0.0657 0.1512 1.0679 0.7939 1.4364 0.4344 0.6640 Não
Modelo 2: ICT + Gênero ICTEFFIC_std 0.4523 0.1282 1.5720 1.2226 2.0212 3.5275 0.0004 Sim
Modelo 2: ICT + Gênero generoMasculino 0.1462 0.2877 1.1575 0.6586 2.0342 0.5083 0.6112 Não
Modelo 3: ICT + Gênero + País ICTREG_std 0.1038 0.1586 1.1093 0.8129 1.5139 0.6540 0.5131 Não
Modelo 3: ICT + Gênero + País ICTEFFIC_std 0.4249 0.1298 1.5294 1.1858 1.9727 3.2724 0.0011 Sim
Modelo 3: ICT + Gênero + País generoMasculino 0.1456 0.2898 1.1567 0.6555 2.0412 0.5023 0.6154 Não 
comp_scie <- data.frame(
  Modelo = c("M1", "M2", "M3"),
  AIC = c(if(!is.null(mods_scie$m1)) AIC(mods_scie$m1) else NA,
          if(!is.null(mods_scie$m2)) AIC(mods_scie$m2) else NA,
          if(!is.null(mods_scie$m3)) AIC(mods_scie$m3) else NA),
  BIC = c(if(!is.null(mods_scie$m1)) BIC(mods_scie$m1) else NA,
          if(!is.null(mods_scie$m2)) BIC(mods_scie$m2) else NA,
          if(!is.null(mods_scie$m3)) BIC(mods_scie$m3) else NA)
)

kable(comp_scie, digits = 2, caption = "Comparação de Modelos - Ciências")
Comparação de Modelos - Ciências
Modelo AIC BIC
M1 380.60 394.42
M2 382.34 400.78
M3 387.23 428.71

Robustez e Diagnósticos

# VIF no melhor modelo com variáveis ICT + gênero (sem ESCS)
if(!is.null(mods_math$m2)){
  cat("VIF - Modelo 2 (Matemática):
")
  print(vif(mods_math$m2))
}
## VIF - Modelo 2 (Matemática):
##   ICTREG_std ICTEFFIC_std       genero 
##     1.013990     1.004369     1.009601
# Resíduos do modelo com efeitos fixos por país (se existir)
modelo_final <- mods_math$m3
if(!is.null(modelo_final)){
  residuos <- residuals(modelo_final, type = "deviance")
  fit <- fitted(modelo_final)
  ggplot(data.frame(fit = fit, res = residuos), aes(fit, res)) +
    geom_point(alpha = 0.4) +
    geom_smooth(method = "loess", se = FALSE, color = "red") +
    geom_hline(yintercept = 0, linetype = "dashed") +
    theme_minimal() +
    labs(title = "Resíduos vs Ajustados - Matemática (M3)", x = "Ajustados", y = "Resíduos (Deviance)")
}

Visualização Comparativa

all_res <- bind_rows(
  res_math %>% mutate(disciplina = "Matemática"),
  res_read %>% mutate(disciplina = "Leitura"),
  res_scie %>% mutate(disciplina = "Ciências")
)

p <- ggplot(all_res, aes(x = modelo, y = OR, color = term)) +
  geom_point(position = position_dodge(width = 0.3), size = 3) +
  geom_errorbar(aes(ymin = OR_lower, ymax = OR_upper), position = position_dodge(width = 0.3), width = 0.2) +
  geom_hline(yintercept = 1, linetype = "dashed", color = "red") +
  facet_wrap(~disciplina, scales = "free_y") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title = "Efeitos de ICTREG e ICTEFFIC (Q1 ESCS)", x = "Modelo", y = "OR (IC 95%)", color = "Variável") +
  scale_color_manual(values = c("ICTREG_std" = "blue", "ICTEFFIC_std" = "green"), labels = c("ICTREG", "ICTEFFIC"))

print(p)

# Tabela resumo do modelo final (M3)
resumo_final <- all_res %>% filter(modelo == "Modelo 3: ICT + Gênero + País") %>%
  arrange(disciplina, term)

kable(resumo_final, digits = 4, caption = "Resumo dos Efeitos (Modelo Final, M3)")
Resumo dos Efeitos (Modelo Final, M3)
modelo term estimate std.error OR OR_lower OR_upper statistic p.value sig disciplina
Modelo 3: ICT + Gênero + País ICTEFFIC_std 0.4249 0.1298 1.5294 1.1858 1.9727 3.2724 0.0011 Sim Ciências
Modelo 3: ICT + Gênero + País ICTREG_std 0.1038 0.1586 1.1093 0.8129 1.5139 0.6540 0.5131 Não Ciências
Modelo 3: ICT + Gênero + País generoMasculino 0.1456 0.2898 1.1567 0.6555 2.0412 0.5023 0.6154 Não Ciências
Modelo 3: ICT + Gênero + País ICTEFFIC_std 0.3936 0.1238 1.4824 1.1629 1.8896 3.1790 0.0015 Sim Leitura
Modelo 3: ICT + Gênero + País ICTREG_std 0.0695 0.1467 1.0719 0.8040 1.4291 0.4735 0.6358 Não Leitura
Modelo 3: ICT + Gênero + País generoMasculino -0.1185 0.2764 0.8882 0.5167 1.5268 -0.4289 0.6680 Não Leitura
Modelo 3: ICT + Gênero + País ICTEFFIC_std 0.3231 0.1331 1.3815 1.0643 1.7931 2.4283 0.0152 Sim Matemática
Modelo 3: ICT + Gênero + País ICTREG_std 0.2406 0.1648 1.2720 0.9209 1.7569 1.4600 0.1443 Não Matemática
Modelo 3: ICT + Gênero + País generoMasculino 0.1288 0.2903 1.1375 0.6439 2.0092 0.4437 0.6572 Não Matemática

Interpretação

Resultados referem-se exclusivamente ao subconjunto de estudantes no 1º quartil do ESCS (vulnerabilidade social). O efeito de ICTREG e ICTEFFIC é interpretado como variação por desvio-padrão (padronização z-score). Modelos com efeitos fixos de país controlam por diferenças sistemáticas entre países. P-values referem-se a testes Wald e ORs incluem intervalos de 95%.