mixed.model.fnc

OBJETIVO

Lleva cabo la estimación de modelos de efectos mixtos (fijos y aleatorios) Mix model y Mix Logit model (para variables binomiales como aciertos, etc).

El modelo Mix-Model asume la estimación de 2 tipos de efectos: los fijos y los aleatorios. En el contexto del paradigma de presentar una muestra de items a una muestra de sujetos significa asumir que junto a los efectos esperados por la manipulación experimental existe una fuente de varianza situada en la variación aleatoria de los sujetos y otra perteneciente a los items. Además el modelo asume la existencia cruzada de ambas fuentes de varianza (los sujetos que fluctúan aleatoriamente respecto a una media común y los items que también lo hacen pero obviamente en general los sujetos pasan por todos los items). Los mix-model vienen a mejorar sustancialmente la pérdida de potencia que se produce en el análisis de la minF para conseguir la generalización de los efectos fijos encontrados de forma simultánea a los sujetos y los items).

mixed.model.fnc

Como ejemplo utilizaremos el ya conocido utilizado en la estimación de un modelo Anova F1, F2 y minF. Acude a esa página si necesitas comprender la naturaleza de los datos y el diseño.

datos.limpios=lee.archivo.fnc('datos_limpios.Rdata')

#--------------------------------------------------------------------------------

# Lectura de archivo binario Rdata

#--------------------------------------------------------------------------------

*** Se ha leido correctamente el archivo datos_limpios.Rdata

*** Esta es la cabecera:

El objeto es una matriz que pertenece a la clase data.frame

y tiene 5160 filas y 8 columnas (variables)

vd item sujeto mrA mrB condicion vd.original vd.eliminado

1 631.94 item1 suj1 low short low#short 631.94 0

2 573.98 item1 suj2 low short low#short 573.98 0

3 553.66 item1 suj3 low short low#short 553.66 0

4 532.54 item1 suj4 low short low#short 532.54 0

5 624.92 item1 suj5 low short low#short 624.92 0

6 716.63 item1 suj6 low short low#short 716.63 0

En el siguiente ejemplo solicitamos 2 modelos: 1 con solo el intercepto aleatorio para sujetos e items y un segundo al que añadimos pendiente aleatoria en sujetos para todos los efectos (A,B y AxB). Dado que dichos factores son intergrupo por items solo incluiremos intercepto aleatorio en item.

Mediante el argumento modelo el usuario debe incluir de forma explícita el modelo que desea estimar.

Estimaremos además los contrastes de efectos simples de la interacción de ambos factores incluyendo el argumento poshoc.

modelo1= 'vd ~ mrA*mrB + (1|sujeto) + (1|item)'

modelo2= 'vd ~ mrA*mrB + (mrA+mrB|sujeto) + (1|item)'

res1=mixed.model.fnc(datos.limpios, modelo=modelo1, poshoc='mrA:mrB')

res2=mixed.model.fnc(datos.limpios, modelo=modelo2, poshoc='mrA:mrB', grafica=T)

Si observas la última figura de la de la primera gráfica (Residuales), puedes ver que los residuales no tienen la distribución normal esperada por el modelo. En este caso deberíamos reestimar el modelo en dos fases: en la primera incluiremos el argumento boxcox=T para estimar el parámetro de potencia (o logarítmico) adecuado para transformar la variable criterio y mejorar o solucionar el problema de la incorrecta distribución de los residuos. Posteriormente debemos crear la nueva variable a partir de la potencia estimada y re-difiniremos el modelo sustituyendo la variable dependiente por la nueva creada.

res2=mixed.model.fnc(datos.limpios, modelo=modelo2, boxcox=T)

Extraemos del almacen res2 solo el apartado boxcox de la salida

res2$boxcox

bcPower Transformation to Normality

Est.Power Std.Err. Wald Lower Bound Wald Upper Bound

vd -0.9482 0.0562 -1.0585 -0.838

