1 Problema

O objetivo da modelagem é prever o resultado no exame normalizado (normexam) de acordo com as características dos alunos e das escolas.

Quais os tipos de escola (gênero misto, só de homens ou só de mulheres) é mais eficiente na formação de seus alunos de acordo com o resultado do exame normalizado?

2 Planejamento da Solução

Carregar as bases de dados individuais.
Fazer uma limpeza dos dados de cada uma das bases de dados incluindo:

Mudança de tipo de variável;
Identificação e tratamento de duplicidade;
Mudar formato da base de “longo” para “largo”, ou seja, transformar em variáveis as observações nas linhas, mantendo alguma variável de identificação.

Fazer uma análise exploratória dos dados para verificar o comportamento dos dados e buscar insights dos mesmos.
Divisão da base em treino e teste.
Ajustar um modelo de regressão inicial (baseline) com o qual faremos comparações ao criar novos modelos.
Fazer predições com a base de teste.
Extrair Métricas de desempenho do modelo e comparar com o modelo baseline.
Salvamento do modelo final para ser usado em produção.

3 Pacotes

library(readr)
library(tidyverse)
library(reshape2)
library(skimr)
library(broom)
library(tidymodels)
library(gglm)
library(broom.mixed) 
library(dotwhisker) 
library(yardstick)
library(vip)
library(ggpubr)
library(modeldata)

4 Base de dados

cat_school_data <- read.csv("base_de_dados/cat_school_data.csv", sep = ';')

cat_student_data <- read.csv("base_de_dados/cat_student_data.csv", sep = ';')

num_school_data <- read.csv("base_de_dados/num_school_data.csv", sep = ';')

num_student_data <- read.csv("base_de_dados/num_student_data.csv", sep = ';')

5 Limpeza de dados

# CAT SCHOOL ------------------------
cat_school_df <- cat_school_data %>% 
  spread(key = variable, value = value)

# CAT STUDENT --------------------------------------
# Mantendo somente linhas distintas
cat_student_dist <- cat_student_data %>% 
  distinct()

# Checando se ainda temos duplicidade:
cat_student_dist %>%
  dplyr::group_by(school, student, variable) %>%
  dplyr::summarise(n = dplyr::n(), .groups = "drop") %>%
  dplyr::filter(n > 1L)

## # A tibble: 2 × 4
##   school student variable     n
##    <int>   <int> <chr>    <int>
## 1     43      86 intake       2
## 2     43      86 sex          2

# Removendo duplicidade remanescente 
cat_student_dist_aux <-  cat_student_dist %>% 
  filter(school != 43 & student != 86)

# Salvando em formato wide
cat_student_df <- cat_student_dist_aux %>% 
  spread(key = variable, value = value)

# NUM SCHOOL --------------------------------------
num_school_df <-  num_school_data %>% 
  spread(key = variable, value = value) %>% 
  mutate(schavg = as.numeric(gsub(",", ".", schavg)))

# NUM STUDENT --------------------------------------
# checando duplicidade
num_student_data %>%
  dplyr::group_by(school, student, variable) %>%
  dplyr::summarise(n = dplyr::n(), .groups = "drop") %>%
  dplyr::filter(n > 1L)

## # A tibble: 34 × 4
##    school student variable     n
##     <int>   <int> <chr>    <int>
##  1     14      NA normexam     2
##  2     14      NA standLRT     2
##  3     17      NA normexam     2
##  4     17      NA standLRT     2
##  5     22      NA normexam     2
##  6     22      NA standLRT     2
##  7     35      NA normexam     3
##  8     35      NA standLRT     3
##  9     40      NA normexam     2
## 10     40      NA standLRT     2
## # … with 24 more rows

# removendo NA de `student` e mantendo linhas distintas
num_student_data_nona <- num_student_data %>% 
  filter(!is.na(student)) %>% 
  distinct()

# duplicidades remanescentes 
duplicate <- num_student_data_nona %>%
  dplyr::group_by(school, student, variable) %>%
  dplyr::summarise(n = dplyr::n(), .groups = "drop") %>%
  dplyr::filter(n > 1L) 

