EST583 – Modelagem de Relações Usando Regressão Linear

UFMG – Especialização em Estatística Computacional Aplicada - Prof. Guilherme Lopes de Oliveira

Author

Igor Souza

Published

03/10/2023

Introdução

O relatório apresentado a seguir foi produzido para fins acadêmicos. Trata-se do trabalho prático final da disciplina de Modelagem de Relações Usando Regressão Linear. O objetivo é aplicar os conceitos referentes a análise de regressão linear simples e múltipla, estimativa e interpretação dos coeficientes de regressão, predição, análise de resíduos e multicolinearidade e seleção de modelos; todos conceitos vistos em sala de aula.

Dados

Para o trabalho final foi escolhido os microdados do ENEM 2022 divulgados pelo Instituto Nacional de Estudos e Pesquisas Educacionais Anísio Teixeira (Inep).

A finalidade da modelagem é predizer o desempenho de um aluno na prova objetiva do ENEM dado suas características individuais, escolares e também dado o contexto socioeconômico em que sua família se encontra. Toda a análise foi realizada na linguagem R (v. 4.3.1).1

Antes de qualquer análise, um trabalho de seleção, filtro e adequação de variáveis foi realizado. Para o estudo, foram considerados somente alunos que realizaram a prova na cidade de Belo Horizonte. Outros critérios de seleção incluiram apenas os alunos que:

  1. Eram considerados jovens (20 anos de idade ou menos) na data da prova;

  2. Informaram o estado civil;

  3. Declaram ter cor branca, preta ou parda;

  4. Possuem nacionalidade brasileira;

  5. Declararam o tipo de escola onde cursaram o Ensino Médio (pública ou privada);

  6. Cursavam o Ensino Médio em instituições de ensino regular na data da prova;

  7. Estudavam em escolas em atividade na data da prova;

  8. Estavam presentes em todas as 4 provas objetivas;

  9. Não zeraram nenhuma das 4 provas objetivas;

  10. Entregaram uma redação considerada “sem problemas”;

  11. Declararam o nível de instrução da mãe.

library(tidyverse)
library(car)
library(lmtest)
library(performance)
library(gt)
library(gtsummary)
library(skimr)

As primeiras 10 observações do data.frame.

alu_solteiro alu_homem alu_branco escol_publica soc_mae_grad soc_renda_baixa soc_renda_media soc_renda_alta soc_internet prov_nota_natur prov_nota_human prov_nota_ling prov_nota_mat prov_nota_redacao prov_nota_objetiva prov_nota_total
0 1 1 0 1 1 1 1 0 509.1 518.3 533.8 561.4 520 2122.6 2642.6
0 1 1 0 1 1 1 1 0 530.3 550.4 585.9 594.7 900 2261.3 3161.3
0 0 0 0 1 1 1 1 0 548.8 560.2 543.7 632.7 920 2285.4 3205.4
0 1 0 0 1 1 1 1 0 533.4 527.8 538.0 578.3 540 2177.5 2717.5
0 1 0 0 1 0 1 1 0 485.9 408.2 579.7 592.4 840 2066.2 2906.2
0 0 1 0 1 1 0 1 0 468.8 556.7 587.0 619.3 620 2231.8 2851.8
0 0 1 0 1 1 1 1 0 593.7 642.9 614.5 689.1 600 2540.2 3140.2
0 1 1 0 0 1 0 1 0 586.4 596.3 575.7 714.3 860 2472.7 3332.7
0 0 1 0 1 1 1 1 0 405.3 476.3 373.6 366.5 380 1621.7 2001.7
0 1 1 0 1 1 1 1 0 459.3 475.8 427.5 411.5 640 1774.1 2414.1

Visão geral da base de dados

