Introdução

O presente trabalho tem por objetivo analisar os fatores que podem determinar a nota média de uma escola no Exame Nacional de Ensino Médio (ENEM). Para isso, utilizamos os microdados por escola do exame de 2015 cruzado com os microdados por escola do censo escolar do mesmo ano. Utilizamos o modelo de regressão linear múltipla com erros normais e variância constante para analisar a relação entre a variável resposta (média das notas das provas objetivos) e os regressores selecionados e testamos a adequacidade com as hipóteses iniciais e precisão do modelo.

Metodologia

Tratamento de dados e ajuste

Para este trabalho, utilizamos o banco de dados dos microdados do enem por escola de 2005 a 2015 (disponível aqui), selecionando os dados do ano de 2015 (os mais recentes disponíveis). Os dados faltantes foram eliminados.
Para a construção do modelo, selecionamos as seguintes variáveis:

  • A variável contínua Nota média da prova objetiva (NU_MEDIA_OBJ);
  • A variável contínua limitada de \(0\) a \(100\) Taxa de Adequação Docente (PC_FORMACAO_DOCENTE), isto é, o percentual de docentes da escola que lecionam disciplinas nos quais são licenciados;
  • A variável contínua limitada de \(0\) a \(100\) Taxa de Reprovação (NU_TAXA_REPROVACAO);
  • A variável contínua limitada de \(0\) a \(100\) Taxa de Permanência (NU_TAXA_PERMANENCIA), isto é, o percentual de alunos que fizeram todo o ensino médio na mesma escola;
  • A variável categórica INSE (INSE), que representa o grupo socioeconômico no qual se encaixa o corpo discente da escola. No qual o grupo 1 é o grupo de maior vulnerabilidade socioeconômica e o grupo 6 o de menor vulnerabilidade e maior poder aquisitivo. Mais detalhes aqui.
  • A variável categórica Localização da Escola (TP_LOCALIZACAO_ESCOLA), isto é, se a escola se encontra no meio urbano ou no meio rural;
  • A variável categórica Dependencia Administrativa (TP_DEPENDENCIA_ADM_ESCOLA), isto é, a qual ente público (união, estado ou município) a escola é subordinada ou se é uma escola de administração privada. Como uma de nossas variáveis é a média aritmética das notas daquela escola, escolhemos apenas as escolas em que houve mais de 30 participantes no Enem, para garantir a representatividade da média das notas como estimativa da média real da escola. Construiremos então o modelo de regressão pretendido. Faremos o modelo linear: \[Y_{i}=\beta _{0}+\beta _{1}x_{1i}+\beta _{2}x_{2i}+\beta _{3}x_{3i}+\beta _{4}l_{i}+\beta _{5}g_{1i}+\beta _{6}g_{2i}+\beta _{7}g_{3i}+\beta _{8}g_{4i}+\beta _{9}g_{5i}+\beta _{10}a_{fi}+\beta _{11}a_{ei}+\beta _{12}a_{mi}+\varepsilon _{i}\]

No qual:

  • \(Y_{i}\) é a variável resposta Nota média da prova objetiva;
  • \(x_{1i}\) é o regressor Taxa de Inadequação Docente, isto é o percentual de professores da escola que ensinam disciplinas no qual não são licenciados;
  • \(x_{2i}\) é o regressor Taxa de Reprovação;
  • \(x_{3i}\) é o regressor Taxa de Não-Permanência, isto é, o percentual de alunos que não fizeram todo o ensino médio na escola em que o concluiram;
  • \(l_{i}\) é a variável dummy Localização, que vale 1 se a escola se localizar na área rural e 0 na área urbana;
  • \(g_{1i}\) é a variável dummy INSE no grupo 1, que vale 1 se a escola pertence ao grupo 1 e 0 caso contrário;
  • \(g_{2i}\) é a variável dummy INSE no grupo 2, que vale 1 se a escola pertence ao grupo 2 e 0 caso contrário;
  • \(g_{3i}\) é a variável dummy INSE no grupo 3, que vale 1 se a escola pertence ao grupo 3 e 0 caso contrário;
  • \(g_{4i}\) é a variável dummy INSE no grupo 4, que vale 1 se a escola pertence ao grupo 4 e 0 caso contrário;
  • \(g_{5i}\) é a variável dummy INSE no grupo 5, que vale 1 se a escola pertence ao grupo 5 e 0 caso contrário;
  • \(a_{fi}\) é a variável dummy Administração federal, que vale 1 se a escola é federal e 0 caso contrário;
  • \(a_{ei}\) é a variável dummy Administração estadual, que vale 1 se a escola é estadual e 0 caso contrário;
  • \(a_{mi}\) é a variável dummy Administração municipal, que vale 1 se a escola é municipal e 0 caso contrário;
  • \(\varepsilon _{i}\) representa o erro.