# criando uma sub-base somente com os valores duplicados 
duplicate_join <- semi_join(num_student_data_nona,
          duplicate) 

# excluindo os duplicados
num_student_data_nona <- anti_join(num_student_data_nona,
          duplicate_join)

# Salvando em formato wide
num_student_df <- num_student_data_nona %>% 
  spread(key = variable, value = value) %>% 
  mutate(normexam = as.numeric(gsub(",", ".", normexam)),
         standLRT = as.numeric(gsub(",", ".", standLRT)))

# JUNTANDO AS BASES DE DADOS -------------------------------
school_df_raw <- num_school_df %>% 
  left_join(cat_school_df)

student_df_raw <- num_student_df %>% 
  left_join(cat_student_df, by = c("school", "student"))

final_data_raw <- student_df_raw %>% 
  left_join(school_df_raw, by = "school") %>% 
  mutate(school = factor(school),
         student = factor(student))

# salvando em csv
write.csv(final_data_raw, "sch_std_data.csv", row.names = FALSE)

6 Análise Exploratória dos Dados

Aqui, carregamos a base de dados final resultante da limpeza feita anteriormente:

data_raw <- read_csv("sch_std_data.csv")
glimpse(data_raw)

## Rows: 4,006
## Columns: 10
## $ school   <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ student  <dbl> 1, 4, 6, 7, 13, 14, 16, 17, 19, 22, 27, 28, 31, 36, 43, 56, 5…
## $ normexam <dbl> 1.5061852, -0.5551120, -1.3353150, -0.5551120, -0.1976110, -0…
## $ standLRT <dbl> 0.7843622, -1.0339700, -0.9513180, -2.3563930, -0.2901070, 0.…
## $ intake   <chr> "bottom 25%", "mid 50%", "mid 50%", "top 25%", "mid 50%", "mi…
## $ sex      <chr> "F", "F", "M", "M", "M", "M", "F", "F", NA, "M", "M", "M", "M…
## $ vr       <chr> "mid 50%", "mid 50%", "mid 50%", "mid 50%", "mid 50%", "mid 5…
## $ schavg   <dbl> 0.1661752, 0.1661752, 0.1661752, 0.1661752, 0.1661752, 0.1661…
## $ schgend  <chr> "mixed", "mixed", "mixed", "mixed", "mixed", "mixed", "mixed"…
## $ type     <chr> "Mxd", "Mxd", "Mxd", "Mxd", "Mxd", "Mxd", "Mxd", "Mxd", "Mxd"…

skim(data_raw)

Data summary
Name	data_raw
Number of rows	4006
Number of columns	10
_______________________
Column type frequency:
character	5
numeric	5
________________________
Group variables	None

Variable type: character

skim_variable	n_missing	complete_rate	min	max	n_unique
intake	147	0.96	7	10	3
sex	173	0.96	1	1	2
vr	134	0.97	7	10	3
schgend	0	1.00	4	5	3
type	0	1.00	3	4	2

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
school	0	1.00	30.90	18.93	1.00	14.00	29.00	47.00	65.00	▆▇▅▆▅
student	0	1.00	134.40	174.61	1.00	39.00	78.00	138.00	913.00	▇▁▁▁▁
normexam	70	0.98	0.00	1.00	-3.67	-0.70	0.00	0.68	3.67	▁▃▇▃▁
standLRT	29	0.99	0.01	0.99	-2.93	-0.62	0.04	0.62	3.02	▁▃▇▃▁
schavg	0	1.00	0.00	0.32	-0.76	-0.15	-0.02	0.27	0.64	▂▃▇▆▃

Num primeiro momento, vamos excluir os NAs. Num próximo ciclo, podemos fazer um imputação usando mediana ou moda (ou outro tipo de imputação), a depender do tipo de variável.

data <- na.omit(final_data_raw)
GGally::ggpairs(data %>% 
                  select(-c(student, school))
                )