Rows: 7,377
Columns: 16
$ alu_solteiro       <fct> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ alu_homem          <fct> 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, …
$ alu_branco         <fct> 1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, …
$ escol_publica      <fct> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, …
$ soc_mae_grad       <fct> 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, …
$ soc_renda_baixa    <fct> 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, …
$ soc_renda_media    <fct> 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 1, 0, 1, 1, 0, 1, …
$ soc_renda_alta     <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, …
$ soc_internet       <fct> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ prov_nota_natur    <dbl> 509.1, 530.3, 548.8, 533.4, 485.9, 468.8, 593.7, 58…
$ prov_nota_human    <dbl> 518.3, 550.4, 560.2, 527.8, 408.2, 556.7, 642.9, 59…
$ prov_nota_ling     <dbl> 533.8, 585.9, 543.7, 538.0, 579.7, 587.0, 614.5, 57…
$ prov_nota_mat      <dbl> 561.4, 594.7, 632.7, 578.3, 592.4, 619.3, 689.1, 71…
$ prov_nota_redacao  <int> 520, 900, 920, 540, 840, 620, 600, 860, 380, 640, 5…
$ prov_nota_objetiva <dbl> 2122.6, 2261.3, 2285.4, 2177.5, 2066.2, 2231.8, 254…
$ prov_nota_total    <dbl> 2642.6, 3161.3, 3205.4, 2717.5, 2906.2, 2851.8, 314…

Dicionário de variáveis

alu_solteiro: variável dummy que indica se o aluno é solteiro (0) ou não (1).

alu_homem: variável dummy que indica se o aluno é homem (0) ou mulher (1).

alu_branco: variável dummy que indica se o aluno é branco (0) ou negro (1).

escol_publica: variável dummy que indica se a escola em que o aluno estuda é pública (0) ou privada (1).

soc_mae_grad: variável dummy que indica se a mãe do aluno concluiu, no mínimo, a faculdade (0) ou não (1).

soc_renda_baixa*: variável dummy que indica se a renda familiar do aluno é baixa (≥ R$ 2.424,01 e ≤ R$ 3.636,00) [0] ou não [1].

soc_renda_media: variável dummy que indica se a renda familiar do aluno é média (≥ R$ 3.636,01 e ≤ R$ 9.696,00) [0] ou não [1].

soc_renda_alta: variável dummy que indica se a renda familiar do aluno é alta (≥ R$ R$ 9.696,01) [0] ou não [1].

soc_internet: variável dummy que indica se o aluno tem acesso à internet (0) ou não (1).

prov_nota_natur: variável numérica referente à nota da prova de ciências da natureza (não utilizada nos modelos).

prov_nota_human: variável numérica referente à nota da prova de ciências humanas (não utilizada nos modelos).

prov_nota_ling: variável numérica referente à nota da prova de linguagens e códigos (não utilizada nos modelos).

prov_nota_mat: variável numérica referente à nota da prova de matemática (não utilizada nos modelos).

prov_nota_redacao: variável numérica referente à nota da redação (não utilizada nos modelos).

prov_nota_objetiva: variável numérica referente ao total da nota das provas objetivas.

prov_nota_total: variável numérica referente ao total da nota das provas - incluindo a redação (não utilizada nos modelos).

*A renda familiar foi divida em \(k\) = 4 categorias (muito baixa, baixa, média e alta). No modelo, foram incluídas \(k - 1\) covariáveis: soc_renda_baixa, soc_renda_media e soc_renda_alta. Portanto, renda muito baixa (≤ R$ 2.424,00) é o grupo de referência.

Explorar

Análise descritiva

Sumário das variáveis numéricas (mesmo aquelas que não fazem parte diretamente dos modelos).

Data summary
Name dados_num
Number of rows 7377
Number of columns 7
Key NULL
_______________________
Column type frequency:
numeric 7
________________________
Group variables None

Variable type: numeric

skim_variable média desvio padrão máximo mínimo
prov_nota_natur 546.95 78.15 853.4 361.1
prov_nota_human 583.79 72.38 827.2 315.8
prov_nota_ling 567.79 66.20 778.2 302.8
prov_nota_mat 635.02 124.53 983.9 336.8
prov_nota_redacao 748.27 150.45 980.0 220.0
prov_nota_objetiva 2333.55 299.87 3289.3 1475.2
prov_nota_total 3081.82 416.30 4249.3 1817.3

Correlação de Pearson.

Variável prov_nota_natur prov_nota_human prov_nota_ling prov_nota_mat prov_nota_redacao prov_nota_objetiva prov_nota_total
prov_nota_natur 1.000 0.702 0.638 0.693 0.551 0.858 0.817
prov_nota_human 0.702 1.000 0.769 0.708 0.617 0.888 0.862
prov_nota_ling 0.638 0.769 1.000 0.639 0.590 0.838 0.817
prov_nota_mat 0.693 0.708 0.639 1.000 0.604 0.908 0.872
prov_nota_redacao 0.551 0.617 0.590 0.604 1.000 0.673 0.846
prov_nota_objetiva 0.858 0.888 0.838 0.908 0.673 1.000 0.964
prov_nota_total 0.817 0.862 0.817 0.872 0.846 0.964 1.000