Ao colocar o grupo 6 da variável INSE, a localização na zona urbana e administração escolar privada como defaults dos respectivos conjuntos de variáveis dummy, e pela própria construção do modelo, desejamos analisar a influência de cada variável no intercepto, isto é, o efeito marginal de cada variável no resultado da nota média da escola.
Inicialmente, suporemos o erro com distribuição \(\varepsilon _{i}\sim N(0,\sigma^{2} )\ \forall i\), isto é, erros com distribuição normal de média 0 e variância constante. Além disso suporemos os erros sendo dois-a-dois não correlacionados, isto é \(cov(\varepsilon _{i},\varepsilon _{j})=0\ \forall i,j\).

No R Studio, teremos no tratamento da base de dados:

library(readr)
MICRODADOS_ENEM_ESCOLA <- read.csv("~/enem/DADOS/MICRODADOS_ENEM_ESCOLA.csv", sep=";")

dados <- (MICRODADOS_ENEM_ESCOLA)
dados <- dados[dados$NU_ANO == 2015 & dados$NU_PARTICIPANTES >30,]
dados$NU_MEDIA_OBJ <- (dados$NU_MEDIA_CH + dados$NU_MEDIA_CN + dados$NU_MEDIA_LP + dados$NU_MEDIA_MT)/4
dados$NU_MEDIA_TOT <- (dados$NU_MEDIA_CH + dados$NU_MEDIA_CN + dados$NU_MEDIA_LP + dados$NU_MEDIA_MT + dados$NU_MEDIA_RED)/5
dados <- na.omit(dados)
dados <- dados[dados$NU_TAXA_PERMANENCIA >= 0 & dados$NU_TAXA_PERMANENCIA <= 100,]

Agora construindo as variáveis:

y = as.double(dados$NU_MEDIA_OBJ)
x1 = as.double(100 - dados$PC_FORMACAO_DOCENTE) 
x2 = as.double(dados$NU_TAXA_REPROVACAO) 
x3 =  as.double((100 - dados$NU_TAXA_PERMANENCIA)) 
l = as.numeric(dados$TP_LOCALIZACAO_ESCOLA == 2)  ##dummy 
g1 = as.numeric(dados$INSE == 'Grupo 1') #dummy
g2 = as.numeric(dados$INSE == 'Grupo 2') #dummy
g3 = as.numeric(dados$INSE == 'Grupo 3') #dummy
g4 = as.numeric(dados$INSE == 'Grupo 4') #dummy
g5 = as.numeric(dados$INSE == 'Grupo 5') #dummy
af = as.numeric(dados$TP_DEPENDENCIA_ADM_ESCOLA == 1) #dummy
ae = as.numeric(dados$TP_DEPENDENCIA_ADM_ESCOLA == 2) #dummy
am = as.numeric(dados$TP_DEPENDENCIA_ADM_ESCOLA == 3) #dummy

base <- as.data.frame(cbind(y,x1,x2,x3,l,g1,g2,g3,g4,g5,af,ae,am))

Ajustando agora o modelo:

fit = lm(y ~ x1 + x2 + x3 + l + g1 + g2 + g3 + g4 + g4 + g5 + af + ae + am, base)

