Modelos multinivel para más de dos niveles

Correspondiente a la sesión del jueves, 5 de junio de 2025

Cargar paquetes

pacman::p_load(car, # función recode
                stargazer, # tablas latex
                lme4, # multilevel
                lattice, # dotplot
                sjPlot, # multilevel plots
                dplyr,
                doBy, # agregar datos
                VGAM, # predict logit
                texreg,
                nlme,
                dplyr,
                sjmisc,
                sjPlot)

options(scipen = 999) # para desactivar notacion cientifica

Modelos de 3 niveles

Ejemplo Finch Cap. 4, pag. 80

Lectura y limpieza datos (Achieve, base de datos de estudiantes en clases en escuelas)

“geread”: puntaje en lectura
“gevocab”: puntaje en vocabulario
“gender”: genero
“clenroll”: tamaño de los cursos
“cenroll”: tamaño del colegio
“class”: id cursos
“school”: id escuelas
“corp”: id distritos

Achieve <- read.csv("https://github.com/cursos-metodos-facso/multinivel-facso/raw/refs/heads/main/assignment/data/Achieve.csv")
names(Achieve)

myvars <- c("geread", # puntaje en lectura
            "gevocab", #puntaje en vocabulario
            "gender", # genero
            "clenroll", # tamaño de los cursos
            "cenroll", # tamaño del colegio
            "class", # id cursos
            "school", # id escuelas
            "corp" # id distritos
            ) # seleccionar variables
Achieve <- Achieve[myvars] # definir data con variables seleccionadas
names(Achieve)
attach(Achieve)

# para descriptivos por nivel
length(unique(school)) # Número escuelas
length(unique(class)) # Número (máximo) salas clases por escuela
Achieve$class2=(school*10000)+class # id único por clase
length(unique(Achieve$class2)) 

AchieveL2 = Achieve %>% group_by(class2) %>% 
            summarise_all(funs(mean))

Warning: `funs()` was deprecated in dplyr 0.8.0.
ℹ Please use a list of either functions or lambdas:

# Simple named list: list(mean = mean, median = median)

# Auto named with `tibble::lst()`: tibble::lst(mean, median)

# Using lambdas list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))

AchieveL2 = as.data.frame(AchieveL2) # para stargazer
AchieveL3 = Achieve %>% group_by(school) %>% 
            summarise_all(funs(mean))

Warning: `funs()` was deprecated in dplyr 0.8.0.
ℹ Please use a list of either functions or lambdas:

# Simple named list: list(mean = mean, median = median)

# Auto named with `tibble::lst()`: tibble::lst(mean, median)

# Using lambdas list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))

AchieveL3 = as.data.frame(AchieveL3)

Algunos descriptivos

stargazer(Achieve[c("geread", "gevocab","gender")],title="Estadísticos descriptivos L1", type = "html")

**Estadísticos descriptivos L1**

Statistic	N	Mean	St. Dev.	Min	Max

geread	10,320	4.341	2.332	0.000	12.000
gevocab	10,320	4.494	2.368	0.000	11.200
gender	10,320	1.502	0.500	1	2

stargazer(AchieveL2[c("clenroll","class")],title="Estadísticos descriptivos L2", type = "html")

**Estadísticos descriptivos L2**

Statistic	N	Mean	St. Dev.	Min	Max

clenroll	568	20.944	3.080	12	31
class	568	2.555	1.418	1	8

stargazer(AchieveL3[c("cenroll","school")],title="Estadísticos descriptivos L3", type = "html")

**Estadísticos descriptivos L3**

Statistic	N	Mean	St. Dev.	Min	Max

cenroll	160	7,494.138	8,656.893	110	41,359
school	160	4,522.288	2,533.673	767	9,173

table(Achieve$school, Achieve$class) # cursos por escuela y N

Modelo Nulo

Model4.7 <- lmer(geread~1+(1|school/class), data = Achieve)
sjPlot::tab_model(Model4.7, show.ci = FALSE)

	geread
Predictors	Estimates	p
(Intercept)	4.31	<0.001
Random Effects
σ²	4.85
τ₀₀ _class:school	0.27
τ₀₀ _school	0.31
ICC	0.11
N _class	8
N _school	160
Observations	10320
Marginal R² / Conditional R²	0.000 / 0.108

Correlación intra-clase

0.273/(0.273+4.847) # Correlación intra clase de class
0.312/(0.312 + 0.273 + 4.847) # Correlación intra clase de schools

Modelos con predictores

Model4.8 <- lmer( geread~gevocab+clenroll+cenroll+(1|school/class), data = Achieve)

Warning: Some predictor variables are on very different scales: consider
rescaling

Cuatro niveles

Model4.10 <- lmer( geread~1+(1|corp/school/class), data = Achieve)

Pendiente aleatoria

Model4.11  <-  lmer( geread~gevocab+gender+(gender|school/class), data = Achieve)

Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Model failed to converge with max|grad| = 0.0625298 (tol = 0.002, component 1)

sjPlot::tab_model(Model4.8, Model4.10, Model4.11,
                  show.ci = FALSE,
                  dv.labels = c("Modelo con predictores", "Cuatro niveles", "Pendiente aleatoria"))

	Modelo con predictores		Cuatro niveles		Pendiente aleatoria
Predictors	Estimates	p	Estimates	p	Estimates	p
(Intercept)	1.68	<0.001	4.33	<0.001	2.02	<0.001
gevocab	0.51	<0.001			0.51	<0.001
clenroll	0.02	0.047
cenroll	-0.00	0.307
gender					0.02	0.658
Random Effects
σ²	3.70		4.85		3.69
τ₀₀	0.09 _class:school		0.28 _{class:(school:corp)}		0.13 _class:school
	0.08 _school		0.09 _school:corp		0.04 _school
			0.18 _corp
τ₁₁					0.01 _{class:school.gender}
					0.01 _{school.gender}
ρ₀₁					-0.56 _class:school
					0.20 _school
ICC	0.04		0.10		0.04
N	8 _class		8 _class		8 _class
	160 _school		160 _school		160 _school
			59 _corp
Observations	10320		10320		10320
Marginal R² / Conditional R²	0.275 / 0.306		0.000 / 0.100		0.273 / 0.306

Modelos de 3 niveles

Foro