Correlação ponto-bisserial (prov_nota_objetiva x variáveis qualitativas/dicotômicas).

Variável Ponto-bisserial
soc_mae_grad 0.373
alu_branco 0.292
soc_renda_baixa 0.120
alu_homem 0.091
soc_internet 0.081
alu_solteiro −0.013
soc_renda_media −0.022
soc_renda_alta −0.411
escol_publica −0.444

Visualizar

Histograma

Boxplot

Regressão linear múltipla

Seleção de covariáveis

Dois métodos de seleção de covariáveis foram empregados: manual e stepwise. A estratégia é aplicar mais de um método seleção, fazer comparações dos resultados dos métodos e então escolher aquele modelo que melhor se adequa ao objetivo proposto pelo trabalho.

Manual

Nesta abordagem foi criado um modelo com todas as variáveis disponíveis na base de dados. Em seguida, foram removidas aquelas variáveis que não são estatisticamente significativas (valor-p > 0,05).

Code 
# Regressão com todas as variáveis
rlm_todas_var <- lm(prov_nota_objetiva ~ ., data = dados_num_cat)

# -alu_solteiro
rlm_todas_var_upd_1 <- update(rlm_todas_var, ~ . -alu_solteiro)

Seleção manual de variáveis

k = 9
Characteristic Beta 95% CI1 p-value VIF1
alu_solteiro 1.0
    0 — —
    1 -17 -128, 93 0.8
alu_homem 1.0
    0 — —
    1 -40 -52, -29 <0.001
alu_branco 1.2
    0 — —
    1 -67 -79, -54 <0.001
escol_publica 1.6
    0 — —
    1 125 110, 139 <0.001
soc_mae_grad 1.4
    0 — —
    1 -76 -90, -63 <0.001
soc_renda_baixa 1.4
    1 — —
    0 62 42, 81 <0.001
soc_renda_media 1.8
    1 — —
    0 102 85, 119 <0.001
soc_renda_alta 2.5
    1 — —
    0 204 185, 223 <0.001
soc_internet 1.0
    0 — —
    1 -62 -120, -3.8 0.037
1 CI = Confidence Interval, VIF = Variance Inflation Factor

k = 8
Characteristic Beta 95% CI1 p-value VIF1
alu_homem 1.0
    0 — —
    1 -40 -52, -29 <0.001
alu_branco 1.2
    0 — —
    1 -67 -79, -54 <0.001
escol_publica 1.6
    0 — —
    1 125 110, 139 <0.001
soc_mae_grad 1.4
    0 — —
    1 -76 -90, -63 <0.001
soc_renda_baixa 1.4
    1 — —
    0 62 42, 81 <0.001
soc_renda_media 1.8
    1 — —
    0 102 85, 119 <0.001
soc_renda_alta 2.5
    1 — —
    0 204 185, 223 <0.001
soc_internet 1.0
    0 — —
    1 -62 -120, -3.8 0.037
1 CI = Confidence Interval, VIF = Variance Inflation Factor

Seguindo o critério escolhido, a variável alu_solteiro foi removida do modelo. soc_renda_media também não pode ser considerada estatisticamente significativa. Entretanto, devido ao efeito aditivo das variáveis de renda, a variável foi mantida. Na tabela acima pode-se observar que todas as demais variáveis são estatisticamente significativas - ao nível de 95% de confiança - e não apresentam multicolinearidade.

Não foi possível realizar o teste de Shapiro devido ao tamanho da amostra (n > 5.000).

Sobre os resultados dos testes de Breusch-Pagan e de Durbin-Watson, tem-se:

\(H_{0}:\) homocedasticidade está presente (os resíduos possuem variância constante).

\(H_{A}:\) heterocedasticidade está presente (os resíduos não possuem variância constante).


    studentized Breusch-Pagan test

data:  rlm_todas_var_upd_1
BP = 63.206, df = 8, p-value = 1.091e-10

\(H_{0}:\) os resíduos não estão correlacionados.

\(H_{A}:\) os resíduos estão correlacionados.


    Durbin-Watson test

data:  rlm_todas_var_upd_1
DW = 1.8816, p-value = 3.644e-07
alternative hypothesis: true autocorrelation is not 0

Stepwise

O algoritmo de seleção automática stepwise foi escolhido. Nenhuma alteração foi feita no modelo resultante.