Nossa variável dependente normexam tem ditribuição normal e apresenta correlação significativa com as variáveis standLRT (r = 0.588, correlação moderada) e schavg (r = 0.285, correlação fraca).
Quanto ao tipo e gênero da (type e schgend) não parece haver muita diferança, na média, no normexam.
Quanto às variáveis dos alunos, não parece haver diferença, na média, do normexam entre os sexos.
Para a variável intake: a classe bottom25% apresenta, na média, um maior normexam .
Para a variável vr: a classe top25% apresenta, na média, um maior normexam.

6.1 Testando hipóteses

O resultado do exame normalizado depende do tipo de escola?
- Sim, há uma diferença estatística significativa entre os exames normalizados para escolas mistas e escolas de um único gênero. As escolas de um único gênero tem, em média 0.113 na nota do exame normalizado, enquanto que as de gênero misto têm -0.099.

library(rstatix)
data %>% 
  group_by(type) %>% 
  get_summary_stats(normexam, type = "full")

## # A tibble: 2 × 14
##   type  variable     n   min   max median     q1    q3   iqr   mad   mean    sd
##   <chr> <chr>    <dbl> <dbl> <dbl>  <dbl>  <dbl> <dbl> <dbl> <dbl>  <dbl> <dbl>
## 1 Mxd   normexam  1915 -3.67  2.92 -0.129 -0.776 0.578  1.35 0.998 -0.099  0.98
## 2 Sngl  normexam  1717 -3.06  3.67  0.134 -0.555 0.747  1.30 1.02   0.113  1.01
## # … with 2 more variables: se <dbl>, ci <dbl>

# Teste
data %>%
  t_test(normexam ~ type) %>% 
  add_significance()

## # A tibble: 1 × 9
##   .y.      group1 group2    n1    n2 statistic    df        p p.signif
## * <chr>    <chr>  <chr>  <int> <int>     <dbl> <dbl>    <dbl> <chr>   
## 1 normexam Mxd    Sngl    1915  1717     -6.40 3560. 1.76e-10 ****

# Gráfico
data %>% 
  ggplot(aes(x = type, y = normexam)) +
  geom_boxplot() +
  stat_summary(fun.y = mean, shape = 15, color = "blue")

O resultado do exame normalizado depende do gênero da escola?
- Escolas somente de garotas, possui a maior média no exame normalizado, 0.144, enquanto que a média para escola só de garotos é de 0.029, e para escolas mistas, de -0.099.
- Há uma diferença significativa no exame normalizado para o gênero da escola (ANOVA, p < 0.0001), ou seja, pelo menos um gênero de escola se sobresai.
- Na comparação dois a dois, não foi encontrada diferença significativa entre escolas só de garotos ou só de garotas. Mas ambas as escolas de gênero único apresentaram diferença significativa quando comparadas à escolas mistas.

library(rstatix)
data %>% 
  group_by(schgend) %>% 
  get_summary_stats(normexam, type = "full")

## # A tibble: 3 × 14
##   schgend variable     n   min   max median     q1    q3   iqr   mad   mean
##   <chr>   <chr>    <dbl> <dbl> <dbl>  <dbl>  <dbl> <dbl> <dbl> <dbl>  <dbl>
## 1 boys    normexam   467 -2.75  2.53  0.074 -0.738 0.822  1.56 1.15   0.029
## 2 girls   normexam  1250 -3.06  3.67  0.134 -0.493 0.747  1.24 0.929  0.144
## 3 mixed   normexam  1915 -3.67  2.92 -0.129 -0.776 0.578  1.35 0.998 -0.099
## # … with 3 more variables: sd <dbl>, se <dbl>, ci <dbl>

# Teste
data %>%
  anova_test(normexam ~ schgend) %>% 
  add_significance()

## ANOVA Table (type II tests)
## 
##    Effect DFn  DFd      F       p p<.05   ges p.signif
## 1 schgend   2 3629 22.834 1.4e-10     * 0.012     ****

# Post-hoc
data %>% 
  tukey_hsd(normexam ~ schgend)