Likelihood ratio tests about transformation parameters

LRT df pval

LR test, lambda = (0) 292.6936 1 0

LR test, lambda = (1) 1271.2201 1 0

Vemos que rechazamos la hipótesis nula de la necesidad de una transformación logarítmica (lambda 0) como la relativa a que el exponente adecuado para los datos sea la constante 1. Por ello debemos aplicar la transformación de la vd mediante el procedimiento de boxcox utilizando para ello la constante de potencia estimada -0.9482:

datos.limpios$vd.p=(datos.limpios$vd^-0.9482 -1)/-0.9482

En relación a los dos modelos anteriormente estimados, podemos compararlos y valorar si es necesaria o no la pendiente aleatoria para ambos factores por sujeto.

compara.modelos.fnc(mod1=res1$estimado, res2$estimado)

#------------------------------------------------------------------

# COMPARACION DE MODELOS MIXED

#------------------------------------------------------------------

refitting model(s) with ML (instead of REML)

Data: dat.st

Models:

object: vd ~ mrA * mrB + (1 | sujeto) + (1 | item)

..1: vd ~ mrA * mrB + (mrA * mrB | sujeto) + (1 | item)

Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)

object 7 58211 58256 -29098 58197

..1 16 58206 58309 -29087 58174 23.425 9 0.005309 **

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Podemos ver que el segundo modelo con intercepto y pendiente aleatoria para ambos factores, genera la menor deviance, AIC y BIC por lo que a pesar de consumir 9 parámetros más y se por tanto menos parsimoniosos debemos seleccionar ese modelo frente al de intercepto random como único efecto aleatorio a considerar.

En el caso de que desees llevar a cabo la estimación del modelo sobre los errores, deberás incluir el argumento familia='logistica'. Por supuesto deberás modificar el modelo con el nombre de la variable con los errores y aciertos (0,1). Supongamos que tu variable se llama aciertos el modelo sería:

modelo= 'aciertos ~ mrA*mrB + (mrA|sujeto) + (mrA|item)'

mixed.model.fnc(datos.limpios, modelo=modelo, poshoc='mrA:mrB',

familia='logistica')

Evidentemente este último modelo de ejemplo propuesto, no puede ser estimado dado que no disponemos en la base de datos datos.limpios de esa variable de aciertos/error

Modelo de regresión multinivel

La función mixed.model.fnc será también la que utilizaremos para estimar modelos jerárquicos clásico y también los de curvas de crecimiento.

En el ambos tipos se asume la existencia de una jerarquía en los datos. El ejemplo más común en la literatura asume dos niveles jerárquicos: primer nivel jerárquico alumnos y segundo nivel escuelas. En su modelo mas simple se asume que la media global de la variable criterio (Intercepto) es diferente en cada colegio (segundo nivel jerárquico). Al descomponer la varianza total de la variable criterio en los dos niveles jerárquicos el usuario puede estudiar y tratar de explicar ambas fuentes de varianza incluyendo variables pertenecientes a ambos niveles (variables de los alumnos y de los colegios) así como su efecto de interacción cruzada.

Para este ejemplo de modelo jerárquico o multinivel, partiremos de las conocidas bases de datos de alumnos y colegios de Brik y Raudenbush(1992) (Extraido de la página de J.Fox apéndice de R compainon to applied regression) La base de datos alumnos contiene los valores de 7165 niños en 160 colegios diferentes. La base de datos colegios contiene por otra parte 6 variables características de cada uno de los 160 colegios. Debemos por tanto fundir ambas bases de datos, lo cual haremos con la función fundir.objetos.fnc.

colegios= lee.archivo.fnc('jerarquico_colegios.Rdata')

alumnos= lee.archivo.fnc('jerarquico.alumnos.Rdata')

head(alumnos)

School Minority Sex SES MathAch MEANSES

1 1224 No Female -1.528 5.876 -0.428

2 1224 No Female -0.588 19.708 -0.428