Code 
# Stepwise (BIC)
rlm_spw <- step(
  rlm_todas_var, 
  direction = "both", 
  k         = log(nrow(dados_num_cat)),
  trace     = F
)
Characteristic Beta 95% CI1 p-value VIF1
alu_homem 1.0
    0 — —
    1 -41 -52, -29 <0.001
alu_branco 1.2
    0 — —
    1 -67 -80, -54 <0.001
escol_publica 1.6
    0 — —
    1 125 111, 139 <0.001
soc_mae_grad 1.4
    0 — —
    1 -77 -90, -63 <0.001
soc_renda_baixa 1.3
    1 — —
    0 63 44, 82 <0.001
soc_renda_media 1.8
    1 — —
    0 103 86, 120 <0.001
soc_renda_alta 2.4
    1 — —
    0 205 186, 224 <0.001
1 CI = Confidence Interval, VIF = Variance Inflation Factor

\(H_{0}:\) homocedasticidade está presente (os resíduos possuem variância constante).

\(H_{A}:\) heterocedasticidade está presente (os resíduos não possuem variância constante).


    studentized Breusch-Pagan test

data:  rlm_spw
BP = 60.484, df = 7, p-value = 1.209e-10

\(H_{0}:\) os resíduos não estão correlacionados.

\(H_{A}:\) os resíduos estão correlacionados.


    Durbin-Watson test

data:  rlm_spw
DW = 1.8809, p-value = 3.126e-07
alternative hypothesis: true autocorrelation is not 0

Heterocedasticidade

Como pode ser visto acima, rejeitamos \(H_{0}\) ao nível de 95% de confiança em ambos os modelos no teste de Breusch-Pagan (valor-p = 1.208509e-10). O gráfico abaixo mostra o problema.

Box-Cox

Para tentar solucionar o problema, a transformação de Box-Cox foi utilizada.

Code 
# Box-Cox
bxcx <- boxCox(rlm_spw, plotit = F)

# Obtendo o valor de lambda
lambda <- bxcx$x[which.max(bxcx$y)]

O resultado de \(\lambda\) é 0,9. Como o valor está próximo de 1 e dado que recomendações relacionadas ao método dizem que quando \(\lambda = 1\) então \(Y^{1} = Y\), nenhuma transformação na variável resposta foi realizada. A transformação de Box-Cox não foi suficiente para remediar o problema. A próxima tentativa é fazer uso do método de estimação por mínimos quadrados ponderados.

Mínimos Quadrados Ponderados

Devido a identificação de heterocedasticidade no modelo durante a fase de análise de resíduos e a não resolução do problema pela transformação de Box-Cox, o método de estimação por mínimos quadrados ponderados foi aplicado. Os pesos foram criados de forma que quanto maior a variância entre as notas das provas que dão origem à variável dependente prov_nota_objetiva, menor o peso. Ou seja, 1 / variância das notas das provas.

Code 
# Pesos = 1 / var(notas provas)
pesos <- dados_num |> 
  select(
    prov_nota_natur,
    prov_nota_human,
    prov_nota_ling,
    prov_nota_mat,
    prov_nota_redacao
  ) |> 
  rowwise() |> 
  mutate(
    var  = var(c_across(everything())),
    peso = 1 / var
  ) |>
  ungroup() |> 
  pull(peso)

# Regressão linear múltipla com as variáveis resultantes da abordagem stepwise
rlm_spw_mqp <- lm(
  prov_nota_objetiva ~ 
  alu_homem + 
  alu_branco + 
  escol_publica + 
  soc_mae_grad + 
  soc_renda_baixa + 
  soc_renda_media +
  soc_renda_alta,
  weights = pesos, 
  data    = dados_num_cat
)

O resultado é:

Modelo de regressão e tabela ANOVA

Modelo
Characteristic Beta 95% CI1 p-value VIF1
alu_homem 1.1
    0 — —
    1 -34 -44, -23 <0.001
alu_branco 1.1
    0 — —
    1 -48 -59, -37 <0.001
escol_publica 1.5
    0 — —
    1 120 105, 134 <0.001
soc_mae_grad 1.3
    0 — —
    1 -80 -92, -68 <0.001
soc_renda_baixa 1.4
    1 — —
    0 94 79, 109 <0.001
soc_renda_media 1.5
    1 — —
    0 80 65, 95 <0.001