summary(fit)
## 
## Call:
## lm(formula = y ~ x1 + x2 + x3 + l + g1 + g2 + g3 + g4 + g4 + 
##     g5 + af + ae + am, data = base)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -86.420 -13.886  -2.112  10.838 186.102 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  622.09645    0.98157 633.777  < 2e-16 ***
## x1            -0.18974    0.01530 -12.402  < 2e-16 ***
## x2            -0.11448    0.03207  -3.569 0.000359 ***
## x3            -0.17108    0.01387 -12.333  < 2e-16 ***
## l             -5.74795    1.63554  -3.514 0.000443 ***
## g1          -110.13800    1.40761 -78.244  < 2e-16 ***
## g2          -101.16446    1.32187 -76.531  < 2e-16 ***
## g3           -88.93046    1.13394 -78.426  < 2e-16 ***
## g4           -72.23101    1.09635 -65.883  < 2e-16 ***
## g5           -43.95537    1.00131 -43.898  < 2e-16 ***
## af            23.83086    1.58380  15.047  < 2e-16 ***
## ae           -35.38710    0.77297 -45.781  < 2e-16 ***
## am           -26.76981    2.76420  -9.684  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 23.01 on 9742 degrees of freedom
## Multiple R-squared:  0.7834, Adjusted R-squared:  0.7832 
## F-statistic:  2937 on 12 and 9742 DF,  p-value: < 2.2e-16

Obtemos então aproximadamente o modelo preliminar \(\widehat{Y_{i}}=618.62-0.19x_{1i}-0.14x_{2i} -0.031x_{3i}-5.49l_{i} -110.14g_{1i}-101.54g_{2i} -89.3g_{3i}-72.47g_{4i}-44.14g_{5i}+26.28a_{fi}-34.38a_{ei} -25.22a_{mi}\) com \(R^{2}=0,7811\), isto é, aproximadamente 78,1% da variabilidade de \(Y\) pode ser explicada pela variabilidade dos regressores. Á primeira vista o modelo parece ter boa adequação aos dados.

Seleção de variáveis

Na seleção dos modelos, utilizamos o critério de seleção backward, em que ajustando-se inicialmente o modelo com todas as variáveis disponíveis, vamos eliminando as variáveis menos significativas e reajustando o modelo sem as variáveis excluídas em cada passo. No nosso caso, temos:

 step(fit, trace = FALSE, direction = 'backward')
## 
## Call:
## lm(formula = y ~ x1 + x2 + x3 + l + g1 + g2 + g3 + g4 + g4 + 
##     g5 + af + ae + am, data = base)
## 
## Coefficients:
## (Intercept)           x1           x2           x3            l  
##    622.0965      -0.1897      -0.1145      -0.1711      -5.7480  
##          g1           g2           g3           g4           g5  
##   -110.1380    -101.1645     -88.9305     -72.2310     -43.9554  
##          af           ae           am  
##     23.8309     -35.3871     -26.7698

Portanto, pela seleçao através do método backward, nenhuma variável original do problema foi eliminada.

Diagnóstico

Normalidade

Uma das primeiras hipóteses admitidas foi a da normalidade dos erros. Analisando graficamente os resíduos do modelo, temos:

Analisando também o gráfico normal com envelope: 

## 
## Banda de  90 % de confianca, obtida por  500  simulacoes.

Nota-se que os resíduos apresentam uma curtose muito maior que o de uma distribuição normal típica. Além de várias observações se distanciarem dos quantis teóricos, fora do envelope. Utilizando o teste de Anderson-Darling para testar a normalidade dos resíduos (sendo a hipótese nula a normalidade).

## 
##  Anderson-Darling normality test
## 
## data:  residuals(fit)
## A = 92.579, p-value < 2.2e-16

Logo, ao nível \(\alpha = 0,05\), rejeitamos a hipotése da normalidade dos resíduos. O modelo assim viola uma de suas hipóteses básicas.

Homocedasticidade

