Abstract
Aula 3 da disciplina de Estatistica Aplicada do curso de Especializacao em Big Data, Data Science e Data Analytics. Esta aula contem um exemplo da analise de regressao multipla onde busca-se avaliar uma possivel dependencia de Y em relacao a variaveis independentes (preditoras), considerando a populacao do estudo, com base em dados amostrais. Alem disso, expressa-se matematicamente esta relacao por meio de uma equacao, bem como se faz previsoes (estimativas) da variavel Y com base em valores das vaiaveis independentes X1, X2, X3, … Xn.library(descr)
library(kableExtra)
library(car)
library(carData)
library(sciplot)
library(tidyverse)
library(dplyr)
library(ggplot2)
library(GGally)
library(plotly)
library(ExpDes.pt)
library(DT)
library(readxl)
library(xlsx)
library(relaimpo)
library(hrbrthemes)
Etapas principais para a modelagem de regressao:
Definir a variavel dependente;
Selecionar variaveis independentes;
Coleta dos dados pareados das respectivas variaveis;
Ajustar o modelo obtendo os coeficientes da equacao e valor-p;
Reavaliacao das variaveis que permanecerao no modelo;
Ajustar o modelo novamente com as variaveis definitivas, obtendo os coeficientes da equacao e valor-p;
Calculo dos residuos para avaliar o ajuste do modelo;
Revisao para redefinicao e possivel transformacao de variaveis para melhoria do ajuste do modelo;
Interpretacao do modelo ajustado (equacao) e calculos de previsao;
Utilizacao da equacao para predicao da base que nao apresenta a variavel resposta;
Exportacao dos dados gerados no modelo.
Uma instituicao de ensino superior deseja avaliar se e possivel identificar fatores associados ao desempenho de seus alunos formandos em uma avaliacao do MEC. As variaveis explicativas testadas foram nota media em disciplinas do curso, idade, sexo e distancia da residencia ao campus (local da prova), para uma amostra aleataria de 127 alunos.
Dados do exemplo1.
bd <- read_excel("C:/LEONARDO/Rpubs/estatistica_aplicada_aula3_leonardo.xlsx")
datatable(
bd, extensions = 'Buttons', options = list(
dom = 'Bfrtip',
buttons = c('copy', 'csv', 'excel', 'pdf', 'print')
)
)
Avaliacao_mec
Como e um exemplo para calculos de previsao, algumas observacoes (alunos) estao sem valores da variavel Avaliacao_mec, portanto, incialmente devemos separar a parte da base de dados com os casos que apresentem essa variavel resposta, formando uma nova base de dados para treino do modelo.
bd_modelo <- bd %>%
na.omit()
bd_teste <- setdiff(bd, bd_modelo)
#opcoes para retirar NA
#bd_modelo<- subset.data.frame(bd, Avaliacao_mec>=0)
#bd_modelo<- subset.data.frame(bd, Avaliacao_mec!="NA")
str(bd_modelo)
## tibble [127 x 14] (S3: tbl_df/tbl/data.frame)
## $ X1_CURSO : chr [1:127] "Administração" "Administração" "Administração" "Administração" ...
## $ X2_DISC : chr [1:127] "Planejamento contábil" "Planejamento contábil" "Planejamento contábil" "Planejamento contábil" ...
## $ X3_EXIG_MAT : num [1:127] 2 2 2 2 2 2 2 2 2 2 ...
## $ X4_MODAL : chr [1:127] "1" "1" "1" "1" ...
## $ X5_NOTAS_ANT : num [1:127] 5.8 6.4 9 9.2 7.3 9.3 7.5 7 9.4 6.4 ...
## $ X6_NUM_EVASAO: num [1:127] 3 2 0 0 0 0 1 2 0 2 ...
## $ X7_NUM_ATRAS : num [1:127] 0 0 1 0 2 1 0 0 0 0 ...
## $ X8_SEXO : num [1:127] 0 0 0 0 1 1 1 0 1 0 ...
## $ X9_IDADE : num [1:127] 30 27 23 25 30 28 25 25 23 25 ...
## $ X10_DIST : num [1:127] 1 3 2 2 3 1 2 1 1 1 ...
## $ X11_SEMESTRE : chr [1:127] "1" "1" "1" "1" ...
## $ X12_TAM_TURMA: num [1:127] 25 30 30 33 30 43 25 35 42 20 ...
## $ Avaliacao_mec: num [1:127] 62 71 87 93.2 82.4 81 74.6 71.8 68.3 72.2 ...
## $ Y_SITUACAO : num [1:127] 0 0 0 0 0 0 0 0 0 0 ...
## - attr(*, "na.action")= 'omit' Named int [1:243] 13 14 15 16 17 18 19 20 21 22 ...
## ..- attr(*, "names")= chr [1:243] "13" "14" "15" "16" ...
Visualizacao de um histograma da variavel resposta Avaliacao_mec na base de dados de modelagem.
hist(bd_modelo$Avaliacao_mec, col="gray", las=1, xlab="Avaliacao_mec", ylab="Frequência",
main="histograma", cex.lab=1.25, cex.axis=1)
abline(v=mean(bd_modelo$Avaliacao_mec), col="red", lty=3)
X5_NOTAS_ANT; X9_IDADE; X8_SEXO; X10_DIST
par(mfrow=c(2,2))
#X5_NOTAS_ANT
plot(bd_modelo$Avaliacao_mec,bd_modelo$X5_NOTAS_ANT, pch=16, col="black", ylab="NOTAS_ANT", xlab="Avaliacao_mec", las=1)
#X9_IDADE
plot(bd_modelo$Avaliacao_mec,bd_modelo$X9_IDADE, pch=16, col="black", ylab="IDADE", xlab="Avaliacao_mec", las=1)
#X8_SEXO
plot(bd_modelo$Avaliacao_mec,bd_modelo$X8_SEXO, pch=16, col="black", ylab="SEXO", xlab="Avaliacao_mec", las=1)
#X10_DIST
plot(bd_modelo$Avaliacao_mec,bd_modelo$X10_DIST, pch=16, col="black", ylab="DIST", xlab="Avaliacao_mec", las=1)
#X5_NOTAS_ANT
correlacao.X5_NOTAS_ANT <- cor(bd_modelo$Avaliacao_mec,bd_modelo$X5_NOTAS_ANT)
#X9_IDADE
correlacao.X9_IDADE <- cor(bd_modelo$Avaliacao_mec,bd_modelo$X9_IDADE)
#X8_SEXO
correlacao.X8_SEXO <- cor(bd_modelo$Avaliacao_mec,bd_modelo$X8_SEXO)
#X10_DIST
correlacao.X10_DIST <- cor(bd_modelo$Avaliacao_mec,bd_modelo$X10_DIST)
Os valores de correlacao entre a variavel resposta Avaliacao_mec e cada uma das variaveis explicativas sao:
X5_NOTAS_ANT: 0.4062551.
X9_IDADE: 0.5934118.
X8_SEXO: 0.0138624.
X10_DIST: -0.026586.
bd_modelo.ggpairs <- bd_modelo[, c(5, 8, 9, 10, 13)]
ggpairs(bd_modelo.ggpairs)
resultado.rlm <- lm(Avaliacao_mec ~ X5_NOTAS_ANT + X9_IDADE + X8_SEXO + X10_DIST,
data=bd_modelo)
summary(resultado.rlm)
##
## Call:
## lm(formula = Avaliacao_mec ~ X5_NOTAS_ANT + X9_IDADE + X8_SEXO +
## X10_DIST, data = bd_modelo)
##
## Residuals:
## Min 1Q Median 3Q Max
## -24.3539 -4.4799 -0.8912 3.0237 29.2744
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.47420 6.26408 -0.395 0.694
## X5_NOTAS_ANT 4.42930 0.58765 7.537 9.36e-12 ***
## X9_IDADE 1.72169 0.16592 10.377 < 2e-16 ***
## X8_SEXO -2.70870 1.27075 -2.132 0.035 *
## X10_DIST 0.05424 0.81972 0.066 0.947
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.971 on 122 degrees of freedom
## Multiple R-squared: 0.559, Adjusted R-squared: 0.5446
## F-statistic: 38.67 on 4 and 122 DF, p-value: < 2.2e-16
r2.resultado.rlm <- summary(resultado.rlm)$r.squared
r2.resultado.rlm <- round(r2.resultado.rlm*100, digits = 2)
adj.r2.resultado.rlm <- summary(resultado.rlm)$adj.r.squared
O modelo utilizado com as 4 variaveis explicou (r.squared) 55.9% da variacao dos dados de Avaliacao_mec. O valor de r2 ajustado ficou em 0.5445919.
A seguir vemos a equacao gerada do modelo de regressao linear multiplo com seus coeficientes estao descritos:
equacao.resultado.rlm <- as.formula(
paste0("Avaliacao_meq ~ ", round(coefficients(resultado.rlm)[1],4), " + ",
paste(sprintf("%.2f * %s",
coefficients(resultado.rlm)[-1],
names(coefficients(resultado.rlm)[-1])),
collapse=" + ")
)
)
equacao.resultado.rlm
## Avaliacao_meq ~ -2.4742 + 4.43 * X5_NOTAS_ANT + 1.72 * X9_IDADE +
## -2.71 * X8_SEXO + 0.05 * X10_DIST
Os coeficientes de cada uma das variaveis explicativas utilizadas no modelo estao destacadas a seguir.
resultado.rlm$coefficients
## (Intercept) X5_NOTAS_ANT X9_IDADE X8_SEXO X10_DIST
## -2.47419518 4.42930473 1.72169305 -2.70869542 0.05424299
p.value.coef.X5_NOTAS_ANT <- coef(summary(resultado.rlm))["X5_NOTAS_ANT","Pr(>|t|)"]
p.value.coef.X9_IDADE <- coef(summary(resultado.rlm))["X9_IDADE","Pr(>|t|)"]
p.value.coef.X8_SEXO <- coef(summary(resultado.rlm))["X8_SEXO","Pr(>|t|)"]
p.value.coef.X10_DIST <- coef(summary(resultado.rlm))["X10_DIST","Pr(>|t|)"]
O valor-p para o coeficiente de X5_NOTAS_ANT resultou em: 9.364891610^{-12}.
O valor-p para o coeficiente de X9_IDADE resultou em:1.813273610^{-18}.
O valor-p para o coeficiente de X8_SEXO resultou em:0.0350484.
O valor-p para o coeficiente de X10_DIST resultou em:0.9473487.
Para avaliacao do modelo utilizamos a verificacao dos valores de p das variaveis explicativas. Observa-se que a variavel X10_DIST apresentou 0.9473487 nao significativo, portanto sera retirada do modelo.
Busca-se tambem fazer uma visualizacao geral do modelo e verificar a importancia de cada uma das variaveis explicativas (X1, X2, X3, … , Xn) sobre a variavel resposta (Y).
avaliacao.resultado.rlm <- calc.relimp(resultado.rlm, type = "lmg")
avaliacao.resultado.rlm
## Response variable: Avaliacao_mec
## Total response variance: 106.6941
## Analysis based on 127 observations
##
## 4 Regressors:
## X5_NOTAS_ANT X9_IDADE X8_SEXO X10_DIST
## Proportion of variance explained by model: 55.9%
## Metrics are not normalized (rela=FALSE).
##
## Relative importance metrics:
##
## lmg
## X5_NOTAS_ANT 0.1833746094
## X9_IDADE 0.3688272024
## X8_SEXO 0.0064981174
## X10_DIST 0.0003493333
##
## Average coefficients for different model sizes:
##
## 1X 2Xs 3Xs 4Xs
## X5_NOTAS_ANT 3.8883631 4.0268083 4.2065801 4.42930473
## X9_IDADE 1.6255199 1.6486793 1.6806107 1.72169305
## X8_SEXO 0.2859671 -0.6292728 -1.6268172 -2.70869542
## X10_DIST -0.3619761 -0.2482600 -0.1128168 0.05424299
#dataframe com os resultados da avaliacao.resultado.rlm
impo.variaveis <- data.frame(avaliacao.resultado.rlm$lmg, avaliacao.resultado.rlm$lmg.rank, avaliacao.resultado.rlm$namen[-1])
#rank de $lmg
impo.variaveis <- arrange(impo.variaveis, desc(avaliacao.resultado.rlm$lmg))
impo.variaveis$p.values <- c(round(p.value.coef.X5_NOTAS_ANT, digits = 5),
round(p.value.coef.X9_IDADE, digits = 5),
round(p.value.coef.X8_SEXO, digits = 5),
round(p.value.coef.X10_DIST, digits = 5))
datatable(head(impo.variaveis), colnames = c('lmg', 'rank', 'Variaveis', "valor-p"),
caption = 'Tabela 1: Rank de importancia das variedades explicativas no modelo de regressao linear multipla.') %>%
formatStyle('p.values',
backgroundColor = styleInterval(0.05, c('green', 'red'))
)
As variaveis que permanecerao no modelo sao: X5_NOTAS_ANT, X9_IDADE e X8_SEXO.
Retiramos do modelo X10_DIST.
resultado.rlm.ajust <- lm(Avaliacao_mec ~ X5_NOTAS_ANT + X9_IDADE + X8_SEXO, data=bd_modelo)
summary(resultado.rlm.ajust)
##
## Call:
## lm(formula = Avaliacao_mec ~ X5_NOTAS_ANT + X9_IDADE + X8_SEXO,
## data = bd_modelo)
##
## Residuals:
## Min 1Q Median 3Q Max
## -24.3283 -4.4596 -0.9163 2.9915 29.2384
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.3745 6.0555 -0.392 0.6956
## X5_NOTAS_ANT 4.4290 0.5852 7.568 7.72e-12 ***
## X9_IDADE 1.7212 0.1651 10.427 < 2e-16 ***
## X8_SEXO -2.7060 1.2649 -2.139 0.0344 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.942 on 123 degrees of freedom
## Multiple R-squared: 0.559, Adjusted R-squared: 0.5483
## F-statistic: 51.98 on 3 and 123 DF, p-value: < 2.2e-16
r2.resultado.rlm.ajust <- summary(resultado.rlm.ajust)$r.squared
r2.resultado.rlm.ajust <- round(r2.resultado.rlm.ajust*100, digits = 2)
adj.r2.resultado.rlm.ajust <- summary(resultado.rlm.ajust)$adj.r.squared
O modelo utilizado com as 3 variaveis explicou (r.squared) 55.9% da variacao dos dados de Avaliacao_mec. O valor de r2 ajustado ficou em 0.5482782.
A seguir vemos a equacao gerada do modelo de regressao linear multiplo com seus coeficientes estao descritos:
equacao.resultado.rlm.ajust <- as.formula(
paste0("Avaliacao_meq ~ ", round(coefficients(resultado.rlm.ajust)[1],3), " + ",
paste(sprintf("%.2f * %s",
coefficients(resultado.rlm.ajust)[-1],
names(coefficients(resultado.rlm.ajust)[-1])),
collapse=" + ")
)
)
equacao.resultado.rlm.ajust
## Avaliacao_meq ~ -2.374 + 4.43 * X5_NOTAS_ANT + 1.72 * X9_IDADE +
## -2.71 * X8_SEXO
Avaliacao da equacao final, com as variaveis explicativas ajustadas.
avaliacao.resultado.rlm.ajust <- calc.relimp(resultado.rlm.ajust, type = "lmg")
avaliacao.resultado.rlm.ajust
## Response variable: Avaliacao_mec
## Total response variance: 106.6941
## Analysis based on 127 observations
##
## 3 Regressors:
## X5_NOTAS_ANT X9_IDADE X8_SEXO
## Proportion of variance explained by model: 55.9%
## Metrics are not normalized (rela=FALSE).
##
## Relative importance metrics:
##
## lmg
## X5_NOTAS_ANT 0.183382480
## X9_IDADE 0.369150704
## X8_SEXO 0.006500253
##
## Average coefficients for different model sizes:
##
## 1X 2Xs 3Xs
## X5_NOTAS_ANT 3.8883631 4.096111 4.428951
## X9_IDADE 1.6255199 1.660340 1.721192
## X8_SEXO 0.2859671 -1.094494 -2.705966
#dataframe com os resultados da avaliacao.resultado.rlm
impo.variaveis.ajust <- data.frame(avaliacao.resultado.rlm.ajust$lmg, avaliacao.resultado.rlm.ajust$lmg.rank, avaliacao.resultado.rlm.ajust$namen[-1])
#rank de $lmg
impo.variaveis.ajust <- arrange(impo.variaveis.ajust, desc(avaliacao.resultado.rlm.ajust$lmg))
p.value.coef.X5_NOTAS_ANT.ajust <- coef(summary(resultado.rlm.ajust))["X5_NOTAS_ANT","Pr(>|t|)"]
p.value.coef.X9_IDADE.ajust <- coef(summary(resultado.rlm.ajust))["X9_IDADE","Pr(>|t|)"]
p.value.coef.X8_SEXO.ajust <- coef(summary(resultado.rlm.ajust))["X8_SEXO","Pr(>|t|)"]
impo.variaveis.ajust$p.values <- c(round(p.value.coef.X5_NOTAS_ANT.ajust, digits = 5),
round(p.value.coef.X9_IDADE.ajust, digits = 5),
round(p.value.coef.X8_SEXO.ajust, digits = 5))
datatable(head(impo.variaveis.ajust), colnames = c('lmg', 'rank', 'Variaveis', "valor-p"),
caption = 'Tabela 2: Rank de importancia das variedades explicativas no modelo de regressao linear multipla ajustado.') %>%
formatStyle('p.values',
backgroundColor = styleInterval(0.05, c('green', 'red'))
)
par(mfrow=c(2,2))
plot(fitted(resultado.rlm.ajust),residuals(resultado.rlm.ajust),xlab="Valores Ajustados",ylab="Residuos", pch=16)
abline(h=0, col="red")
plot(bd_modelo$X5_NOTAS_ANT,residuals(resultado.rlm.ajust),xlab="X5_NOTAS_ANT",ylab="Residuos", pch=16)
abline(h=0, col="red")
plot(bd_modelo$X9_IDADE,residuals(resultado.rlm.ajust),xlab="X9_IDADE",ylab="Residuos", pch=16)
abline(h=0, col="red")
plot(bd_modelo$X8_SEXO,residuals(resultado.rlm.ajust),xlab="X8_SEXO",ylab="Residuos", pch=16)
abline(h=0, col="red")
Avaliacao da suposicao de Normalidade dos erros.
par(mfrow=c(1,2))
qqnorm(resultado.rlm.ajust$residuals, ylab="Residuos",xlab="Quantis teoricos",main="")
qqline(resultado.rlm.ajust$residuals)
hist(resultado.rlm.ajust$residuals, main="histogramas residuos")
shapiro.test(resultado.rlm.ajust$residuals)
##
## Shapiro-Wilk normality test
##
## data: resultado.rlm.ajust$residuals
## W = 0.93053, p-value = 6.01e-06
O modelo indicou que 55.9% da variacao da nota dos alunos junto ao MEC foi explicada pela variacao da idade, sexo e media das notas finais de disciplinas conforme a equacao:
equacao.resultado.rlm.ajust
## Avaliacao_meq ~ -2.374 + 4.43 * X5_NOTAS_ANT + 1.72 * X9_IDADE +
## -2.71 * X8_SEXO
Intepretacao dos coeficientes B do modelo:
coef.resultado.rlm.ajust.X5_NOTAS_ANT <- resultado.rlm.ajust$coefficients[2]
coef.resultado.rlm.ajust.X9_IDADE <- resultado.rlm.ajust$coefficients[3]
coef.resultado.rlm.ajust.X8_SEXO <- resultado.rlm.ajust$coefficients[4]
X1 (X5_NOTAS_ANT): A cada um ponto a mais na nota media de disciplinas concluidas pelo aluno, observou-se em media, 4.4289509 pontos a mais na avaliacao do MEC.
X2 (X9_IDADE): A cada um ano a mais de idade do aluno, observou-se em media,1.7211919 pontos a mais na avaliacao do MEC.
X3 (X8_SEXO): Em media, as mulheres tem a nota na avaliacao do MEC, -2.7059656 pontos superiores, em relacao aos homens.
Formacao da base de dados com as variaveis utilizadas na equacao.
bd_new <- subset(bd, select = c('X5_NOTAS_ANT','X9_IDADE','X8_SEXO'))
bd_nova <- bd %>%
dplyr::select("X5_NOTAS_ANT" | 'X9_IDADE' | 'X8_SEXO' )
Aplicar o modelo de predicao de valores na base de dados.
predicao <- predict(resultado.rlm.ajust, newdata=bd, se.fit = TRUE, type="response")
bd <- data.frame(bd,predicao)
Avaliacao dos Residuos (Avaliacao_meq - fit) na aplicacao do modelo.
bd$residuos <- bd$Avaliacao_mec - bd$fit
resultado.modelo <-lm(bd$fit~bd$Avaliacao_mec)
summary(resultado.modelo)
##
## Call:
## lm(formula = bd$fit ~ bd$Avaliacao_mec)
##
## Residuals:
## Min 1Q Median 3Q Max
## -18.4559 -3.4715 -0.0137 3.0987 16.1492
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 31.10932 3.16608 9.826 <2e-16 ***
## bd$Avaliacao_mec 0.55903 0.04441 12.588 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.149 on 125 degrees of freedom
## (243 observations deleted due to missingness)
## Multiple R-squared: 0.559, Adjusted R-squared: 0.5555
## F-statistic: 158.5 on 1 and 125 DF, p-value: < 2.2e-16
plot.resultado.modelo <- ggplot(bd, aes(x=bd$fit, y=bd$Avaliacao_mec)) +
geom_point() +
geom_smooth(method=lm , color="red", fill="#69b3a2", se=TRUE) +
theme_ipsum()
plot.resultado.modelo + ggtitle("Dispersao de Avaliacao_mec e valores do modelo") +
xlab("Predicao do modelo") + ylab("Avaliacao_mec")
Erro de estimativa percentual (EEP) de cada observacao da previsao do modelo sobre o valor real da variavel resposta (Avaliacao_mec).
bd$erro_estimativa_relativo <- (bd$fit-bd$Avaliacao_mec)/bd$Avaliacao_mec
bd$erro_estimativa_percentual <- round(((bd$fit-bd$Avaliacao_mec)/bd$Avaliacao_mec)*100, digits = 2)
par(mfrow=c(1,2))
hist(bd$erro_estimativa_percentual , border=F , col=rgb(0.1,0.8,0.3,0.5) , xlab="distribuicao do erro de estimativa percentual" , main="Histograma de EEP")
boxplot(bd$erro_estimativa_percentual , xlab="erro de estimativa percentual" , col=rgb(0.8,0.8,0.3,0.5) , las=2, main="Boxplot de EEP")
Calculo do Erro Relativo Medio das Previsoes (MAPE mean absolute percentage error):
Esta e uma medida de acuracia de predicao do modelo estatistico, quanto menor o valor, melhor este indice.
MAPE <- mean(abs(bd$fit-bd$Avaliacao_mec)/bd$Avaliacao_mec,na.rm = TRUE)*100
O erro relativo medio das previsoes do modelo resultou em 6.967906% (MAPE).
bd_final <- write.xlsx(bd, file="estatistica_aplicada_aula3_leonardo_exemplo_finalizado.xls")