Modelos Hierárquicos - Lab 1
Modelos Hierárquicos - Lab 1
- 1 Preparação
- 2 Analise inicial
- 3 Analise grafica inicial
- 4 Correlaçõees
- 5 Outras visualizações
- 6 Ajuste
- 7 Verificação das hipoteses
- 8 Análise de nível
- 9 Criando varíável artifical de nível 2
- 10 Modelo de intercepto aleatório
- 11 Modelo com interação
- 12 Modelo de coeficiente aleatório
- 13 Verificando se a proporção de homens é significante
- 14 Modelo final
- 15 Análise de resíduos
- 16 Equação
- 17 Dicionário de variáveis
- 18 Interpretação
1 Preparação
load("C:/Users/raphael.fernandez/Downloads/ArquivoProvaBR2007trab.RData")
data=dadosap
attach(data)2 Analise inicial
str(data)## 'data.frame': 26097 obs. of 29 variables:
## $ id_aluno : chr " 679298" " 679298" " 679298" " 679299" ...
## $ id_turma : chr "RJ01697" "RJ01697" "RJ01697" "RJ01697" ...
## $ id_serie : chr "8" "8" "8" "8" ...
## $ pk_cod_entidade : chr "3307633" "3307633" "3307633" "3307633" ...
## $ id_dependencia_adm: chr "3" "3" "3" "3" ...
## $ id_localizacao : chr "1" "1" "1" "1" ...
## $ sigla_uf : chr "RJ" "RJ" "RJ" "RJ" ...
## $ cod_uf : chr "33" "33" "33" "33" ...
## $ no_municipio : chr "RIO DE JANEIRO" "RIO DE JANEIRO" "RIO DE JANEIRO" "RIO DE JANEIRO" ...
## $ cod_municipio : chr "3304557" "3304557" "3304557" "3304557" ...
## $ y1 : num 150 260 172 155 140 ...
## $ y2 : num 194 230 174 152 223 ...
## $ AcertosLP : int 4 14 6 8 2 4 4 5 4 9 ...
## $ AcertosM : int 4 8 4 5 9 9 6 3 4 4 ...
## $ N_alu_pop : int 35 35 35 35 35 35 35 35 35 35 ...
## $ qst1sexo : Factor w/ 2 levels "Feminin","Masculi": 2 1 2 1 1 1 1 1 1 2 ...
## $ qst2raca : Factor w/ 5 levels "AM","Br","IN",..: 5 5 2 3 3 5 5 5 2 2 ...
## $ qst17numpess : Factor w/ 6 levels ">8","m1","m2",..: 5 2 4 6 5 6 5 5 5 5 ...
## $ qst18moramae : Factor w/ 3 levels "Não","Não outra",..: 3 3 3 3 1 3 3 3 3 3 ...
## $ qst19educmae : Factor w/ 6 levels "11a14","15m",..: 3 3 1 3 6 3 1 5 3 6 ...
## $ qst22morapai : Factor w/ 3 levels "Não","Não outra",..: 1 2 3 1 3 3 3 3 3 3 ...
## $ qst23educpai : Factor w/ 6 levels "11a14","15m",..: 3 6 1 6 6 5 1 4 3 6 ...
## $ qst35trabfora : Factor w/ 2 levels "Não","Sim": 1 1 1 1 1 1 1 1 1 1 ...
## $ qst38reprov : Factor w/ 3 levels "Não","Sim 1vez",..: 1 2 1 1 1 2 2 1 1 1 ...
## $ qst39abando : Factor w/ 3 levels "Não","Sim 1vez",..: 3 1 1 1 1 1 1 1 1 1 ...
## $ r16conspatio : Factor w/ 4 levels "Adequado","Inadequad",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ r36ncompesc : Factor w/ 6 levels "11a15","16a20",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ r52conscomp : Factor w/ 3 levels "Bom","Regular",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ r58percbiblio : Factor w/ 5 levels "Até 25","De26a50",..: 4 4 4 4 4 4 4 4 4 4 ...
## - attr(*, "datalabel")= chr ""
## - attr(*, "time.stamp")= chr " 5 Oct 2017 19:24"
## - attr(*, "formats")= chr "%9s" "%7s" "%1s" "%9s" ...
## - attr(*, "types")= int 7 7 1 7 1 1 2 2 14 7 ...
## - attr(*, "val.labels")= chr "" "" "" "" ...
## - attr(*, "var.labels")= chr "código do aluno na prova brasil" "código da turma na prova brasil" "código da série" "código da escola" ...
## - attr(*, "version")= int 12
## - attr(*, "label.table")=List of 14
## ..$ r58percbiblio: Named int 1 2 3 4 5
## .. ..- attr(*, "names")= chr "Até 25" "De26a50" "De51a75" "Mais75" ...
## ..$ r52conscomp : Named int 1 2 3
## .. ..- attr(*, "names")= chr "Bom" "Regular" "Ruim"
## ..$ r36ncompesc : Named int 1 2 3 4 5 6
## .. ..- attr(*, "names")= chr "11a15" "16a20" "1a5" "21a30" ...
## ..$ r16conspatio : Named int 1 2 3 4
## .. ..- attr(*, "names")= chr "Adequado" "Inadequad" "Inexisten" "Regular"
## ..$ qst39abando : Named int 1 2 3
## .. ..- attr(*, "names")= chr "Não" "Sim 1vez" "Sim duas"
## ..$ qst38reprov : Named int 1 2 3
## .. ..- attr(*, "names")= chr "Não" "Sim 1vez" "Sim duas"
## ..$ qst35trabfora: Named int 1 2
## .. ..- attr(*, "names")= chr "Não" "Sim"
## ..$ qst23educpai : Named int 1 2 3 4 5 6
## .. ..- attr(*, "names")= chr "11a14" "15m" "4a7" "8a10" ...
## ..$ qst22morapai : Named int 1 2 3
## .. ..- attr(*, "names")= chr "Não" "Não outra" "Sim"
## ..$ qst19educmae : Named int 1 2 3 4 5 6
## .. ..- attr(*, "names")= chr "11a14" "15m" "4ª7" "8a10" ...
## ..$ qst18moramae : Named int 1 2 3
## .. ..- attr(*, "names")= chr "Não" "Não outra" "Sim"
## ..$ qst17numpess : Named int 1 2 3 4 5 6
## .. ..- attr(*, "names")= chr ">8" "m1" "m2" "m3" ...
## ..$ qst2raca : Named int 1 2 3 4 5
## .. ..- attr(*, "names")= chr "AM" "Br" "IN" "PR" ...
## ..$ qst1sexo : Named int 1 2
## .. ..- attr(*, "names")= chr "Feminin" "Masculi"summary(data)## id_aluno id_turma id_serie
## Length:26097 Length:26097 Length:26097
## Class :character Class :character Class :character
## Mode :character Mode :character Mode :character
##
##
##
## pk_cod_entidade id_dependencia_adm id_localizacao
## Length:26097 Length:26097 Length:26097
## Class :character Class :character Class :character
## Mode :character Mode :character Mode :character
##
##
##
## sigla_uf cod_uf no_municipio
## Length:26097 Length:26097 Length:26097
## Class :character Class :character Class :character
## Mode :character Mode :character Mode :character
##
##
##
## cod_municipio y1 y2 AcertosLP
## Length:26097 Min. :112.9 Min. :138.7 Min. : 1.00
## Class :character 1st Qu.:205.8 1st Qu.:213.2 1st Qu.: 9.00
## Mode :character Median :238.5 Median :242.8 Median :13.00
## Mean :238.6 Mean :246.5 Mean :12.81
## 3rd Qu.:270.8 3rd Qu.:276.4 3rd Qu.:16.00
## Max. :386.2 Max. :415.1 Max. :26.00
## AcertosM N_alu_pop qst1sexo qst2raca qst17numpess
## Min. : 1.0 Min. : 22 Feminin:13754 AM: 1196 >8 : 500
## 1st Qu.: 7.0 1st Qu.: 60 Masculi:12343 Br: 8106 m1 :1342
## Median :10.0 Median : 91 IN: 1024 m2 :4698
## Mean :10.4 Mean : 97 PR: 4648 m3 :7934
## 3rd Qu.:13.0 3rd Qu.:121 Pa:11123 m4a5:9020
## Max. :26.0 Max. :257 m6a8:2603
## qst18moramae qst19educmae qst22morapai qst23educpai
## Não : 2025 11a14 :7014 Não : 8865 11a14 :6397
## Não outra: 618 15m :1963 Não outra: 1605 15m :2291
## Sim :23454 4ª7 :5925 Sim :15627 4a7 :3988
## 8a10 :5168 8a10 :4330
## ate 3 :1928 ate 3 :1631
## naosabe:4099 naosabe:7460
## qst35trabfora qst38reprov qst39abando r16conspatio
## Não:23321 Não :17756 Não :24984 Adequado :18826
## Sim: 2776 Sim 1vez: 6389 Sim 1vez: 943 Inadequad: 1947
## Sim duas: 1952 Sim duas: 170 Inexisten: 350
## Regular : 4974
##
##
## r36ncompesc r52conscomp r58percbiblio
## 11a15:11304 Bom :21765 Até 25 :2390
## 16a20: 2162 Regular: 3491 De26a50:6500
## 1a5 : 8609 Ruim : 841 De51a75:9411
## 21a30: 568 Mais75 :7103
## 6a10 : 2248 Não há : 693
## m30 : 1206prop.table(table(qst1sexo))## qst1sexo
## Feminin Masculi
## 0.5270338 0.4729662prop.table(table(r52conscomp))## r52conscomp
## Bom Regular Ruim
## 0.83400391 0.13377017 0.032225933 Analise grafica inicial
plot(y1,AcertosLP)
plot(y2,AcertosM)
4 Correlaçõees
cor(cbind(y2,AcertosM,qst1sexo,qst35trabfora))## y2 AcertosM qst1sexo qst35trabfora
## y2 1.00000000 0.93309211 0.09473209 -0.07812509
## AcertosM 0.93309211 1.00000000 0.09494714 -0.07575214
## qst1sexo 0.09473209 0.09494714 1.00000000 0.13742237
## qst35trabfora -0.07812509 -0.07575214 0.13742237 1.00000000tapply(y2,qst1sexo,mean)## Feminin Masculi
## 242.1967 251.21165 Outras visualizações
boxplot(y2~qst1sexo,main="sexo")
boxplot(y2~qst2raca,main="raça")
boxplot(y2~qst17numpess,main="mora com quantas pessoas")
boxplot(y2~qst18moramae,main="mora mãe")
boxplot(y2~qst22morapai,main="mora pai")
boxplot(y2~qst23educpai,main="educa pai")
boxplot(y2~qst35trabfora,main="sexo")
boxplot(y2~qst38reprov,main="reprovação")
boxplot(y2~qst39abando,main="abandonou")
boxplot(y2~r16conspatio,main="patio")
boxplot(y2~r36ncompesc,main="num comp")
boxplot(y2~r52conscomp,main="cons comp")
boxplot(y2~r58percbiblio,main="biblio")
6 Ajuste
detach(data)
table(data$qst39abando)##
## Não Sim 1vez Sim duas
## 24984 943 170data$aban = ifelse(data$qst39abando=="Sim duas",1,0)
data$aban = as.factor(data$aban)
levels(data$aban) = c("Não abandonou mais de uma vez","Abandonou mais de uma vez")
attach(data)
mod = lm(y2 ~ AcertosLP + qst1sexo+qst19educmae+qst38reprov+aban)
summary(mod)##
## Call:
## lm(formula = y2 ~ AcertosLP + qst1sexo + qst19educmae + qst38reprov +
## aban)
##
## Residuals:
## Min 1Q Median 3Q Max
## -150.238 -24.832 -0.196 24.011 142.097
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 172.04323 0.90036 191.082 < 2e-16 ***
## AcertosLP 5.74165 0.05057 113.549 < 2e-16 ***
## qst1sexoMasculi 16.08747 0.46194 34.826 < 2e-16 ***
## qst19educmae15m 11.00958 0.93767 11.741 < 2e-16 ***
## qst19educmae4ª7 -6.34010 0.65194 -9.725 < 2e-16 ***
## qst19educmae8a10 -5.40765 0.67303 -8.035 9.77e-16 ***
## qst19educmaeate 3 -8.31340 0.94971 -8.754 < 2e-16 ***
## qst19educmaenaosabe -8.22372 0.72973 -11.270 < 2e-16 ***
## qst38reprovSim 1vez -9.70824 0.54897 -17.685 < 2e-16 ***
## qst38reprovSim duas -9.87577 0.88970 -11.100 < 2e-16 ***
## abanAbandonou mais de uma vez -7.44178 2.83081 -2.629 0.00857 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 36.63 on 26086 degrees of freedom
## Multiple R-squared: 0.406, Adjusted R-squared: 0.4058
## F-statistic: 1783 on 10 and 26086 DF, p-value: < 2.2e-167 Verificação das hipoteses
residualPlot(mod)
hist(mod$residuals)
qqnorm(mod$residuals)
qqline(mod$residuals)
acf(mod$residuals)
vif(mod)## GVIF Df GVIF^(1/(2*Df))
## AcertosLP 1.117128 1 1.056943
## qst1sexo 1.034795 1 1.017249
## qst19educmae 1.062980 5 1.006126
## qst38reprov 1.082549 2 1.020027
## aban 1.008929 1 1.004455jarque.bera.test(mod$residuals)##
## Jarque Bera Test
##
## data: mod$residuals
## X-squared = 15.548, df = 2, p-value = 0.0004204Passa em todos os testes, menos na normalidade.
8 Análise de nível
plot(tapply(y2,pk_cod_entidade, mean)); abline(h=mean(y2))
9 Criando varíável artifical de nível 2
data$sx=ifelse(data$qst1sexo=="Feminin",0,1)
prop=aggregate(data$sx,list(data$pk_cod_entidade),mean)
names(prop) = c("pk_cod_entidade","prop")
dado=merge(data,prop,by="pk_cod_entidade")10 Modelo de intercepto aleatório
modlm0 <- lm(y2 ~ AcertosLP + qst1sexo+qst19educmae+qst38reprov+aban, data=dado)
modmm0 <- lmer(y2 ~ AcertosLP + qst1sexo+qst19educmae+qst38reprov+aban + ( 1 | pk_cod_entidade), data=dado, REML = F)
anova(modmm0,modlm0)## Data: dado
## Models:
## modlm0: y2 ~ AcertosLP + qst1sexo + qst19educmae + qst38reprov + aban
## modmm0: y2 ~ AcertosLP + qst1sexo + qst19educmae + qst38reprov + aban +
## modmm0: (1 | pk_cod_entidade)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## modlm0 12 262011 262109 -130993 261987
## modmm0 13 259745 259851 -129859 259719 2267.6 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1(sigma_2_u<- as.numeric(VarCorr(modmm0)) )## [1] 142.7751(sigma_2_e <- attr(VarCorr(modmm0), "sc")^2 )## [1] 1193.263(rho <- sigma_2_u / (sigma_2_u+sigma_2_e))## [1] 0.106864511 Modelo com interação
mod1 <- lmer(y2 ~ AcertosLP + qst1sexo+qst19educmae+qst38reprov+aban+( 1 | pk_cod_entidade), data=dado)
mod2 <- lmer(y2 ~ AcertosLP*qst1sexo + qst1sexo+qst19educmae+qst38reprov+aban + ( 1 | pk_cod_entidade), data=dado, REML = F)
anova(mod1 ,mod2)## refitting model(s) with ML (instead of REML)## Data: dado
## Models:
## mod1: y2 ~ AcertosLP + qst1sexo + qst19educmae + qst38reprov + aban +
## mod1: (1 | pk_cod_entidade)
## mod2: y2 ~ AcertosLP * qst1sexo + qst1sexo + qst19educmae + qst38reprov +
## mod2: aban + (1 | pk_cod_entidade)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mod1 13 259745 259851 -129859 259719
## mod2 14 259740 259854 -129856 259712 7.1972 1 0.007302 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 112 Modelo de coeficiente aleatório
mod1 <- lmer(y2 ~ AcertosLP*qst1sexo + qst1sexo+qst19educmae+qst38reprov+aban+( 1 | pk_cod_entidade), data=dado)
mod2 <- lmer(y2 ~ AcertosLP*qst1sexo + qst1sexo+qst19educmae+qst38reprov+aban + ( AcertosLP | pk_cod_entidade), data=dado, REML = F)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
## control$checkConv, : Model failed to converge with max|grad| = 0.0238219
## (tol = 0.002, component 1)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?anova(mod1 ,mod2)## refitting model(s) with ML (instead of REML)## Data: dado
## Models:
## mod1: y2 ~ AcertosLP * qst1sexo + qst1sexo + qst19educmae + qst38reprov +
## mod1: aban + (1 | pk_cod_entidade)
## mod2: y2 ~ AcertosLP * qst1sexo + qst1sexo + qst19educmae + qst38reprov +
## mod2: aban + (AcertosLP | pk_cod_entidade)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mod1 14 259740 259854 -129856 259712
## mod2 16 259603 259733 -129785 259571 141.06 2 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 113 Verificando se a proporção de homens é significante
mod1 <- lmer(y2 ~ AcertosLP*qst1sexo + qst1sexo+qst19educmae+qst38reprov+aban+( 1 | pk_cod_entidade), data=dado)
mod2 <- lmer(y2 ~ AcertosLP*qst1sexo + qst1sexo+qst19educmae+qst38reprov+aban + prop + ( 1 | pk_cod_entidade), data=dado, REML = F)
anova(mod1 ,mod2)## refitting model(s) with ML (instead of REML)## Data: dado
## Models:
## mod1: y2 ~ AcertosLP * qst1sexo + qst1sexo + qst19educmae + qst38reprov +
## mod1: aban + (1 | pk_cod_entidade)
## mod2: y2 ~ AcertosLP * qst1sexo + qst1sexo + qst19educmae + qst38reprov +
## mod2: aban + prop + (1 | pk_cod_entidade)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mod1 14 259740 259854 -129856 259712
## mod2 15 259735 259858 -129853 259705 6.3995 1 0.01141 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 114 Modelo final
Observação: Deu unini no modelo com coeficiente aleatório
final <- lmer(y2 ~ AcertosLP*qst1sexo + qst1sexo+qst19educmae+qst38reprov+aban + prop + ( 1 | pk_cod_entidade), data=dado, REML = F)
summary(final)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula:
## y2 ~ AcertosLP * qst1sexo + qst1sexo + qst19educmae + qst38reprov +
## aban + prop + (1 | pk_cod_entidade)
## Data: dado
##
## AIC BIC logLik deviance df.resid
## 259735.4 259857.9 -129852.7 259705.4 26082
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.9441 -0.6734 0.0035 0.6686 3.9891
##
## Random effects:
## Groups Name Variance Std.Dev.
## pk_cod_entidade (Intercept) 139.9 11.83
## Residual 1192.9 34.54
## Number of obs: 26097, groups: pk_cod_entidade, 356
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 173.54081 4.33537 40.029
## AcertosLP 4.76242 0.06750 70.554
## qst1sexoMasculi 12.04662 1.25106 9.629
## qst19educmae15m 1.09266 0.91600 1.193
## qst19educmae4ª7 -4.11864 0.63074 -6.530
## qst19educmae8a10 -3.70878 0.64121 -5.784
## qst19educmaeate 3 -6.04094 0.91505 -6.602
## qst19educmaenaosabe -7.01367 0.69587 -10.079
## qst38reprovSim 1vez -11.21265 0.52596 -21.318
## qst38reprovSim duas -12.21805 0.85384 -14.310
## abanAbandonou mais de uma vez -5.20878 2.70295 -1.927
## prop 22.62283 8.91046 2.539
## AcertosLP:qst1sexoMasculi 0.24724 0.09164 2.698##
## Correlation matrix not shown by default, as p = 13 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it15 Análise de resíduos
residuo=residuals(final)
qqnorm(residuo)
qqline(residuo)
hist(residuo)
jarque.bera.test(residuo)##
## Jarque Bera Test
##
## data: residuo
## X-squared = 1.3854, df = 2, p-value = 0.5002plot(residuo)
16 Equação
\(Y_{ij} = 184,1 + u_{j} + 4,7\times x_{1j} + 12,1\times SexoMasculino -7\times EducMaeNaoSabe\) \(- 6\times EducMae3 - 4,13\times EducMae4a7 - 3,7\times EducMae8a10+ 1,1\times EducMae15\)
\(-11,2\times reprovSim1vez -12,2\times reprovSim2vez -5,2\times Abandonou +\) \(+ 0.24 \times AcertosLP + 22,6\times Prop\)
17 Dicionário de variáveis
\(u_j\) = efeito específico de cada grupo
17.1 Nível 1
17.1.1 Contínua
\(x_{1j}\) = Acertos em língua portuguesa
17.1.2 Categórica
SexoMasculino = Se é do sexo masculino
EducMaeNaoSabe = Não sabe a escolaridade da mãe
EducMae3 = Não completou a quarta série
EducMae4a7 = Completou a quarta série, mas não a oitava.
EducMae8a10 = Completou a oitava série, mas não o ensino médio.
EducMae15 = Completou o ensino médio, mas não a faculdade
ReprovSim1vez = Reprovou uma vez
ReprovSim2vez = Reprovou duas vezes
Abandonou = Abandonou a escola pelo menos duas vezes
17.2 Nível 2
17.2.1 Contínua
Prop = proporção de homens para cada escola
17.2.2 Categórica
Nenhuma
18 Interpretação
O modelo final que obteve o melhor desempenho, utilizando como critério comparativo o método anova, foi o de interceptos aleatórios com interação. O modelo de coeficientes aleatório apresentou “anini”. Foi necessário artificializar uma variável de nível 2, pois nenhuma apresentada na base era capaz de explicar eficientemente.
Tem-se que, consirando as demais variaveis constantes, um acréscimo nos acertos de língua portuguesa, dado que o aluno é homem, verifica-se um aumento de 12,24 na média da nota de matemática se comparado com alunas do sexo feminino. Observa-se que a escolaridade materna tem um impacto negativo na nota de matemática, sendo que não saber tem o pior impacto e quanto maior a escolaridade, menor este impacto será. Em outras palavras, é esperado que quanto mais educação uma mãe tenha, maior seja a nota de matemática do filho, em média.
O mesmo efeito verifica-se para a quantidade de reprovações. Para alunos que reprovaram uma vez, tem-se que, em média, -11 de nota em matemática, já os que reprovaram duas, foram penalizados com -12 de nota.
E tem-se que se um aluno abandonou a escola mais de 1 vez, tem a nota de matemática penalizada em 5,2, na média, se comparado com quem abandanou uma vez ou nenhuma.
Por fim, quanto maior a proporção de homens em uma escola, maior será a nota média dos alunos desta, podendo atingir um máximo de 22 pontos a mais.