3 1224 No Male -0.528 20.349 -0.428

4 1224 No Male -0.668 8.781 -0.428

5 1224 No Male -0.158 17.898 -0.428

6 1224 No Male 0.022 4.583 -0.428

head(colegios)

School Size Sector PRACAD DISCLIM HIMINTY

1224 1224 842 Public 0.35 1.597 0

1288 1288 1855 Public 0.27 0.174 0

1296 1296 1719 Public 0.32 -0.137 1

1308 1308 716 Catholic 0.96 -0.622 0

1317 1317 455 Catholic 0.95 -1.694 1

1358 1358 1430 Public 0.25 1.535 0

Fundimos ambas bases de datos, con School como variable de emparejamiento de ambas bases de datos (que.var='School')

alumnos=fundir.objetos.fnc(alumnos,colegios, que.var='School')

#------------------------------------------------------------------

# FUNDIDO DE ARCHIVOS. ADD VARIABLES

#------------------------------------------------------------------

*** WARNING. Tu unidad de registro: School en la la primera base de

datos ocupa mas de una fila.

*** Dimensiones de entrada y salida ***

dim.datos1 dim.datos2 dim.combinados

filas 7185 160 7185

columnas 6 6 11

Una cuestión de especial relevancia tiene que ver con la naturaleza de la variable que define el segundo nivel jerárquico (School en este ejemplo). Es importante sobre todo a la hora de solicitar gráficos exploratorios que dicha variable sea categórica o factor. Esto lo conseguiremos fácilmente transformando dicha variable mediante:

alumnos$School= factor(alumnos$School)

Queremos modelar el rendimiento en matemáticas de los alumnos (MathAch) a partir de variables características de los alumnos y los colegios. Tenemos dos tipos de colegios recogidos en la variable Sector (públicos y católicos) y deseamos en primer lugar valorar la importancia que el nivel socioeconómico tiene sobre el rendimiento en matemáticas para ambos tipos de colegios.

Valoraremos en primer lugar un supuesto de extraordinaria importancia en los modelos lineales y en especial en los jerárquicos. Hablamos de la necesaria relación lineal que debe existir entre la variable criterio a modelar y las predictoras. Para ello solicitaremos un gráfico de panel que nos permita ver dicha relación para ambos tipos de colegios (incluyendo la variable Sector en el argumento x.panel).

grafica.panel.fnc(alumnos,vd='MathAch', que.factor='SES', x.panel='Sector',

regresion=T)

Vemos que la pendiente de regresión es mas acusada en los colegios de carácter público que en los católicos. Es el momento de plantearnos si esta pendiente pudiera tener un caracter aleatorio (random) a través del las escuelas (segundo nivel jerárquico), o lo que es lo mismo si el efecto del nivel socioeconómico sobre el rendimiento en matemáticas es diferente en cada escuela. Dado que tenemos una gran cantidad de colegios, no tiene sentido solicitar una gráfica de panel genere tantas subgráficas como colegios haya. Resulta altamente conveniente extraer una muestra aleatoria de 20 colegios de ambos tipos de la variable Sector. Utilizaremos para este objetivo la función extrae.muestra.fnc.

muestra = extrae.muestra.fnc(alumnos, que.factor='Sector', n=20, ID='School')

#------------------------------------------------------------------

# EXTRACCION DE MUESTRA ALEATORIA

#------------------------------------------------------------------

*** Se ha extraido una muestra de 20 registros por cada nivel de Sector ***

Por conveniencia, seguidamente crearemos una lista donde separaremos los colegios públicos y católicos utilizando divide.por.factor.fnc.

x.sector=divide.por.factor.fnc(muestra, que.factor='Sector')

#------------------------------------------------------------------

# DIVIDE LA BASE DE DATOS POR FACTOR

#------------------------------------------------------------------

*** Se ha creado una lista con 2 elementos, correspondientes a los niveles del factor Sector