Homocedasticidade é o termo para designar variância constante dos erros. Caso a suposição de homocedasticidade não seja válida, os erros padrões dos estimadores, obtidos pelo Método dos Mínimos Quadrados, são incorretos e portanto a inferência estatística pode não ser válida.
Vale ressaltar que a ausência de homoscedasticidade é chamada de heteroscedasticidade. Com isso, testamos as hipóteses: \[\left\{\begin{matrix} H_{0}:\sigma_{1} ^{2}=...=\sigma_{n} ^{2}\\ H_{1}:\exists i,j;\ \sigma _{i}^{2}\neq \sigma _{j}^{2}\ \end{matrix}\right.\]

Isto é, existe pelo menos um valor de variância diferente das demais.
O gráfico dos resíduos versus valores ajustados (valores preditos) é uma das principais técnicas utilizadas para verificar as suposições dos resíduos. Para o diagnóstico de heteroscedasticidade, tentamos encontrar alguma tendência no gráfico. Por isso, se os pontos estão aleatoriamente distribuídos em torno do 0, sem nenhum comportamento ou tendência, temos indícios de que a variância dos resíduos é homoscedástica.
Analisando o nosso exemplo:

Como Observado, os pontos não estão aleatoriamente distribuídos em torno do 0, seguindo um tipo de tendência. Mostrando um comportamento ou tendência dos resíduos e sugerindo a heterocedasticidade do modelo. Como temos uma amostra suficientemente grande, para corroborar ou rejeitar nossa impressão, vamos proceder com o teste de Goldfeld-Quandt:

## 
##  Goldfeld-Quandt test
## 
## data:  fit
## GQ = 1.265, df1 = 4865, df2 = 4864, p-value < 2.2e-16
## alternative hypothesis: variance increases from segment 1 to 2

E portanto, ao nível \(\alpha = 0,05\), rejeitamos a hipotése da homocedasticidade. O modelo assim viola outra de suas hipóteses básicas.

Independência

Para o diagnóstico de independência também se recorre à análise gráfica e a testes. Uma das análises gráficas para verificar a independência dos resíduos é a feita pela plotagem dos resíduos versus a ordem da coleta dos dados. Se ao avaliar o gráfico, percebemos uma tendência dos pontos, ou seja, se os pontos tiverem um comportamento periódico ou que se repetem em determinadas áreas, temos indícios de dependência dos resíduos. No nosso caso:

Observa-se que os resíduos sugerem um comportamento periódico em relação ao hipótese. Nos dando indícios de que os resíduos são correlatados. Para verificar a hipótese da correlação dos resíduos, usamos o teste de Durbin-Watson, em que testamos se a correlação dos erros no primeiro passo é nula sob a hipótese \(H_{0}\) e maior que zero caso contrário.

## 
##  Durbin-Watson test
## 
## data:  fit
## DW = 1.5208, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0

E portanto, ao nível \(\alpha = 0,05\), rejeitamos a hipotése da independência dos erros.

Multicolinearidade

Preliminarmente, vamos analisar se há alguma forte correlação entre as variáveis explicativas contínuas do modelo.

##            x1         x2         x3
## x1 1.00000000 0.01168165 0.00143391
## x2 0.01168165 1.00000000 0.02894702
## x3 0.00143391 0.02894702 1.00000000

Nenhum par de variáveis apresenta uma correção maior que \(0,5\) (em módulo), além de visualmente não ser evidente que há alguma relação entre elas. Para formalizamos que nenhuma dessas três variáveis apresenta o problema de multicolinearidade, usaremos o conceito de Variance Inflation Factor (VIF), o qual não admitiremos multicolinearidade se algum valor for maior que 10 para alguma variável.

library(faraway)

vif(p)
##       x1       x2       x3 
## 1.000138 1.000974 1.000840

Ou seja, não há problema de multicolinearidade entre as três variáveis em questão.

Medidas de Influência

Inicialmente através de uma análise gráfica, vamos verificar a presença e a quantidade de possíveis outliers no modelo. Temos então:

No qual se observa uma grande quantidade de outliers e pontos influentes, o que provavelmente faz com que o modelo viole as suas hipóteses iniciais. Avaliando as observações pelo conceito da Distância de Cook, temos:

No qual novamente se constata a grande quantidade de outliers. Não se torna viável eliminar todos para um novo ajuste. O que nos sugere buscar hipóteses diferentes para a distribuição dos erros.

Conclusão

O modelo construído sob as suposições de normalidade dos erros, homocedasticidade e ausência de autocorrelação não se mostrou adequado aos dados disponíveis. Necessitando de outras hipóteses para a adequada modelagem.

Bibliografia