## # A tibble: 3 × 9
##   term    group1 group2 null.value estimate conf.low conf.high         p.adj
## * <chr>   <chr>  <chr>       <dbl>    <dbl>    <dbl>     <dbl>         <dbl>
## 1 schgend boys   girls           0    0.115  -0.0117   0.241   0.0843       
## 2 schgend boys   mixed           0   -0.128  -0.249   -0.00808 0.0332       
## 3 schgend girls  mixed           0   -0.243  -0.328   -0.158   0.00000000688
## # … with 1 more variable: p.adj.signif <chr>

# Gráfico
data %>% 
  ggplot(aes(x = schgend, y = normexam)) +
  geom_boxplot() +
  stat_summary(fun.y = mean, shape = 15, color = "blue")

6.2 Alguns gráficos exploratórios adicionais

data %>%
  select(where(is.numeric)) %>%
  cor(use = "pairwise.complete.obs") %>%
  corrplot::corrplot()

g1 <- ggplot(data,
       aes(x = standLRT, 
           y = normexam, 
           group = type, 
           col = type)) + 
  geom_point(alpha = 0.2) + 
  geom_smooth(method = lm, se = FALSE) +
  scale_color_viridis_d(option = "plasma", end = .7)

g2 <- ggplot(data,
       aes(x = standLRT, 
           y = normexam, 
           group = sex, 
           col = sex)) + 
  geom_point(alpha = 0.2) + 
  geom_smooth(method = lm, se = FALSE) +
  scale_color_viridis_d(option = "plasma", end = .7)

g3 <- ggplot(data,
       aes(x = standLRT, 
           y = normexam, 
           group = vr, 
           col = vr)) + 
  geom_point(alpha = 0.2) + 
  geom_smooth(method = lm, se = FALSE) +
  scale_color_viridis_d(option = "plasma", end = .7)

g4 <- ggplot(data,
       aes(x = standLRT, 
           y = normexam, 
           group = schgend, 
           col = schgend)) + 
  geom_point(alpha = 0.2) + 
  geom_smooth(method = lm, se = FALSE) +
  scale_color_viridis_d(option = "plasma", end = .7)

ggarrange(g1, g2, g3, g4, ncol = 2, nrow = 2)

7 Modelo 1: Baseline

# PRÉ-PROCESSAMENTO DOS DADOS -------------
df <- data %>% 
  select(-c(school, student, type)) %>% 
  mutate_if(is.character, as.factor) %>% 
  na.omit()

# DIVISÃO EM TREINO E TESTE ------------------------
set.seed(123)
data_split <- initial_split(df, prop = 3/4)

train_data <- training(data_split)
test_data <- testing(data_split)

# Especificação do Modelo ------------------------------------
especificacao_modelo <- linear_reg() %>%
  set_engine("lm") %>%
  set_mode("regression")

lm_fit <- especificacao_modelo %>% 
  fit(normexam ~., data = train_data)
tidy(lm_fit)

## # A tibble: 10 × 5
##    term          estimate std.error statistic  p.value
##    <chr>            <dbl>     <dbl>     <dbl>    <dbl>
##  1 (Intercept)     0.665     0.0930     7.14  1.17e-12
##  2 standLRT        0.383     0.0220    17.4   3.20e-64
##  3 intakemid 50%  -0.435     0.0414   -10.5   2.66e-25
##  4 intaketop 25%  -0.783     0.0706   -11.1   5.21e-28
##  5 sexM           -0.117     0.0414    -2.84  4.60e- 3
##  6 vrmid 50%      -0.205     0.0694    -2.96  3.14e- 3
##  7 vrtop 25%      -0.335     0.112     -2.98  2.92e- 3
##  8 schavg          0.620     0.115      5.38  7.97e- 8
##  9 schgendgirls    0.0258    0.0638     0.404 6.86e- 1
## 10 schgendmixed   -0.121     0.0502    -2.41  1.62e- 2

Sumário do modelo:

summary(lm_fit$fit)