*** Estos son los elementos de esa lista: Public Catholic

*** Si deseas acceder a uno de ellos indicalo mediante: nombre.de.la.lista$que.lista

*** donde que.lista se refiere al elemento que deseas extraer.

Ahora podemos solicitar dos gráficos de panel (uno para cada tipo de colegio) donde podamos valorar la pendiente de regresión del rendimiento en matemáticas de los alumnos sobre el nivel socioeconómico de los mismos. A diferencia de la gráfica anterior el argumento x.panel ahora tiene como valor la variable School. Es decir habrá 20 subgráficas o paneles (1 por colegio).

grafica.panel.fnc(x.sector$Public, vd='MathAch',

que.factor='SES',

x.panel='School',

regresion=T,

titulo='Publicos')

grafica.panel.fnc(x.sector$Catholic, vd='MathAch',

que.factor='SES',

x.panel='School',

regresion=T,

titulo='Privados')

Tal y como habíamos constatado anteriormente, parece que la pendiente de regresión es mayor en los colegios públicos que en los privados (católicos). Es decir parece que en los primeros el rendimiento en matemáticas está mas vinculado al nivel socioeconómico de los alumnos que en los segundos.

Seguidamente nos plantearemos un modelo jerárquico multinivel donde el primer nivel jerárquico serían los alumnos y el segundo los colegios a los que estos alumnos asisten.

En primer lugar debemos centrar la variable nivel socioconómico SES. Recuerda que en todo modelo de regresión el intercepto es el valor que adopta la variable criterio y cuando la predictora X vale 0. Debemos por tanto hacer que el modelo de regresión a estimar pase por cero introduciendo la variable nivel socioeconómico centrada (puntuaciones diferenciales). El centrado puede realizarse por la media global (gran media) o por la media de cada colegio.

Centrado por la media de grupo (Group mean centering)

alumnos=centra.variable.fnc(alumnos, variable='SES', que.factor='School')

Centrado por la gran media (Grand mean centering)

alumnos= centra.variable.fnc(alumnos, variable='SES')

Utilizaremos el centrado por la gran media porque el intercepto puede entonces ser interpretado como el valor de rendimiento ajustado por la media global.

#------------------------------------------------------------------

# CENTRADO DE VARIABLE

#------------------------------------------------------------------

*** Se ha centrado con exito la variable SES .El centrado

*** se ha realizado sobre la gran media. Si deseas un centrado

*** en cada nivel de un determinado factor, incluye el argumento

*** que.factor. Ej: que.factor='School'

Para estimar el modelo utilizaremos nuevamente la función mixed.model.fnc aunque habremos obligatoriamente de definir previamente el modelo que deseamos estimar.

model1='MathAch ~ c.SES + (c.SES | School)'

mixed.model.fnc(alumnos, modelo=model1)

Podemos definir además un modelo alternativo que queremos sea contrastado con el propuesto en la misma llamada a la función.

alt='MathAch ~ c.SES + (1 | School)'

mixed.model.fnc(alumnos, modelo=model1, alternativo=alt)

Sobre la base de datos alumnos, queremos estimar un modelo jerárquico multinivel con la variable MathAch como variable dependiente y la variable c.SES como independiente. Esta variable la asumimos con pendiente aleatoria en el segundo nivel jerárquico (una pendiente diferente por colegio) al definirla de esta manera (c.SES |School). El segundo nivel jerárquico (obligatorio) está definido al incluir la variable School a la derecha del caracter barra (|).

En el modelo alternativo planteamos la estimación del efecto fijo del nivel socioeconómico y suponemos una pendiente (efecto) común de esta variable para todos los sujetos.

Un supuesto muy importante de los modelos multinivel tiene que ver con la distribución normal de los factores aleatorios. En este caso podemos ver en el QQplot de los valores aleatorios estimados para el intercepto y la pendiente de regresión se encuentran dentro del intervalo de confianza para dicha normalidad.

Subir ->