Cuando tenemos datos longitudinales o medidas repetidas podemos usar un modelo multinivel donde la medida repetida está anidada dentro de los individuos.1 Aunque en psicología generalmente la unidad de análisis es el individuo también se puede tomar como unidad, por ejemplo, un país respecto del cual se tienen múltiples observaciones en el tiempo Según distintos autores (Bryk & Raudenbush, 2002; Hox, 2010) hay múltiples ventajas asociadas al uso de modelos multinivel para diseños de medidas repetidas:
Para cada individuo existe un modelo de crecimiento lo que es consistente con la manera en que se conceptualiza el desarrollo individual
El número de medidas y el espaciamiento entre dichas medidas puede ser distinto para cada individuo
La covarianza entre las medidas repetidas se puede especificar
Es posible estimar un modelo incluso habiendo datos perdidos2 Un ANOVA de medidas repetidas (RM-ANOVA) sólo correría para un mismo número de medidas por individuo. RM-ANOVA es un caso especial de un modelo multinivel usando restricted maximum likelihood
Es muy sencillo añadir más niveles al modelo (e.g., individuos anidados en familias, escuelas, etc.)
Es fácil añadir variables que varían con el tiempo (time-varying covariates) o variables explicativas constantes (time-constant explanatory variables)
La flexibilidad en el modo de codificar la variable tiempo permite evaluar modelos que no están disponibles en RM-ANOVA, particularmente cuando se plantean modelos discontinuos3 Este detalle lo propongo yo como una ventaja adicional
La data
Para este ejemplo usaremos la data asociada al libro Multilevel Analysis: Techniques and Applications (2010) de Hox. Un archivo zip que contiene la data de Curran (1997) está disponible aquí. Específicamente usaremos el archivo CurranLong.sav. Como este es un archivo SPSS lo importaremos usando el paquete haven.
Un análisis preliminar de los datos4 Esta análisis lo hice usando la base CurranData.sav que está en formato wide indica que hay un número importante de datos perdidos. Esto es usual en diseños longitudinales donde la muestra inicial se va reduciendo con cada nueva ola de medición, o donde algunos casos tiene mediciones en un momento pero no en otro.
donde \(T{ti}\) es un indicador para la variable tiempo y \(X_{ti}\) es otra variable predictora que va cambiando en el tiempo. Veremos más adelante un modelo que usa una variable tiempo fija (e.g., tiempo 0, 1, 2, 3) y otro que usa una variable de tiempo variable. La variación entre individuos se especifica añadiendo variables de Nivel 2, que en este caso son variables individuales y efectos cruzados de Nivel 2 con el Nivel 1. Esto es, cada individuo tiene su intercepto pero además se puede incluir una pendiente, una variable predictora de nivel individual y las interacciones entre todos los componentes del modelo.
A continuación intentaremos reproducir los resultados de la Tabla 14.2 de Hox. El modelo inicial sólo tiene un intercepto y una variable de agrupación que es el id del individuo. Para todos los modelos usaremos maximum likelihood ya que se deduce, dado los resultados obtenidos, que ese es el método de estimación utilizado para construir la Tabla 14.2.
Es importante consignar que este modelo NO es el más indicado para evaluar el efecto de los predictores sustantivos del análisis. Más adelante usaremos como benchmark el modelo que incorpora la variable tiempo.
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: read ~ 1 + (1 | id)
## Data: d1_long
##
## AIC BIC logLik deviance df.resid
## 5057.8 5073.4 -2525.9 5051.8 1322
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.21973 -0.79677 0.02322 0.72395 2.34508
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.2978 0.5458
## Residual 2.3903 1.5461
## Number of obs: 1325, groups: id, 405
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 4.11298 0.05087 80.86
Means-as-outcomes adding occasion as predictor
En el segundo modelo incluimos la variable occasion como la variable predictora, que en este caso es el tiempo en que se tomó la medición y que va de 0 (primera medición) a 3 (cuarta medición). Es crítico consignar que el valor del intercepto corresponde a la media de read en la primera medición (occasion = 0). Esta codificación es recomendable para interpretar el valor del intercepto con facilidad.
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
## control$checkConv, : Model failed to converge with max|grad| = 0.00211389
## (tol = 0.002, component 1)
summary(model.2)
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: read ~ occasion + (occasion | id)
## Data: d1_long
##
## AIC BIC logLik deviance df.resid
## 3383.8 3414.9 -1685.9 3371.8 1319
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.7183 -0.5209 -0.0229 0.4803 4.1920
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## id (Intercept) 0.57102 0.7557
## occasion 0.07406 0.2721 0.29
## Residual 0.34590 0.5881
## Number of obs: 1325, groups: id, 405
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.69613 0.04524 59.59
## occasion 1.11907 0.02165 51.69
##
## Correlation of Fixed Effects:
## (Intr)
## occasion -0.155
## convergence code: 0
## Model failed to converge with max|grad| = 0.00211389 (tol = 0.002, component 1)
Cross-level effects (time constant variable)
En el modelo 4 de la Tabla 14.2 Hox incluye la edad del niño kidage, la presencia de estimulación cognitiva homecog y su interacción respectiva con la variable tiempo occasion. Existe una pequeña discrepancia entre los resultados obtenidos aquí y los resultados reportados por Hox en relación al efecto de occasion que trasunta en una diferencia en Deviance.
En los modelo 4 de la Tabla 14.3 Hox compara los resultados de occasion con kidage6. La variable kidage6 resulta de dos transformaciones que no están explicitadas en el capítulo pero que Hox describe en su libro. La trasformación consiste en sumarle 2, 4, y 8 años a la edad inicial del niño para que refleje la diferencia de tiempo entre administraciones de la prueba (hubiese sido mejor tener la edad exacta). En la data, el resultado de esta transformación aparece como kidagetv. La segunda transformación consiste en restarle 6 a la edad resultante para que la edad mínima sea 0; así obtenemos kidage6.
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
## control$checkConv, : Model failed to converge with max|grad| = 0.0193925
## (tol = 0.002, component 1)
summary(model.4a)
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: read ~ kidage6 + (kidage6 | id)
## Data: d1_long
##
## AIC BIC logLik deviance df.resid
## 3238.8 3269.9 -1613.4 3226.8 1319
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.5910 -0.5586 -0.0356 0.5002 5.2980
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## id (Intercept) 0.16659 0.4081
## kidage6 0.01761 0.1327 0.87
## Residual 0.35948 0.5996
## Number of obs: 1325, groups: id, 405
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.16429 0.03724 58.12
## kidage6 0.56432 0.01067 52.88
##
## Correlation of Fixed Effects:
## (Intr)
## kidage6 -0.240
## convergence code: 0
## Model failed to converge with max|grad| = 0.0193925 (tol = 0.002, component 1)
Por último, podemos comparar el ajuste de todos los modelos en un sólo paso mediante la función anova. Como señala Hox, sólo podemos usar la desviación para comparar modelos cuando uno de los modelos se obtiene añadiendo (o extrayendo) predictores del modelo de comparación. Cuando este no es el caso (como aquí) podemos usar el Akaike Information Criterion (AIC) que es la desviación más dos veces el número de parámetros estimados de un modelo. Por ejemplo, el modelo 2b tiene 8 parámetros (1 intercepto, 1 efecto de occasion, 1 residuo, 1 efecto de individuo, y 4 medidas repetidas) por lo que AIC = 3477.1 + 8 = 3485.1. Valore menores de AIC reflejan un mejor ajuste.
En comparación con el model.1 el modelo model.2 incluye la edad de la madre, la estimulación cognitiva y el soporte emocional. Para obtener los mismos resultados de Hox es imprescindible centrar estas variables.
#Estoy usando el promedio de la data wide porque si no las medias estarían influidas por la presencia de medidas repetidas perdidas en la data long
mean(d1_wide$momage)
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: read ~ kidage6 + kidage6sq + momagec + homecogc + homeemoc +
## (1 | id)
## Data: d1_long
##
## AIC BIC logLik deviance df.resid
## 3232.2 3273.7 -1608.1 3216.2 1317
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.8173 -0.5510 -0.0589 0.4919 3.8358
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.6016 0.7756
## Residual 0.3858 0.6212
## Number of obs: 1325, groups: id, 405
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.744493 0.060144 29.005
## kidage6 0.924936 0.028773 32.146
## kidage6sq -0.049201 0.003642 -13.510
## momagec 0.053691 0.023635 2.272
## homecogc 0.050708 0.018031 2.812
## homeemoc 0.039233 0.019960 1.966
##
## Correlation of Fixed Effects:
## (Intr) kidag6 kdg6sq momagc homcgc
## kidage6 -0.652
## kidage6sq 0.549 -0.961
## momagec 0.013 -0.017 0.015
## homecogc 0.001 -0.002 0.000 -0.176
## homeemoc -0.001 -0.001 0.003 -0.150 -0.324
En el modelo 3 se permite que la relación entre edad y read varíe por sujeto.5 Nótese que no incluyo en este modelo el término cuadrático como efecto aleatorio ya que su variabilidad es muy pequeña 0.0009338
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: read ~ kidage6 + kidage6sq + momagec + homecogc + homeemoc +
## (1 + kidage6 | id)
## Data: d1_long
##
## AIC BIC logLik deviance df.resid
## 3035.4 3087.3 -1507.7 3015.4 1315
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1340 -0.5312 -0.0466 0.4488 5.3381
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## id (Intercept) 0.21245 0.4609
## kidage6 0.01859 0.1363 0.57
## Residual 0.27807 0.5273
## Number of obs: 1325, groups: id, 405
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.748517 0.045791 38.185
## kidage6 0.920582 0.026082 35.295
## kidage6sq -0.048115 0.003252 -14.796
## momagec 0.026514 0.019333 1.371
## homecogc 0.037024 0.014866 2.490
## homeemoc 0.003082 0.016599 0.186
##
## Correlation of Fixed Effects:
## (Intr) kidag6 kdg6sq momagc homcgc
## kidage6 -0.638
## kidage6sq 0.619 -0.921
## momagec 0.042 -0.046 0.043
## homecogc 0.009 -0.015 0.017 -0.175
## homeemoc -0.022 0.025 -0.025 -0.145 -0.333
Por último incluimos la interacción entre la edad de la madre y el soporte emocional recibido en casa con la edad del niño para predecir read6 Al parecer hay un error en el libro de Hox ya que los valores fijos no coinciden con los resultados obtenidos en este paso