soc_renda_alta 1.9
    1 — —
    0 178 159, 197 <0.001
1 CI = Confidence Interval, VIF = Variance Inflation Factor
ANOVA
Df Sum Sq Mean Sq F value Pr(>F)
alu_homem 1 1069 1069 87 0
alu_branco 1 6907 6907 562 0
escol_publica 1 13381 13381 1089 0
soc_mae_grad 1 5483 5483 446 0
soc_renda_baixa 1 417 417 34 0
soc_renda_media 1 24 24 2 0
soc_renda_alta 1 4249 4249 346 0
Residuals 7369 90553 12 NA NA

\(H_{0}:\) homocedasticidade está presente (os resíduos possuem variância constante).

\(H_{A}:\) heterocedasticidade está presente (os resíduos não possuem variância constante).


    studentized Breusch-Pagan test

data:  rlm_spw_mqp
BP = 0.4022, df = 7, p-value = 0.9997

\(H_{0}:\) os resíduos não estão correlacionados.

\(H_{A}:\) os resíduos estão correlacionados.

Não é possível fazer o dwtest() em regressões com M.Q.P.

Após a introdução dos pesos, o modelo apresenta homocedasticidade. Como consequência, os resultados dos testes de estatística são confiáveis.

Comparação

Abaixo, uma comparação do desempenho dos modelos utilizados. Foi escolhido o último modelo, resultado da seleção stepwise juntamente com a estimação por mínimos quadrados ponderados. A escolha se deve a significância estatística, a não multicolinearidade e a homocedasticidade da variância dos resíduos.

Code 
modelo_desemp <- compare_performance(
  rlm_todas_var_upd_1, 
  rlm_spw,
  rlm_spw_mqp
) 

modelo_desemp |> 
  gt()
Name Model AIC AIC_wt AICc AICc_wt BIC BIC_wt R2 R2_adjusted RMSE Sigma
rlm_todas_var_upd_1 lm 102426.8 0.7655165 102426.9 0.7650284 102495.9 0.09364705 0.3040834 0.3033278 250.1385 250.291235
rlm_spw lm 102429.2 0.2344835 102429.2 0.2349716 102491.3 0.90635295 0.3036714 0.3030099 250.2125 250.348330
rlm_spw_mqp lm 105907.7 0.0000000 105907.7 0.0000000 105969.9 0.00000000 0.2582629 0.2575583 255.5193 3.505482

Interpretação

No modelo final, a seguinte regressão linear foi ajustada aos dados:

\[ \hat{Y} \approx 2.208 - 34F_{1} - 48F_{2} + 120F_{3} - 80F_{4} + 94F_{5} + 80F_{6} + 178F_{7} + \epsilon \]

Isso significa que, partindo de uma pontuação inicial \(\beta_{0}\) de 2.208, as alunas sofrem uma redução de 34 pontos na prova. Fixados os parâmetros \(F_{1}\), \(F_{3}\), \(F_{4}\), \(F_{5}\), \(F_{6}\) e \(F_{7}\), os alunos de cor negra apresentam uma baixa de 48 pontos. Seguindo essa linha de raciocínio (mantendo os parâmetros fixos), os alunos provenientes de escolas privadas têm um desempenho superior em 120 pontos. Os alunos cujas mães não concluíram o ensino superior registram um decréscimo de 80 pontos. Em relação aos grupos de renda, em comparação com o grupo de renda considerada muito baixa, os alunos de renda baixa apresentam um aumento de 94 pontos, os de renda média um aumento de 80 pontos e os de renda alta um aumento de 178 pontos. Cabe ressaltar que todos esses valores são estimativas médias.

O modelo é estatisticamente significativo mas explica apenas 25,76% da variação nos dados.

Predição

Dois exemplos para a predição da nota total da prova objetiva do ENEM 2022:

  • Aluno 1) Mulher, negra, estuda em uma escola pública, a mãe não concluiu a graduação e possui uma renda familiar inferior a R$ 2.424,00.

  • Aluno 2) Homem, branco, estuda em uma escola privada, a mãe concluiu a graduação e possui uma renda familiar superior a R$ 9.696,01.

A suposta nota do aluno 1 seria 2.046 enquanto que a nota do aluno 2 seria 2.505.

Footnotes

  1. R Core Team (2023). _R: A Language and Environment for Statistical Computing_. R Foundation for Statistical Computing, Vienna, Austria. <https://www.R-project.org/>.↩︎