## 
## Call:
## stats::lm(formula = normexam ~ ., data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.75368 -0.50948  0.00797  0.52507  2.71049 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    0.66455    0.09304   7.143 1.17e-12 ***
## standLRT       0.38263    0.02202  17.378  < 2e-16 ***
## intakemid 50% -0.43483    0.04142 -10.499  < 2e-16 ***
## intaketop 25% -0.78314    0.07058 -11.096  < 2e-16 ***
## sexM          -0.11747    0.04142  -2.836  0.00460 ** 
## vrmid 50%     -0.20532    0.06944  -2.957  0.00314 ** 
## vrtop 25%     -0.33503    0.11246  -2.979  0.00292 ** 
## schavg         0.62027    0.11524   5.382 7.97e-08 ***
## schgendgirls   0.02577    0.06377   0.404  0.68617    
## schgendmixed  -0.12077    0.05018  -2.407  0.01616 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7803 on 2714 degrees of freedom
## Multiple R-squared:  0.3991, Adjusted R-squared:  0.3971 
## F-statistic: 200.3 on 9 and 2714 DF,  p-value: < 2.2e-16

# Métricas de performance na base de treino
glance(lm_fit)

## # A tibble: 1 × 12
##   r.squared adj.r.squared sigma statistic   p.value    df logLik   AIC   BIC
##       <dbl>         <dbl> <dbl>     <dbl>     <dbl> <dbl>  <dbl> <dbl> <dbl>
## 1     0.399         0.397 0.780      200. 2.28e-292     9 -3184. 6391. 6456.
## # … with 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>

Visualizando os coeficientes das variáveis:

# Coeficientes das variáveis 
tidy(lm_fit) %>% 
  dwplot(dot_args = list(size = 2, color = "black"),
         whisker_args = list(color = "black"),
         vline = geom_vline(xintercept = 0, colour = "grey50", linetype = 2))

Gráficos de diagnóstico do modelo treinado:

gglm(lm_fit$fit)

No primeiro gráfico vemos que os resíduos estão espalhados de forma aleatória em torno da linha horizontal, sem seguir nenhum padrão específico.
O Q-Q plot confirma a suposição de que os resíduos são normalmente distribuídos.
O gráfico de Scale-Location também mostram os resíduos espalhados aleatoriamente, sem padrão específico.
No gráfico de Residual vs. Leverage, não observa-se pontos de alavancagem ou alta influência nos resíduos.
De maneira geral, o modeloé bom para ajustar o normexam em função das variáveis dos alunos e das escolas. Mas, temos que ver como ele se comporta na base de teste.

7.1 Importância das variáveis

vip(lm_fit)

O standLRT e o intake são as variáveis mais importantes para o modelo, de acordo com os scores atribuídos, muito mais que as variáveis ligadas ao sexo dos alunos. O tipo de escola, se de um único sexo, ou mista, não aparece na lista de importância do modelo.

7.2 Avaliando o desempenho do modelo na base de teste

Agora, vamos usar o modelo treinado para fazer predição na base de teste.

test_results <- predict(lm_fit, new_data = test_data) %>% 
  bind_cols(test_data)

7.3 Calculando as métricas: RMSE e \(R^2\) na base de teste

# RMSE
rmse(test_results,
     truth = normexam,
     estimate = .pred)

## # A tibble: 1 × 3
##   .metric .estimator .estimate
##   <chr>   <chr>          <dbl>
## 1 rmse    standard       0.747

# R^2
rsq(test_results,
     truth = normexam,
     estimate = .pred)

## # A tibble: 1 × 3
##   .metric .estimator .estimate
##   <chr>   <chr>          <dbl>
## 1 rsq     standard       0.423

7.4 Gráfico dos preditos vs observados

ggplot(data = test_results,
       mapping = aes(x = .pred, y = normexam)) +
  geom_point(color = '#006EA1', alpha = 0.5) +
  geom_abline(intercept = 0, slope = 1, color = 'orange') +
  labs(title = 'Linear Regression Results - Modelo 1',
       x = 'Predicted normexam',
       y = 'Actual normexam')

# Salvando o modelo final
modelo_final <- especificacao_modelo %>% 
  fit(normexam ~., data = df)

saveRDS(modelo_final, file = "final_model_normexam.rds")

# Uso em predições ---------------------------------- 
# <- readRDS("final_model_normexam.rds")
# new_data <- test_data[1, -1]
# predict(reg_lin_mdl, new_data = new_data)

