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

  1. Carregar as bases de dados individuais.
  2. 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.
  1. Fazer uma análise exploratória dos dados para verificar o comportamento dos dados e buscar insights dos mesmos.
  2. Decidir qual a melhor metodologia para responder à pergunta relacionada ao nosso objetivo
  3. Ajustar um modelo apropriado para prever as notas normalizadas dos alunos e responder às perguntas de interesse.
  4. Salvamento do modelo final para ser usado em predições de novos dados.

3 Pacotes utilizados

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 para explorar os resultados.

  • Carregamento dos dados já limpos:
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"…
  • Olhada geral nos dados:s
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 empty n_unique whitespace
intake 147 0.96 7 10 0 3 0
sex 173 0.96 1 1 0 2 0
vr 134 0.97 7 10 0 3 0
schgend 0 1.00 4 5 0 3 0
type 0 1.00 3 4 0 2 0

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

Deste gráfico, podemos tirar alguns insights:

  • 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 escola (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 gêneros.

  • 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 Tabela descritiva para os dados da nossa base:

library(gtsummary)
data %>% 
  dplyr::select(-c(school, student)) %>% 
  tbl_summary()
Characteristic N = 3,6321
normexam 0.00 (-0.70, 0.68)
standLRT 0.04 (-0.62, 0.62)
intake
bottom 25% 1,033 (28%)
mid 50% 2,124 (58%)
top 25% 475 (13%)
sex
F 2,156 (59%)
M 1,476 (41%)
vr
bottom 25% 581 (16%)
mid 50% 2,047 (56%)
top 25% 1,004 (28%)
schavg -0.02 (-0.15, 0.21)
schgend
boys 467 (13%)
girls 1,250 (34%)
mixed 1,915 (53%)
type
Mxd 1,915 (53%)
Sngl 1,717 (47%)

1 Median (IQR); n (%)

7 Perguntas

Antes de ajustar um modelo aos dados, vamos levantar algumas perguntas adicionais que nos ajudar a compreender os dados e gerar insights para uma melhor estratégia de modelagem.

7.1 O resultado do exame normalizado depende do gênero da escola?

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

Ao analisar o gráfico, verificamos que a média nas notas dos exames é maior para escolas só de garotas (nota média normalizada = 0.144), seguida pela escola só de garotos (nota média normalizada = 0.029), e finalmente pela escola de gênero misto (nota média normalizada = -0.099). Isso nos dá indícios de que estudantes de escolas de gênero único parece perfomar melhor em seus exames. Esta hipótese só será validada ao ajustarmos o modelo final, onde a adição das outras variáveis da nossa base será feita.

7.2 O resultado do exame normalizado é diferente para garotos e garotas, de maneira geral?

data %>% 
  group_by(sex) %>% 
  get_summary_stats(normexam, type = "full")
## # A tibble: 2 × 14
##   sex   variable     n   min   max median     q1    q3   iqr   mad   mean    sd
##   <chr> <chr>    <dbl> <dbl> <dbl>  <dbl>  <dbl> <dbl> <dbl> <dbl>  <dbl> <dbl>
## 1 F     normexam  2156 -3.06  3.67  0.074 -0.555 0.679  1.23 0.932  0.089 0.975
## 2 M     normexam  1476 -3.67  2.53 -0.129 -0.853 0.611  1.46 1.07  -0.126 1.02 
## # … with 2 more variables: se <dbl>, ci <dbl>
# Gráfico
data %>% 
  ggplot(aes(x = sex, y = normexam)) +
  geom_boxplot() +
  stat_summary(fun.y = mean, shape = 15, color = "blue")

Considerando somento o gênero dos estudantes, as garotas tiram, em média, notas melhores no exame do que os garotos (nota normalizada média = 0.089 para as garotas vs. 0.-126 para os garotos). Isto nos leva às próxima perguntas, se realmente o gênero é importante para a performance no teste, ou o tipo da escola é que se sobressai.

7.3 Os garotos de escolas de gênero único tem notas diferentes dos garotos de escola mista?

data %>% 
  # Filtrando só para garotos
  dplyr::filter(sex == "M") %>% 
  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  1009 -3.67  2.41 -0.198 -0.853 0.478  1.33  1.00 -0.199 0.991
## 2 Sngl  normexam   467 -2.75  2.53  0.074 -0.738 0.822  1.56  1.15  0.029 1.07 
## # … with 2 more variables: se <dbl>, ci <dbl>
# Gráfico
data %>% 
  # Filtrando só para garotos
  dplyr::filter(sex == "M") %>% 
  ggplot(aes(x = type, y = normexam)) +
  geom_boxplot() +
  stat_summary(fun.y = mean, shape = 15, color = "blue") +
  ggtitle("Boys")

Os resultados nos mostram que, em média, os garotos de escolas de gênero único tiram melhores notas normalizadas que os garotos de escolas mistas. A nota média normalizada para escola mista é -0.199, enquanto que para escolas só de garotos é de 0.029.

7.4 As garotas de escolas de gênero único tem notas diferentes das garotas de escola mista?

data %>% 
  # Filtrando só para garotas
  dplyr::filter(sex == "F") %>% 
  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   906 -2.75  2.92  0.004 -0.623 0.679  1.30 0.93  0.012 0.955
## 2 Sngl  normexam  1250 -3.06  3.67  0.134 -0.493 0.747  1.24 0.929 0.144 0.987
## # … with 2 more variables: se <dbl>, ci <dbl>
# Gráfico
data %>% 
  # Filtrando só para garotas
  dplyr::filter(sex == "F") %>% 
  ggplot(aes(x = type, y = normexam)) +
  geom_boxplot() +
  stat_summary(fun.y = mean, shape = 15, color = "blue") +
  ggtitle("Girls")

Os resultados nos mostram mais uma vez que, em média, estudantes de escola de gênero único parecem tirar notas melhores. As garotas de escolas de gênero único tiram melhores notas normalizadas que as garotas de escolas mistas, notas média de 0.144 vs. 0.012, repectivamente.

Muitas outras perguntas podem ser levantadas, incluindo as outras variáveis na base de dados, porém, focamos na nossa variável dependente normexam e no gênero da escola schgend para saber seu impacto nas notas normalizadas dos alunos. A seguir, vamos ajustar um modelo incluindo as outras variáveis para verificar seus respectivos pesos em predizer a nota normalizada.

8 Modelo

Como queremos prever uma variál quantitativa, podemos pensar em usar uma análise de regressão linear tradicional. Porém, nossos dados não atendem à premissa de independência das observações. Isso ocorre por que as informações da nossa base de dados estão aninhadas, isto é, estudantes podem ser do gênero feminino ou masculino, os mesmos podem estudar numa escola do tipo gênero misto e gênero único; o tipo da escola, por sua vez é subdividido em escolas só de garotos, escolas só de garotas e escolas de gênero misto. Assim, por exemplo, o conjunto de notas dos estudantes pode variar por fatores como a escola em que estuda, de que tipo é a escola, entre outros. Então podemos perguntar, os estudantes são realmente independentes uns dos outros na mesma escola ou num mesmo tipo de escola (só garotos, só garotas, gênero misto)? A resposta é, provavelmente não!

Os gráficos abaixo nos dão uma intuição sobre as regressões corrigidas pelo gênero das escolas. Aqui, utilizamos apenas variáveis numéricas para mostrar como os interceptos e as inclinações das retas se modificam para as classes da variável schgend.

g1 <- ggplot(data,
       aes(x = schavg, 
           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)

g2 <- 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, ncol = 1, nrow = 2)

Sendo, assim, vamos utilizar uma regressão linear de efeitos mistos, para modelar a nota normalizada normexam usando todas as variáveis preditoras como efeito fixo, e a variável schgend como efeito aleatório (por que queremos levar em conta a sua variabilidade no resultado do modelo).

library(lme4)
library(broom.mixed)

lmer_normexam = lmer(normexam ~ standLRT + intake + sex + vr + schavg + 0 + (1 | schgend), data = data)

lmer_normexam
## Linear mixed model fit by REML ['lmerMod']
## Formula: normexam ~ standLRT + intake + sex + vr + schavg + 0 + (1 | schgend)
##    Data: data
## REML criterion at convergence: 8463.781
## Random effects:
##  Groups   Name        Std.Dev.
##  schgend  (Intercept) 0.09224 
##  Residual             0.77182 
## Number of obs: 3632, groups:  schgend, 3
## Fixed Effects:
##         standLRT  intakebottom 25%     intakemid 50%     intaketop 25%  
##           0.3738            0.6129            0.1683           -0.1947  
##             sexM         vrmid 50%         vrtop 25%            schavg  
##          -0.1299           -0.1728           -0.2374            0.5357
summary(lmer_normexam)
## Linear mixed model fit by REML ['lmerMod']
## Formula: normexam ~ standLRT + intake + sex + vr + schavg + 0 + (1 | schgend)
##    Data: data
## 
## REML criterion at convergence: 8463.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.5660 -0.6468  0.0167  0.6636  3.4745 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  schgend  (Intercept) 0.008508 0.09224 
##  Residual             0.595711 0.77182 
## Number of obs: 3632, groups:  schgend, 3
## 
## Fixed effects:
##                  Estimate Std. Error t value
## standLRT          0.37382    0.01889  19.786
## intakebottom 25%  0.61292    0.08631   7.101
## intakemid 50%     0.16829    0.08386   2.007
## intaketop 25%    -0.19470    0.09314  -2.090
## sexM             -0.12990    0.03431  -3.786
## vrmid 50%        -0.17283    0.05940  -2.910
## vrtop 25%        -0.23738    0.09674  -2.454
## schavg            0.53574    0.09983   5.367
## 
## Correlation of Fixed Effects:
##             stnLRT intkb25% int50% intkt25% sexM   vrm50% vrt25%
## intkbttm25% -0.185                                              
## intakemd50%  0.040  0.914                                       
## intaketp25%  0.268  0.778    0.875                              
## sexM         0.026 -0.236   -0.245 -0.228                       
## vrmid 50%   -0.020 -0.662   -0.698 -0.629    0.060              
## vrtop 25%   -0.023 -0.655   -0.687 -0.628    0.061  0.898       
## schavg      -0.065  0.577    0.609  0.564   -0.056 -0.784 -0.901

8.1 Extraindo os coeficientes

lmer_coef <- tidy(lmer_normexam, conf.int = TRUE)
print(lmer_coef)
## # A tibble: 10 × 8
##    effect   group    term        estimate std.error statistic conf.low conf.high
##    <chr>    <chr>    <chr>          <dbl>     <dbl>     <dbl>    <dbl>     <dbl>
##  1 fixed    <NA>     standLRT      0.374     0.0189     19.8   0.337      0.411 
##  2 fixed    <NA>     intakebott…   0.613     0.0863      7.10  0.444      0.782 
##  3 fixed    <NA>     intakemid …   0.168     0.0839      2.01  0.00393    0.333 
##  4 fixed    <NA>     intaketop …  -0.195     0.0931     -2.09 -0.377     -0.0122
##  5 fixed    <NA>     sexM         -0.130     0.0343     -3.79 -0.197     -0.0626
##  6 fixed    <NA>     vrmid 50%    -0.173     0.0594     -2.91 -0.289     -0.0564
##  7 fixed    <NA>     vrtop 25%    -0.237     0.0967     -2.45 -0.427     -0.0478
##  8 fixed    <NA>     schavg        0.536     0.0998      5.37  0.340      0.731 
##  9 ran_pars schgend  sd__(Inter…   0.0922   NA          NA    NA         NA     
## 10 ran_pars Residual sd__Observ…   0.772    NA          NA    NA         NA
# Efeito fixo
fixef(lmer_normexam)
##         standLRT intakebottom 25%    intakemid 50%    intaketop 25% 
##        0.3738171        0.6129210        0.1682862       -0.1946998 
##             sexM        vrmid 50%        vrtop 25%           schavg 
##       -0.1299049       -0.1728312       -0.2373789        0.5357355
# Efeito aleatório
ranef(lmer_normexam)
## $schgend
##       (Intercept)
## boys   0.05310138
## girls  0.04773605
## mixed -0.10083743
## 
## with conditional variances for "schgend"
  • Escolas só de garotos ou só de garotas apresentam um intercepto similar, cerca de 0.05, enquanto escolas mistas apresentam um intercepto menor, -0.1. Isso quer dizer que, se todas as outras variáveis são mantidas constantes, esperamos que a nota do estudante seja maior em escolas de gênero único.
# Plot results
lmer_coef %>%
    filter(effect == "fixed" & term != "(Intercept)") %>%
    ggplot(., aes(x = term, y = estimate,
                  ymin = conf.low, ymax = conf.high)) +
    geom_hline(yintercept = 0, color = 'red') + 
    geom_point() +
    geom_linerange() +
    coord_flip() +
    theme_bw() +
    ylab("Coefficient estimate and 95% CI") +
    xlab("Regression coefficient")

Todas as variáveis listadas no gráfico acima são importantes para a predição do normexam. As que estão do lado esquerdo da reta vermelha, tem um efeito negativo no normexam, ou seja, quando elas aumentam o normexam diminui; e as que estão do lado direito têm um efeito positivo, ou seja, quando elas aumentam, o normexam também aumenta. Para as variáveis categóricas, funciona de maneira similar, por exemplo, se todas as outras variáveis forem mantidas constantes, um estudante com vrmid50% terá uma nota melhor que um estudante no vrtop25%.

8.2 Predição

Podemos usar o modelo gerado para fazer predições:

pred_df <- tibble(
   actual = data$normexam,
  .pred = predict(lmer_normexam))

Plotando os valores observados vs. os valores preditos pelo modelo, temos:

ggplot(data = pred_df,
       mapping = aes(x = .pred, y = actual)) +
  geom_point(color = '#006EA1', alpha = 0.5) +
  geom_abline(intercept = 0, slope = 1, color = 'orange') +
  labs(title = 'Model Result',
       x = 'Predicted normexam',
       y = 'Actual normexam')

9 Próximos passos

  • Explorar outros efeitos aleatórios e fixos.

  • Gerar outros modelos para comparação de eficiência.

10 Conclusão

De acordo com o modelo de regressão linear de efeitos mistos gerado, escolas de gênero único são mais eficientes na formação de seus alunos de acordo com o exame normalizado. Outras variáveis também tem forte importância na determinação do sucesso do estudante no exame, com destaque para schavg , standLRT e intake.