Introducción a modelos de coeficientes aleatorios (RCM)
RCM se refiere a una técnica de análisis para modelar datos agrupados u organizados jerárquicamente donde hay una única variable dependiente medida al nivel más bajo, y variables predictoras en todos los niveles.
Diagnosticando Depresión: Una Ilustración
Un investigador quiere saber si el sexo del evaluador influye en el puntaje que le asignan a sus clientes en una escala de depresión. Los datos se presentan a continuación.
La variable predictora (sexo) tiene dos niveles (Hombre = 1; Mujer = 0). En este caso, hay 6 evaluadores (3 hombres y 3 mujeres) pero no vamos a considerar esa información por el momento. A cada evaluador se le asignaron aleatoriamente 4 clientes.
Podemos evaluar la hipótesis nula de que hombres y mujeres asignan (en promedio) los mismos puntajes usando una prueba t o ANOVA. En ambos casos, el análisis sugiere que la diferencia entre ambos grupos es estadísticamente significativa
#Diferencia de medias usando una prueba t
t.test(rating~sex,data=d1,var.equal=TRUE)
##
## Two Sample t-test
##
## data: rating by sex
## t = 4.4749, df = 22, p-value = 0.0001891
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 7.511767 20.488233
## sample estimates:
## mean in group 0 mean in group 1
## 60 46
#Diferencia de medias usando ANOVA
aov1<-aov(rating~sex,data=d1)
summary(aov1)
## Df Sum Sq Mean Sq F value Pr(>F)
## sex 1 1176 1176.0 20.02 0.000189 ***
## Residuals 22 1292 58.7
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Sin embargo, este análisis no considera los efectos del evaluador (trainee). Los clientes en este caso están anidados en los evaluadores, lo que representan una fuente de varianza que debemos considerar si queremos extrapolar los resultados a otros evaluadores. Esto es, los evaluadores son un factor aleatorio ya que podríamos haber tomado un grupo distintos de personas para llevar a cabo las evaluaciones, y en cada réplica los resultados se van a ver afectados no sólo por las variaciones en las medias muestrales de hombres y mujeres (efectos fijos) sino también por variaciones entre evaluadores.
Modelo Nivel 1: \[Y_{ij}=\beta_{0_j}+e_{ij}\] Modelo Nivel 2: \[\beta_{0_j}=\gamma_{00}+u_{0j}\] Modelo combinado: \[Y_{ij}=\gamma_{00}+u_{0j}+e_{ij}\] Lo que quiere decir este modelo es que el puntaje del cliente i evaluado por el evaluador j corresponde a la gran media, \(\gamma_{00}\), más el efecto del evaluador, \(U_{0j}\), más un error, \(e_{ij}\).
En la Figura de más abajo se observa que los evaluadores (trainees) asignan puntajes distintos a sus clientes.
## Loading required package: ggplot2
## Registered S3 methods overwritten by 'ggplot2':
## method from
## [.quosures rlang
## c.quosures rlang
## print.quosures rlang
Como un paso preparatorio inicial vamos simplemente a determinar si existen diferencias entre evaluadores que sea pertinente incluir en el modelo.
Random intercept model (One-way ANOVA with random effect)
¿Qué proporción de la varianza de los puntajes de los clientes está asociada con diferencias entre evaluadores?
Para el caso de un evaluador podemos generar un modelo usando una regresión. Como se ve a continuación el intercepto para los evaluadores 1 y 2 es 40 y 50, respectivamente. El intercepto para cada uno es promedio de los puntajes que asignarion a sus clientes.
#Para el evaluador 1
lm.e1<-lm(rating~1, data=subset(d1,trainee==1))
summary(lm.e1)
Vovamos a la fórmula de nuestro modelo combinado: \[Y_{ij}=\gamma_{00}+u_{0j}+e_{ij}\] A partir de esta fórmula podemos determinar qué proporción de la varianza de Y depende del efecto del evaluador (i.e., la variable de agrupación). La varianza de \(u_{0j}\) corresponde a la varianza de los evaluadores, y se designa mediante \(\sigma^2_{u0}\) o \(\tau_{00}\). La varianza de \(e_{ij}\) corresponde a la varianza residual, y se designa mediante \(\sigma^2_e\) o simplemente \(\sigma^2\).
La varianza total del modelo corresponde aproximadamente a la suma de la varianza entre evaluadores y la varianza residual, 1 El valor preciso de estos componentes se obtiene mediante maximum likelihood, no mediante OLS
\[SS_T=SS_B+SS_W\approx\tau_{00}+\sigma^2_e\]
Correlación intraclase (ICC)
A partir de la decomposición de la varianza en el paso anterior podemos estimar qué proporción de la varianza total está asociada a la variable de agrupación (i.e., evaluadores) mediante:
\[\rho=\frac{\tau_{00}}{\tau_{00}+\sigma^2_e}\]
El valor de ICC se puede interpretar como la correlación promedio entre cualquier par observaciones (clientes) dentro de un mismo grupo (evaluador). En el caso extremo, si ICC = 1.0 todos las observaciones dentro de un grupo tienen el mismo valor.
Efecto del diseño
Es posible evaluar si es necesario usar un análisis multinivel mediante la siguiente fórmula (Muthen & Satorra, 1995):
\[Design effect = 1+(\overline{n}-1)*ICC\] donde \(\overline{n}\) corresponde al tamaño promedio de los grupos. Valores menores que 2 indican que un diseño uninivel no produciría resultados sesgados.
Estimación mediante lme4
Usaremos el paquete lme4 para estimar el modelo inicial. Como en este paso sólo nos interesa modelar el intercepto usaremos la fórmula rating~1 para indicar que no se incluyen variables fijas. Sí incluiremos el factor aleatorio evaluador (trainee). Por último, el método de estimación será restricted maximum likelihood o REML.
## Linear mixed model fit by REML ['lmerMod']
## Formula: rating ~ 1 + (1 | trainee)
## Data: d1
##
## REML criterion at convergence: 166.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.59965 -0.49688 -0.01505 0.60318 1.50207
##
## Random effects:
## Groups Name Variance Std.Dev.
## trainee (Intercept) 71.01 8.427
## Residual 45.56 6.749
## Number of obs: 24, groups: trainee, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 53.000 3.706 14.3
De acuerdo a los resultados del modelo, la varianza explicada por el evaluador (trainee) es 71.01 y la varianza residual es 45.56. La gran media es 53, que corresponde al intercepto del grupo completo.
Ahora podemos obtener el valor para ICC mediante:
\[\rho=\frac{\tau_{00}}{\tau_{00}+\sigma^2_e}=\frac{71.01}{71.01+45.56}=0.61\] El valor de ICC es bastante alto. El efecto del diseño en este caso es:
\[Design effect = 1+(\overline{n}-1)\times ICC\]
\[Design effect = 1+(4-1)\times 0.61=2.44\]
Means-as-outcoms Model
¿Existe una diferencia entre los puntajes que hombres y mujeres le asignan a sus clientes en la escala de depresión?
Ahora que hemos establecido que la variable de agrupación es relevante, procederemos a evaluar la hipótesis sustantiva del estudio mediante un modelo multinivel incluyendo la variable sexo como factor en el Nivel 2 (También podríamos haber incluido variable adicionales en el Nivel 1, como la edad de los clientes, pero veremos eso más adelante)
Modelo Nivel 1: \[Y_{ij}=\beta_{0_j}+e_{ij}\]
Modelo Nivel 2: \[\beta_{0_j}=\gamma_{00}+\gamma_{01}Sexo+u_{0j}\]
## Linear mixed model fit by REML ['lmerMod']
## Formula: rating ~ sex + (1 | trainee)
## Data: d1
##
## REML criterion at convergence: 155.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8248 -0.4452 -0.1247 0.6843 1.4141
##
## Random effects:
## Groups Name Variance Std.Dev.
## trainee (Intercept) 18.11 4.256
## Residual 45.56 6.749
## Number of obs: 24, groups: trainee, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 60.000 3.136 19.134
## sex1 -14.000 4.435 -3.157
##
## Correlation of Fixed Effects:
## (Intr)
## sex1 -0.707
El efecto de la variable sexo en el puntaje asignado es -14, lo que quiere decir que, en promedio, los hombres le asignan 14 puntos menos a sus clientes que las mujeres. El valor t asociado es de -3.157 cuyo valor p obviamente será significativo pero la función lmer no entrega valores p exactos.
Para obtener el valor p exacto debemos instalar una extensión del paquete lme4 que se llama lmerTest.
require(lmerTest)
## Loading required package: lmerTest
##
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
##
## lmer
## The following object is masked from 'package:stats':
##
## step
Una vez instalado el paquete lmerTest simplemente volvemos a correr el modelo anterior y podemos ver ahora que el valor p es 0.0343.
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: rating ~ sex + (1 | trainee)
## Data: d1
##
## REML criterion at convergence: 155.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8248 -0.4452 -0.1247 0.6843 1.4141
##
## Random effects:
## Groups Name Variance Std.Dev.
## trainee (Intercept) 18.11 4.256
## Residual 45.56 6.749
## Number of obs: 24, groups: trainee, 6
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 60.000 3.136 4.000 19.134 4.4e-05 ***
## sex1 -14.000 4.435 4.000 -3.157 0.0343 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## sex1 -0.707
El valor p asociado a la variable sexo en el modelo multinivel (p=0.0343) es bastante mayor que el valor p del modelo inicial (p=0.0001891). Siempre que la variable de agrupación tenga alguna asociación con la variable dependiente (ICC>0) el modelo multinivel va a usar un error más grande como estándar de comparación (ya que incluye la varianza de la variable de agrupación). Esto es importante en el sentido que manteniendo todo constante el los modelos multinivel tienen MENOS poder que modelos que no incluyen la información de la variable de agrupación.