8 Modelo 2: Interações entre variáveis

Neste modelo fazemos a inclusão de termos de interação entre variáveis. Especificamente:

standLRT com intake, vr e sex.

# Treino e Teste
set.seed(1234)
data_split <- initial_split(df, prop = 3/4)
train_data <- training(data_split)
test_data <- testing(data_split)

# Feature Engineering
recipe <- recipe(normexam ~., data = train_data) %>% 
  step_dummy(intake, sex, vr, schgend) %>% 
  step_interact(terms = ~ standLRT:starts_with("intake") +
                  standLRT:starts_with("vr") +
                  standLRT:starts_with("sex"))

# Especificando o modelo
lm_model <- linear_reg() %>%
  set_engine("lm") %>%
  set_mode("regression")

# Criando o workflow
exam_wf <- workflow() %>% 
  add_model(lm_model) %>% 
  add_recipe(recipe)

# Executando o workflow
exam_fit <- exam_wf %>% 
  last_fit(split = data_split)

# Métricas na base de teste
exam_fit %>% 
  collect_metrics()

## # A tibble: 2 × 4
##   .metric .estimator .estimate .config             
##   <chr>   <chr>          <dbl> <chr>               
## 1 rmse    standard       0.787 Preprocessor1_Model1
## 2 rsq     standard       0.392 Preprocessor1_Model1

# Predições
predicoes_res <- exam_fit %>% 
  collect_predictions()
predicoes_res

## # A tibble: 908 × 5
##    id                 .pred  .row normexam .config             
##    <chr>              <dbl> <int>    <dbl> <chr>               
##  1 train/test split  0.703      1   1.51   Preprocessor1_Model1
##  2 train/test split -0.268      8  -0.0621 Preprocessor1_Model1
##  3 train/test split  0.565     14   0.822  Preprocessor1_Model1
##  4 train/test split  0.439     18   1.81   Preprocessor1_Model1
##  5 train/test split -0.327     20  -1.12   Preprocessor1_Model1
##  6 train/test split -1.07      33  -1.22   Preprocessor1_Model1
##  7 train/test split  0.639     43   0.261  Preprocessor1_Model1
##  8 train/test split  0.0554    44   0.134  Preprocessor1_Model1
##  9 train/test split  1.04      49   1.90   Preprocessor1_Model1
## 10 train/test split  0.659     52   0.968  Preprocessor1_Model1
## # … with 898 more rows

# Importância das variáveis
train_data_baked <- recipe %>% 
  prep() %>% 
  bake(new_data = train_data)
vip(lm_model %>% 
      fit(normexam ~., data = train_data_baked)
    )

Gráfico do preditos vs observados:

ggplot(data = predicoes_res,
       mapping = aes(x = .pred, y = normexam)) +
  geom_point(color = '#006EA1', alpha = 0.5) +
  geom_abline(intercept = 0, slope = 1, color = 'orange') +
 labs(title = 'Linear Regression Results - Modelo 2',
       x = 'Predicted normexam',
       y = 'Actual normexam')

O Modelo 2 não performa melhor que o Modelo 1, como pode ser visto na tabela abaixo.

Métricas	Modelo 1	Modelo 2
RMSE	0.747	0.787
\(R^2\)	0.423	0.392

O Modelo 2 tem um erro maior (RMSE) e uma menor explicabilidade (\(R^2\) ) das notas normalizadas em termos das variáveis explicativas.

9 Próximos passos

Fazer cross-validation
Explorar mais processos de feature engeneering

Usar o scikitlearn
Usar árvore de decisão

10 Conclusão

De acordo com o modelo de regressão linear gerado, as variáveis de maior relevância para prever o resultado no exame normalizado são standLRT, intake, schavg com maiores scores. Seguidas de vr e sex, com scores um pouco menores. O gênero da escola (schgend) foi a variável que teve menor importância na predição na nota no exame normalizado. Então, podemos dizer que o tipo de escola afeta muito pouco na predição do resultado do exame.

Case Técnico

Processo Seletivo Murabei

Edneide Ramalho

06/04/2022