Modelo de Regressão Linear Multinível
Aplicado nas situações em que a amostra apresenta subgrupos (conglomerados), onde é suposto que as observações dentro de um conglomerado possam apresentar correlações (correlações intraclasse).
\(y_{ij}=\gamma_{00}+\gamma_{10}x_{ij}+u_{0j}+u_{1j}x_{ij}+e_{ij}\)
Exemplo SAEB
Proficiência : A pontuação final do aluno \(i\) na escola \(j\) . O valor base foi definido em 250 (escala Saeb).
\(NSE_{ij}\) (Nível 1): O nível socioeconômico afeta o aluno individualmente.
\(Tipo_{j}\) (Nível 2): Variável dicotômica (0=Pública, 1=Privada).
Abrindo dados
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr 1.1.4 ✔ readr 2.1.5
✔ forcats 1.0.0 ✔ stringr 1.5.1
✔ ggplot2 4.0.2 ✔ tibble 3.2.1
✔ lubridate 1.9.4 ✔ tidyr 1.3.1
✔ purrr 1.2.2
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
<- read_csv ("https://raw.githubusercontent.com/Silvioest/teste/refs/heads/master/dadosSAEB.csv" )
Rows: 1067 Columns: 4
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
dbl (4): id_escola, tipo_escola, nse, proficiencia
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
id_escola tipo_escola nse proficiencia
Min. : 1.00 Min. :0.0000 Min. :-3.24892 Min. :152.9
1st Qu.: 8.00 1st Qu.:0.0000 1st Qu.:-0.60700 1st Qu.:239.6
Median :16.00 Median :0.0000 Median : 0.04871 Median :264.0
Mean :15.48 Mean :0.3018 Mean : 0.03937 Mean :269.5
3rd Qu.:22.00 3rd Qu.:1.0000 3rd Qu.: 0.71188 3rd Qu.:299.7
Max. :30.00 Max. :1.0000 Max. : 3.47607 Max. :385.6
Correlação
nse proficiencia
nse 1.0000000 0.3786929
proficiencia 0.3786929 1.0000000
Histograma da variável Proficiência
hist (dados$ proficiencia,main = "Histograma da v.a Proficiência" ,xlab= "Proficiência" )
Gráfico de dispersão Proficência VS NSE
<- ggplot (dados, aes (x= nse, y= proficiencia))+ geom_point (shape= 20 )+ geom_smooth (method= lm, se= FALSE , size= .5 )+ theme (legend.position= "none" )
Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.
`geom_smooth()` using formula = 'y ~ x'
Gráfico de dispersão Proficência VS NSE por escola
<- ggplot (dados, aes (x= nse, y= proficiencia, color= as.factor (id_escola)))+ geom_point (shape= 20 )+ geom_smooth (method= lm, se= FALSE , size= .5 )+ theme (legend.position= "none" )
`geom_smooth()` using formula = 'y ~ x'
Gráfico de dispersão Proficência VS NSE por escola na forma de painel
<- p1 + facet_wrap (~ id_escola, nrow = 4 )
`geom_smooth()` using formula = 'y ~ x'
Ajustando o modelo nulo
Carregando pacotes exigidos: Matrix
Anexando pacote: 'Matrix'
Os seguintes objetos são mascarados por 'package:tidyr':
expand, pack, unpack
<- lmer (proficiencia ~ 1 + (1 | id_escola), data= dados)summary (ajuste0)
Linear mixed model fit by REML ['lmerMod']
Formula: proficiencia ~ 1 + (1 | id_escola)
Data: dados
REML criterion at convergence: 10028.2
Scaled residuals:
Min 1Q Median 3Q Max
-3.0530 -0.6610 0.0322 0.6620 3.0279
Random effects:
Groups Name Variance Std.Dev.
id_escola (Intercept) 1169.7 34.20
Residual 632.4 25.15
Number of obs: 1067, groups: id_escola, 30
Fixed effects:
Estimate Std. Error t value
(Intercept) 271.988 6.295 43.21
** Calculo do ICC - Efeito da escola **
\[ ICC=\frac{\sigma_{u0}^{2}}{\sigma_{e}^{2}+\sigma_{u0}^{2}}\]
<- 1169.7 / (1169.7+632.4 )
8.Efeito da variância da proficiencia do aluno dentro da escola
Coeficiente do modelo nulo
$id_escola
(Intercept)
1 308.5453
2 309.2242
3 272.9411
4 252.8599
5 309.5960
6 267.9398
7 262.2940
8 311.3494
9 333.1301
10 237.4679
11 316.9333
12 267.5575
13 324.5646
14 254.1270
15 263.0281
16 314.7944
17 252.4165
18 321.5714
19 227.6949
20 234.8530
21 238.3434
22 239.8882
23 262.1959
24 264.3514
25 247.9590
26 240.9145
27 242.0418
28 310.6395
29 242.5376
30 227.8711
attr(,"class")
[1] "coef.mer"
Modelo com NSE no primeiro nível
<- lmer (proficiencia ~ nse + (1 | id_escola), data= dados)summary (ajuste1)
Linear mixed model fit by REML ['lmerMod']
Formula: proficiencia ~ nse + (1 | id_escola)
Data: dados
REML criterion at convergence: 9598.9
Scaled residuals:
Min 1Q Median 3Q Max
-2.8896 -0.6510 0.0032 0.6607 3.5346
Random effects:
Groups Name Variance Std.Dev.
id_escola (Intercept) 1141.6 33.79
Residual 419.2 20.47
Number of obs: 1067, groups: id_escola, 30
Fixed effects:
Estimate Std. Error t value
(Intercept) 271.5954 6.2028 43.79
nse 15.0955 0.6561 23.01
Correlation of Fixed Effects:
(Intr)
nse -0.003
Coeficiente do modelo com NSE no primeiro nível
$id_escola
(Intercept) nse
1 304.7896 15.09552
2 312.0124 15.09552
3 266.7696 15.09552
4 253.3116 15.09552
5 314.1567 15.09552
6 267.5039 15.09552
7 261.7787 15.09552
8 313.0578 15.09552
9 327.3642 15.09552
10 235.8471 15.09552
11 316.2743 15.09552
12 270.1802 15.09552
13 325.6863 15.09552
14 249.0876 15.09552
15 261.9421 15.09552
16 311.8741 15.09552
17 249.9307 15.09552
18 319.3385 15.09552
19 228.0076 15.09552
20 236.5014 15.09552
21 241.3039 15.09552
22 243.0579 15.09552
23 263.4017 15.09552
24 264.2534 15.09552
25 250.6910 15.09552
26 240.5979 15.09552
27 241.8410 15.09552
28 309.7683 15.09552
29 242.0182 15.09552
30 225.5148 15.09552
attr(,"class")
[1] "coef.mer"
Modelo com NSE no primeiro nível e tipo de escola segundo nível
<- lmer (proficiencia ~ nse + factor (tipo_escola) + (1 | id_escola), data= dados)summary (ajuste2)
Linear mixed model fit by REML ['lmerMod']
Formula: proficiencia ~ nse + factor(tipo_escola) + (1 | id_escola)
Data: dados
REML criterion at convergence: 9532.5
Scaled residuals:
Min 1Q Median 3Q Max
-2.8545 -0.6629 -0.0064 0.6638 3.5052
Random effects:
Groups Name Variance Std.Dev.
id_escola (Intercept) 130.9 11.44
Residual 419.2 20.47
Number of obs: 1067, groups: id_escola, 30
Fixed effects:
Estimate Std. Error t value
(Intercept) 249.3898 2.6720 93.33
nse 15.1054 0.6553 23.05
factor(tipo_escola)1 66.5105 4.6501 14.30
Correlation of Fixed Effects:
(Intr) nse
nse -0.003
fctr(tp_s)1 -0.575 -0.006
Coeficiente do modelo com NSE no primeiro nível e tipo de escola segundo nível
$id_escola
(Intercept) nse factor(tipo_escola)1
1 239.3334 15.10542 66.51055
2 246.6309 15.10542 66.51055
3 265.5228 15.10542 66.51055
4 252.6035 15.10542 66.51055
5 248.2876 15.10542 66.51055
6 266.2941 15.10542 66.51055
7 260.7019 15.10542 66.51055
8 247.3278 15.10542 66.51055
9 260.2059 15.10542 66.51055
10 236.4486 15.10542 66.51055
11 250.0924 15.10542 66.51055
12 267.2898 15.10542 66.51055
13 258.7382 15.10542 66.51055
14 248.9420 15.10542 66.51055
15 260.9519 15.10542 66.51055
16 246.0177 15.10542 66.51055
17 249.7240 15.10542 66.51055
18 252.9618 15.10542 66.51055
19 229.0743 15.10542 66.51055
20 237.0570 15.10542 66.51055
21 241.6558 15.10542 66.51055
22 243.2632 15.10542 66.51055
23 262.1155 15.10542 66.51055
24 262.8641 15.10542 66.51055
25 250.3361 15.10542 66.51055
26 241.1187 15.10542 66.51055
27 242.2947 15.10542 66.51055
28 244.3022 15.10542 66.51055
29 242.3801 15.10542 66.51055
30 227.1581 15.10542 66.51055
attr(,"class")
[1] "coef.mer"
Modelo com NSE no primeiro nível e tipo de escola segundo nível e coef aleatório
<- lmer (proficiencia ~ nse + factor (tipo_escola) + (1 + nse | id_escola), data= dados)
boundary (singular) fit: see help('isSingular')
Linear mixed model fit by REML ['lmerMod']
Formula: proficiencia ~ nse + factor(tipo_escola) + (1 + nse | id_escola)
Data: dados
REML criterion at convergence: 9527.8
Scaled residuals:
Min 1Q Median 3Q Max
-2.8069 -0.6605 0.0000 0.6787 3.4903
Random effects:
Groups Name Variance Std.Dev. Corr
id_escola (Intercept) 131.951 11.487
nse 2.144 1.464 -1.00
Residual 417.062 20.422
Number of obs: 1067, groups: id_escola, 30
Fixed effects:
Estimate Std. Error t value
(Intercept) 249.4745 2.6219 95.15
nse 15.0229 0.7058 21.28
factor(tipo_escola)1 66.3658 4.3357 15.31
Correlation of Fixed Effects:
(Intr) nse
nse -0.320
fctr(tp_s)1 -0.547 0.014
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
Coeficiente do modelo com NSE no primeiro nível e tipo de escola segundo nível e coef aleatório
$id_escola
(Intercept) nse factor(tipo_escola)1
1 239.5172 16.29207 66.36578
2 246.7561 15.36936 66.36578
3 266.5831 12.84210 66.36578
4 252.6494 14.61817 66.36578
5 248.2120 15.18379 66.36578
6 266.4730 12.85612 66.36578
7 260.5280 13.61392 66.36578
8 247.8239 15.23326 66.36578
9 260.3690 13.63419 66.36578
10 236.3666 16.69367 66.36578
11 250.1187 14.94075 66.36578
12 267.2951 12.75135 66.36578
13 258.6936 13.84774 66.36578
14 248.8874 15.09770 66.36578
15 261.0683 13.54505 66.36578
16 246.0831 15.45515 66.36578
17 250.2143 14.92856 66.36578
18 252.8483 14.59281 66.36578
19 228.6651 17.67535 66.36578
20 237.2482 16.58129 66.36578
21 242.1275 15.95935 66.36578
22 243.0663 15.83969 66.36578
23 261.7322 13.46043 66.36578
24 262.9657 13.30319 66.36578
25 250.3187 14.91526 66.36578
26 240.3862 16.18131 66.36578
27 242.8595 15.86605 66.36578
28 245.1620 15.57256 66.36578
29 242.4129 15.92297 66.36578
30 226.8049 17.91246 66.36578
attr(,"class")
[1] "coef.mer"
16.Comparação entre os modelos
anova (ajuste0,ajuste1,ajuste2,ajuste3)
refitting model(s) with ML (instead of REML)
Data: dados
Models:
ajuste0: proficiencia ~ 1 + (1 | id_escola)
ajuste1: proficiencia ~ nse + (1 | id_escola)
ajuste2: proficiencia ~ nse + factor(tipo_escola) + (1 | id_escola)
ajuste3: proficiencia ~ nse + factor(tipo_escola) + (1 + nse | id_escola)
npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
ajuste0 3 10039.7 10054.6 -5016.9 10033.7
ajuste1 4 9613.4 9633.3 -4802.7 9605.4 428.3071 1 < 2.2e-16 ***
ajuste2 5 9551.7 9576.6 -4770.9 9541.7 63.6621 1 1.477e-15 ***
ajuste3 7 9550.9 9585.7 -4768.5 9536.9 4.8205 2 0.08979 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Gráfico de análise dos coeficientes com intervalo de confiança
library (lattice)dotplot (ranef (ajuste3, condVar = TRUE ))
Analise da qualidade do modelo
library (performance)check_model (ajuste3, check = "linearity" )check_model (ajuste3, check = "qq" )check_model (ajuste3, check = "reqq" )
AJUSTE OK
