Este capítulo amplía la regresión multinivel básica permitiendo que interceptos y pendientes varíen por grupo y abordando estructuras no anidadas (p. ej., individuos clasificados por demografía y por estado simultáneamente). Se muestra cómo interpretar modelos del tipo \(y_{i}=\alpha_{j[i]}+\beta_{j[i]}x_{i}+\cdots\) como interacciones entre el índice de grupo y predictores individuales, y cómo incluir predictores a nivel de grupo (p. ej., uranio del suelo) e interacciones con \(x\). Se presentan ajustes en R con lmer() y se discute la extracción/interpretación de efectos fijos y aleatorios.
3.1 Interceptos y pendientes que varían por grupo
Siguiente paso: permitir que más de un coeficiente varíe por grupo. Ejemplo (radón) con un predictor individual \(x\) (indicador de primer piso), sin uranio a nivel condado:
Rows: 3194 Columns: 7
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (4): stfips, ctfips, st, cty
dbl (3): lon, lat, Uppm
ℹ 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.
radon_data <- radon_data %>%mutate(fips = stfips *1000+ cntyfips)# Obtener el FIPS a nivel de EE.UU.usa.fips <-1000*as.numeric(cty$stfips) +as.numeric(cty$ctfips)# Encontrar las filas correspondientes a los condados únicos en Minnesotausa.rows <-match(unique(radon_data$fips), usa.fips)uranium <- cty[usa.rows, "Uppm"]u <-log(uranium)county_u_data <-data.frame(county_name =unique(radon_data$county_name),u = u)radon_data <- radon_data %>%left_join(county_u_data, by ="county_name") %>%rename(u = Uppm)# Supóngase u.full = u[county], donde u es el vector a nivel condadoM4 <-lmer(log_radon ~ floor + u + floor:u + (1+ floor | county_name), data = radon_data)
boundary (singular) fit: see help('isSingular')
summary(M4)
Linear mixed model fit by REML ['lmerMod']
Formula: log_radon ~ floor + u + floor:u + (1 + floor | county_name)
Data: radon_data
REML criterion at convergence: 2132.2
Scaled residuals:
Min 1Q Median 3Q Max
-5.2770 -0.6057 0.0507 0.6288 3.4495
Random effects:
Groups Name Variance Std.Dev. Corr
county_name (Intercept) 0.0000 0.0000
floor 0.1392 0.3731 NaN
Residual 0.5735 0.7573
Number of obs: 919, groups: county_name, 85
Fixed effects:
Estimate Std. Error t value
(Intercept) 1.44407 0.02926 49.353
floor -0.64131 0.08915 -7.194
u 0.84628 0.07476 11.320
floor:u -0.41582 0.23894 -1.740
Correlation of Fixed Effects:
(Intr) floor u
floor -0.328
u 0.354 -0.116
floor:u -0.111 0.155 -0.313
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
Interpretación de parámetros: los coeficientes fijos corresponden a \(\gamma_{0}^{\alpha}\) (Intercepto), \(\gamma_{0}^{\beta}\) (coef. de \(x\)), \(\gamma_{1}^{\alpha}\) (coef. de \(u\)), \(\gamma_{1}^{\beta}\) (coef. de \(floor\:u\)).
3.1.4 Coeficientes por grupo y efectos fijos/aleatorios
# Coeficientes por condado (filas) para (Intercepto), x, u.full, x:u.fullcoef(M4)
(Intercept) floor u floor:u
1.4440747 -0.6413080 0.8462757 -0.4158211
# Efectos aleatorios por condado en (Intercepto) y xranef(M4)
$county_name
(Intercept) floor
AITKIN 0 0.055139162
ANOKA 0 -0.314684718
BECKER 0 0.029175313
BELTRAMI 0 0.085028824
BENTON 0 0.063178258
BIG STONE 0 0.000000000
BLUE EARTH 0 0.421386890
BROWN 0 0.074853969
CARLTON 0 0.207639616
CARVER 0 -0.084171171
CASS 0 0.000000000
CHIPPEWA 0 0.000000000
CHISAGO 0 0.000000000
CLAY 0 0.023673273
CLEARWATER 0 -0.014828170
COOK 0 0.000000000
COTTONWOOD 0 -0.256492884
CROW WING 0 0.172484933
DAKOTA 0 -0.184233642
DODGE 0 0.000000000
DOUGLAS 0 0.085657656
FARIBAULT 0 -0.130359390
FILLMORE 0 -0.012711249
FREEBORN 0 0.116455394
GOODHUE 0 0.398568204
HENNEPIN 0 -0.167347834
HOUSTON 0 0.051971277
HUBBARD 0 0.031437565
ISANTI 0 0.000000000
ITASCA 0 0.000000000
JACKSON 0 0.000000000
KANABEC 0 0.000000000
KANDIYOHI 0 0.000000000
KITTSON 0 -0.032697578
KOOCHICHING 0 0.036012969
LAC QUI PARLE 0 0.290673290
LAKE 0 -0.047405686
LAKE OF THE WOODS 0 0.245876238
LE SUEUR 0 0.107158617
LINCOLN 0 0.050488676
LYON 0 0.187285350
MAHNOMEN 0 0.000000000
MARSHALL 0 -0.374677135
MARTIN 0 -0.306134829
MCLEOD 0 -0.240818966
MEEKER 0 0.000000000
MILLE LACS 0 -0.274527905
MORRISON 0 0.065624417
MOWER 0 -0.240617161
MURRAY 0 0.000000000
NICOLLET 0 0.000000000
NOBLES 0 0.000000000
NORMAN 0 -0.087296405
OLMSTED 0 -0.366194184
OTTER TAIL 0 0.172343023
PENNINGTON 0 -0.234888776
PINE 0 -0.107459085
PIPESTONE 0 0.087863261
POLK 0 -0.095081543
POPE 0 0.000000000
RAMSEY 0 0.273006641
REDWOOD 0 0.404260674
RENVILLE 0 0.154554776
RICE 0 -0.002038909
ROCK 0 0.000000000
ROSEAU 0 0.190310619
SCOTT 0 0.404415825
SHERBURNE 0 0.000000000
SIBLEY 0 0.000000000
ST LOUIS 0 -0.251697206
STEARNS 0 -0.155430480
STEELE 0 0.000000000
STEVENS 0 0.000000000
SWIFT 0 0.000000000
TODD 0 0.197113289
TRAVERSE 0 0.009327209
WABASHA 0 -0.150802666
WADENA 0 -0.036031487
WASECA 0 -0.195239440
WASHINGTON 0 -0.231147389
WATONWAN 0 0.375766801
WILKIN 0 0.000000000
WINONA 0 -0.487172809
WRIGHT 0 0.013456688
YELLOW MEDICINE 0 0.000000000
with conditional variances for "county_name"
Los coeficientes para el intercepto y \(floor\)varían por grupo (como se especificó); los de u y floor:u son fijos (porque u es constante dentro de condado).
3.1.5 De la salida de lmer a interceptos y pendientes por grupo
La recta ajustada en el condado \(j\) se expresa como
Cálculo a partir de coef(M4) (ilustrativo, asumiendo vector \(u\) a nivel condado):
# Supóngase u es un vector (longitud J) con el uranio por condado en el mismo ordena.hat.M4 <-coef(M4)$county[, 1] +coef(M4)$county[, 3] * ub.hat.M4 <-coef(M4)$county[, 2] +coef(M4)$county[, 4] * utibble(alpha_hat = a.hat.M4, beta_hat = b.hat.M4) %>%head()
donde \(v_i=u_{j[i]}\). Interpretaciones equivalentes:
Modelo con intercepto y pendiente variables y cuatro predictores individuales \((1, v_i, x_i, v_i x_i)\) con efectos centrados en cero.
Regresión con \(4+2J\) predictores (constante, \(v_i\), \(x_i\), \(v_i x_i\), indicadores de grupo y sus interacciones con \(x\)).
Regresión con cuatro predictores y tres términos de error.
Modelo multinivel original con un predictor de nivel de grupo.
Estas expresiones alternativas ayudan a conectar la formulación jerárquica con la implementación computacional e interpretación de interacciones y variación entre grupos.
3.2 Modelos no anidados
Este apartado aborda los modelos multinivel no anidados, en los cuales los individuos pueden pertenecer simultáneamente a múltiples estructuras de agrupamiento. Se presentan tres ejemplos representativos: (1) un experimento psicológico con pilotos en simuladores de vuelo (tratamientos y aeropuertos), (2) una regresión de ingresos que considera etnicidad, edad y estatura, y (3) un diseño experimental de cuadrado latino aplicado a cultivos.
3.2.1 Ejemplo 1: experimento psicológico con dos factores de agrupación
En el primer ejemplo (pilotos en simuladores de vuelo), los datos corresponden a un experimento psicológico en el que se evaluó el desempeño de pilotos entrenados en simuladores. Características principales:
Número de observaciones:\(n = 40\).
Factores de agrupamiento no anidados:
Tratamientos (group):\(J = 5\) condiciones experimentales distintas de entrenamiento.
Escenarios (scenario):\(K = 8\) aeropuertos distintos en los que se probó a los pilotos.
Variable de respuesta (\(y\)):
Tasa de éxito (proporción de aterrizajes exitosos o de recuperaciones correctas en el simulador).
Construida como:
\[
y = \frac{\text{número de éxitos}}{\text{número de intentos}}
\]
en cada combinación tratamiento × aeropuerto.
Variables en bruto en el archivo original (pilots.dat):
group: identificador de la condición de tratamiento asignada.
scenario: identificador del aeropuerto en el simulador.
recovered: indicador de éxito del intento (\(1 =\) éxito, \(0 =\) fallo, NA si falta el dato).
Estructura del dataset procesado:
Cada fila representa una combinación única de grupo × escenario.
Para cada celda se contabilizan los éxitos y fallos, y a partir de ello se calcula la proporción \(y\).
Esto da lugar a una matriz de tamaño \(8 \times 5\) (o sea, \(K \times J\)), donde cada celda resume el desempeño promedio de los pilotos en ese tratamiento y aeropuerto.
Modelo no anidado con \(J=5\) tratamientos y \(K=8\) aeropuertos (\(n=40\) observaciones):
display(M1) # muestra coeficientes fijos y desviaciones estándar por nivel
lmer(formula = y ~ 1 + (1 | group) + (1 | scenario), data = result)
coef.est coef.se
0.44 0.12
Error terms:
Groups Name Std.Dev.
scenario (Intercept) 0.32
group (Intercept) 0.00
Residual 0.22
---
number of obs: 40, groups: scenario, 8; group, 5
AIC = 20.3, DIC = 7.3
deviance = 9.8
Objetivo del bloque: Ajustar un modelo equivalente con nlme::lme() (estructura alternativa de aleatoriedad).
M1_nlme <- nlme::lme(y ~1, random =~1| group / scenario, data = result)summary(M1_nlme)
Linear mixed-effects model fit by REML
Data: result
AIC BIC logLik
45.53959 52.19384 -18.76979
Random effects:
Formula: ~1 | group
(Intercept)
StdDev: 7.146626e-06
Formula: ~1 | scenario %in% group
(Intercept) Residual
StdDev: 0.3734507 0.001986792
Fixed effects: y ~ 1
Value Std.Error DF t-value p-value
(Intercept) 0.4418155 0.05904858 35 7.482237 0
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-0.006293823 -0.004513152 -0.001248589 0.005916490 0.007951542
Number of Observations: 40
Number of Groups:
group scenario %in% group
5 40
Notas de interpretación: En un ajuste de este tipo, es común que la variación entre aeropuertos sea sustancial en comparación con la variación entre tratamientos, lo que se refleja en las desviaciones estándar estimadas a cada nivel.
3.2.2 Ejemplo 2: regresión de ingresos en etnicidad, edad y estatura
3.2.2.1 Descripción del dataset de ingresos y altura
Unidad de análisis: Individuos encuestados en EE.UU.
Tamaño aproximado: varios miles de observaciones (dependiendo de los filtros aplicados).
Variables principales en el dataset (heights.dta):
earn: ingresos anuales (variable continua, en dólares).
height: estatura en pulgadas.
sex: sexo del individuo (1 = mujer, 2 = hombre).
race: categoría racial (1 = blanca, 2 = negra, otros valores para otras categorías).
hisp: indicador de origen hispano (1 = sí, 0 = no).
ed: nivel educativo alcanzado (categórico/ordinal).
yearbn: año de nacimiento.
Variables derivadas en el preprocesamiento:
age: edad aproximada del individuo, calculada como $90 - $.
age_category: edad agrupada en tres categorías:
1 = joven (<35 años),
2 = mediana edad (35–49 años),
3 = mayor (50+ años).
eth: etnicidad en 4 categorías, combinando race y hisp:
Negra,
Hispana,
Blanca,
Otras.
male: indicador de sexo (1 = hombre, 0 = mujer).
y: logaritmo de los ingresos (log(earn)), usado como respuesta en la regresión.
x: estatura, a veces centrada en la media (x.centered) para mejorar la interpretación del modelo.
El dataset permite explorar cómo los ingresos dependen de:
Objetivo del bloque: Ajustar un modelo con intercepto y pendiente aleatorios por etnicidad.
M1 <-lmer(y ~ x + (1+ x | eth), data = heights_clean)
boundary (singular) fit: see help('isSingular')
display(M1)
lmer(formula = y ~ x + (1 + x | eth), data = heights_clean)
coef.est coef.se
(Intercept) 7.27 1.10
x 0.04 0.02
Error terms:
Groups Name Std.Dev. Corr
eth (Intercept) 1.55
x 0.02 -1.00
Residual 0.90
---
number of obs: 1062, groups: eth, 4
AIC = 2823.3, DIC = 2788.6
deviance = 2800.0
Objetivo del bloque: Ajustar el modelo análogo con nlme (estructura alternativa de covarianza).
M1_nlme <- nlme::lme(fixed = y ~ x, random =~ x | eth, data = heights_clean)summary(M1_nlme)
Linear mixed-effects model fit by REML
Data: heights_clean
AIC BIC logLik
2823.787 2853.583 -1405.893
Random effects:
Formula: ~x | eth
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
(Intercept) 1.77599145 (Intr)
x 0.02742268 -1
Residual 0.90300309
Fixed effects: y ~ x
Value Std.Error DF t-value p-value
(Intercept) 7.380834 1.2082771 1057 6.108560 0.0000
x 0.034270 0.0185317 1057 1.849284 0.0647
Correlation:
(Intr)
x -0.999
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-4.9284272 -0.4286750 0.1752600 0.6456474 2.6066249
Number of Observations: 1062
Number of Groups: 4
Objetivo del bloque:Centrar la estatura para reducir la correlación entre intercepto y pendiente aleatorios.
x.centered <- x -mean(x)M2 <-lmer(y ~ x.centered + (1+ x.centered | eth), data = heights_clean)
boundary (singular) fit: see help('isSingular')
display(M2)
lmer(formula = y ~ x.centered + (1 + x.centered | eth), data = heights_clean)
coef.est coef.se
(Intercept) 9.68 0.05
x.centered 0.03 0.02
Error terms:
Groups Name Std.Dev. Corr
eth (Intercept) 0.07
x.centered 0.03 1.00
Residual 0.90
---
number of obs: 1062, groups: eth, 4
AIC = 2823.3, DIC = 2788.6
deviance = 2800.0
Objetivo del bloque: Extender el modelo para incluir efectos por edad y por la interacción etnicidad×edad.
Notas de interpretación: Los errores residuales pueden ser grandes respecto a las variaciones a nivel de grupo, lo que limita la precisión predictiva a nivel individual. Las varianzas de interceptos y pendientes por bloques (etnicidad, edad e interacción) cuantifican la heterogeneidad sistemática entre grupos.
3.2.3 Ejemplo 3: diseño de cuadrado latino \(5 \times 5\)
3.2.3.1 Descripción del dataset
Unidad experimental: una parcela agrícola de mijo.
Diseño: un cuadrado latino \(5 \times 5\), es decir:
25 parcelas dispuestas en una cuadrícula de 5 filas × 5 columnas.
Se prueban 5 tratamientos distintos (por ejemplo, fertilizantes o métodos de cultivo).
Cada tratamiento aparece exactamente una vez en cada fila y en cada columna.
Variables principales:
row (fila)
Índice de la fila en la cuadrícula (1–5).
Representa una posible fuente de variación sistemática (e.g., gradiente de suelo en dirección norte–sur).
column (columna)
Índice de la columna en la cuadrícula (1–5).
Otra fuente potencial de variación sistemática (e.g., condiciones de riego en dirección este–oeste).
treatment (tratamiento)
Uno de los 5 tratamientos experimentales, asignado de manera que cada uno aparezca una vez por fila y por columna.
Estos tratamientos están ordenados (lo que permite modelar un posible efecto lineal).
yield (rendimiento)
Producción de mijo observada en cada parcela (variable continua).
\(\beta^{\text{row}}_j\): efecto de la fila \(j\).
\(\beta^{\text{column}}_k\): efecto de la columna \(k\).
\(\beta^{\text{treat}}_l\): efecto del tratamiento \(l\).
\(\gamma^{\text{row}}, \gamma^{\text{column}}, \gamma^{\text{treat}}\): tendencias lineales asociadas a cada dimensión.
\(\sigma_y\): variación residual (dentro de las parcelas).
El cuadrado latino permite controlar simultáneamente dos fuentes de variación sistemática (filas y columnas), para así evaluar el efecto de los tratamientos con mayor precisión.
En el modelo multinivel, además de estimar efectos independientes por fila, columna y tratamiento, se pueden incorporar predictores ordenados (lineales) para captar tendencias sistemáticas dentro de cada dimensión.
El centrado (restar 3) facilita interpretar \(\mu\) como media global y a las \(\gamma\) como tendencias lineales por fila, columna y tratamiento. Este ejemplo ilustra cómo, en modelos no anidados, es posible incorporar predictores a nivel de grupo (tendencias) junto con efectos aleatorios por cada factor, sin forzar una elección estricta entre un ajuste lineal y efectos específicos por nivel.
# A tibble: 1 × 3
n mean_yield sd_yield
<int> <dbl> <dbl>
1 25 252. 39.0
# 3) Visualización tipo “heatmap” (rendimiento por celda)ggplot(latin, aes(x = col, y = row, fill = yield)) +geom_tile(color ="white") +scale_fill_viridis_c() +labs(title ="Cuadrado latino 5x5: rendimiento por parcela",x ="Columna", y ="Fila", fill ="Rendimiento") +coord_equal() +theme_minimal()
# 6) Varianzas por nivel (resumen de incertidumbre a nivel fila/columna/tratamiento y residual)VarCorr(m_ls)
Groups Name Std.Dev.
row (Intercept) 26.831
col (Intercept) 0.000
treat (Intercept) 0.000
Residual 28.248
# 7) Efectos aleatorios estimados (desviaciones respecto a la tendencia)ranef_ls <-ranef(m_ls, condVar =TRUE)# 8) Graficar efectos por nivel (con intervalos aproximados)# a) Filas# 1) Extraer efectos aleatoriosranef_list <-ranef(m_ls)# 2) Convertir a tibble para filasranef_row <- ranef_list$row %>%rownames_to_column("grp") %>%as_tibble() %>%rename(est =`(Intercept)`)# 3) Graficar efectos aleatorios de Filaggplot(ranef_row, aes(x = grp, y = est)) +geom_hline(yintercept =0, linetype =2) +geom_point(size =3, color ="blue") +labs(title ="Efectos aleatorios por Fila (desvío del intercepto)",x ="Fila", y ="Efecto aleatorio") +theme_minimal()
# b) Columnas# Columnasranef_col <- ranef_list$col %>%rownames_to_column("grp") %>%as_tibble() %>%rename(est =`(Intercept)`)ggplot(ranef_col, aes(x = grp, y = est)) +geom_hline(yintercept =0, linetype =2) +geom_point(size =3, color ="darkgreen") +labs(title ="Efectos aleatorios por Columna (desvío del intercepto)",x ="Columna", y ="Efecto aleatorio") +theme_minimal()
# Tratamientosranef_trt <- ranef_list$treat %>%rownames_to_column("grp") %>%as_tibble() %>%rename(est =`(Intercept)`)ggplot(ranef_trt, aes(x = grp, y = est)) +geom_hline(yintercept =0, linetype =2) +geom_point(size =3, color ="red") +labs(title ="Efectos aleatorios por Tratamiento (desvío del intercepto)",x ="Tratamiento", y ="Efecto aleatorio") +theme_minimal()
# 9) Predicciones por celda (para diagnóstico rápido)latin %>%mutate(pred =predict(m_ls)) %>%ggplot(aes(x = yield, y = pred)) +geom_abline(slope =1, intercept =0, linetype =2) +geom_point() +labs(title ="Observado vs. Predicción",x ="Rendimiento observado", y ="Predicción") +theme_minimal()
3.3 Modelos clásicos para coeficientes de regresión
En esta sección se abordan estrategias para seleccionar, transformar y combinar predictores en modelos de regresión multinivel. Se discute cómo los modelos multinivel pueden extender los modelos clásicos, aprovechando información estructurada en diferentes niveles para mejorar las estimaciones y reducir la complejidad en problemas con gran número de predictores. También se exploran transformaciones lineales y combinaciones ponderadas de variables, conectándolas con el análisis factorial.
En un modelo multinivel, cada coeficiente tiene asociado un modelo con media y desviación estándar, lo cual generaliza la regresión clásica:
Si un predictor está “incluido”: equivale a \(\sigma_{\alpha}=\infty\), es decir, se estima directamente por mínimos cuadrados.
Si un predictor está “excluido”: equivale a \(\sigma_{\alpha}=0\), y el coeficiente se fija en cero.
De este modo, la regresión clásica puede verse como un caso particular de los modelos multinivel.
3.3.1 Modelos multinivel como alternativa a la selección de predictores
Los modelos multinivel permiten combinar múltiples entradas de manera más eficiente.
Ejemplo: Witte et al. (1994) ajustaron una regresión logística para predecir incidencia de cáncer en un estudio de 362 personas con consumo reportado de 87 alimentos.
