¿Cuál es el número mínimo recomendado de grupos para un factor de efectos aleatorios?

Question

Estoy utilizando un modelo mixto en R ( lme4 ) para analizar algunos datos de medidas repetidas. Tengo una variable de respuesta (contenido de fibra en las heces) y 3 efectos fijos (masa corporal, etc.). Mi estudio sólo tiene 6 participantes, con 16 medidas repetidas para cada uno (aunque dos sólo tienen 12 repeticiones). Los sujetos son lagartos a los que se les dio diferentes combinaciones de alimentos en diferentes "tratamientos".

Mi pregunta es: ¿puedo utilizar el ID del sujeto como efecto aleatorio?

Sé que esta es la forma habitual de actuar en los modelos longitudinales de efectos mixtos, para tener en cuenta la naturaleza de la muestra aleatoria de los sujetos y el hecho de que las observaciones dentro de los sujetos estarán más correlacionadas que las de los sujetos. Pero, tratar el ID del sujeto como un efecto aleatorio implica estimar una media y una varianza para esta variable.

• Dado que sólo tengo 6 sujetos (6 niveles de este factor), ¿es esto suficiente para obtener una caracterización precisa de la media y la varianza?

• ¿Ayuda en este sentido el hecho de que tenga bastantes mediciones repetidas para cada sujeto (no veo qué importancia tiene)?

• Por último, si no puedo utilizar el ID del sujeto como efecto aleatorio, ¿incluirlo como efecto fijo me permitirá controlar el hecho de que tengo medidas repetidas?

Editar: Me gustaría aclarar que cuando digo "puedo" utilizar la identificación del sujeto como efecto aleatorio, quiero decir "es una buena idea hacerlo". Sé que puedo ajustar el modelo con un factor de sólo 2 niveles, pero seguramente esto sería indefendible. Pregunto en qué momento es sensato pensar en tratar a los sujetos como efectos aleatorios. Parece que la literatura aconseja que 5-6 niveles es un límite inferior. Me parece que las estimaciones de la media y la varianza del efecto aleatorio no serían muy precisas hasta que hubiera más de 15 niveles de factores.

Answer 1

En un intento de averiguar el número mínimo de grupos para un modelo multinivel, consulté el libro Análisis de datos mediante modelos de regresión y multinivel/jerárquicos de Gelman y Hill (2007).

Parece que abordan este tema en el capítulo 11, sección 5 (página 247), donde escriben que cuando hay < 5 grupos, los modelos multinivel suelen aportar poco respecto a los modelos clásicos. Sin embargo, parece que escriben que hay poco riesgo al aplicar un modelo multinivel.

Los mismos autores parecen retomar este tema en el capítulo 12, sección 9 (páginas 275-276). Allí escriben que los consejos sobre el número mínimo de grupos para un modelo multinivel son erróneos. Allí vuelven a decir que los modelos multinivel suelen añadir poco respecto a los modelos clásicos cuando el número de grupos es pequeño. Sin embargo, también escriben que los modelos multinivel no deberían ser peores que la regresión sin grupos (donde sin grupos parece significar que los indicadores de grupo se utilizan en la regresión clásica).

En las páginas 275-276 los autores tienen una subsección específica para el caso de uno o dos grupos (por ejemplo, hombres frente a mujeres). Aquí escriben que suelen expresar el modelo en forma clásica. Sin embargo, afirman que la modelización multinivel puede ser útil incluso con uno o dos grupos. Escriben que con uno o dos grupos la modelización multinivel se reduce a la regresión clásica.

Mi impresión es que la regresión clásica es un extremo de un continuo de modelos, es decir, un caso especial de un modelo multinivel.

Basándome en lo anterior, mi impresión es que la regresión clásica y el modelo multinivel devolverán estimaciones casi idénticas cuando sólo hay dos grupos y que utilizar modelos multinivel con sólo uno, dos, tres, cuatro, cinco o seis grupos es correcto.

Intentaré modificar esta respuesta en el futuro con R y un pequeño conjunto de datos que comparan las estimaciones obtenidas con ambos enfoques cuando se utilizan dos grupos.

Answer 2

Respuesta corta: Sí, puede utilizar la ID como efecto aleatorio con 6 niveles.

Respuesta un poco más larga: El Preguntas frecuentes sobre la GLMM de @BenBolker dice (entre otras cosas) lo siguiente bajo el título " ¿Debo tratar el factor xxx como fijo o aleatorio? ":

Un punto de especial relevancia para la estimación de modelos mixtos "modernos (en lugar de la estimación "clásica" del método de los momentos) es que, para para fines prácticos, debe haber un número razonable de niveles de efectos aleatorios (por ejemplo, bloques), más de 5 o 6 como mínimo.

Así que estás en el límite inferior, pero en el lado derecho del mismo.

Answer 3

Por si sirve de algo, hice un pequeño estudio de simulación para ver la estabilidad de la estimación de la varianza para un LMM relativamente simple (usando el sleepstudy conjunto de datos disponible a través de lme4 ). La primera forma genera todas las combinaciones posibles de temas para ngroups número de sujetos, y reajusta el modelo para cada combinación posible. La segunda toma varios subconjuntos aleatorios de sujetos.

library(lme4)
library(ggplot2)
library(tidyr)

m0 <- lmer(Reaction ~ Days + (1|Subject), data = sleepstudy,
           control = lmerControl(optimizer = "nloptwrap"))
# set the number of factor levels
ngroups <- 3:18 
# generate all possible combinations
combos <- lapply(X = ngroups, 
                 FUN = function(x) combn(unique(sleepstudy$Subject), x)) 

# allocate output (sorry, this code is entirely un-optimized)
out <- list(matrix(NA, ncol(combos[[1]]), 1), matrix(NA, ncol(combos[[2]]), 1),
            matrix(NA, ncol(combos[[3]]), 1), matrix(NA, ncol(combos[[4]]), 1),
            matrix(NA, ncol(combos[[5]]), 1), matrix(NA, ncol(combos[[6]]), 1),
            matrix(NA, ncol(combos[[7]]), 1), matrix(NA, ncol(combos[[8]]), 1),
            matrix(NA, ncol(combos[[9]]), 1), matrix(NA, ncol(combos[[10]]), 1),
            matrix(NA, ncol(combos[[11]]), 1), matrix(NA, ncol(combos[[12]]), 1),
            matrix(NA, ncol(combos[[13]]), 1), matrix(NA, ncol(combos[[14]]), 1),
            matrix(NA, ncol(combos[[15]]), 1), matrix(NA, ncol(combos[[16]]), 1))
# took ~ 2.5 hrs on my laptop, commented out for safety
#system.time(for(ii in 1:length(combos)) {
#    for(jj in 1:ncol(combos[[ii]])) {
#    sls <- sleepstudy[sleepstudy$Subject %in% combos[[ii]][,jj],]
#    out[[ii]][jj] <- attr(VarCorr(update(m0, data = sls))$Subject, 'stddev')
#        }
#    })

# pad with zeros, not all were equal
# from http://stackoverflow.com/questions/11148429/r-convert-asymmetric-list-to-matrix-number-of-elements-in-each-sub-list-diffe
max.len <- max(sapply(out, length))
corrected.list <- lapply(out, function(x) {c(x, rep(NA, max.len - length(x)))})
mat <- do.call(rbind, corrected.list)
mat <- data.frame(t(mat))
names(mat) <- paste0('s',3:18)
mat <- gather(mat, run, value)

ggplot(mat, aes(x = value, fill = run)) + 
    geom_histogram(bins = 60) +
    geom_vline(xintercept = 37.12, linetype =  'longdash', 
               aes(colour = 'original')) +
    facet_wrap(~run, scales = 'free_y') +
    scale_x_continuous(breaks = seq(0, 100, by = 20)) + 
    theme_bw() + 
    guides(fill = FALSE)

La línea negra punteada es la estimación puntual original de la varianza, y las facetas representan diferentes números de sujetos ( s3 siendo grupos de tres sujetos, s4 ser cuatro, etc.).

Y la forma alternativa:

ngroups <- 3:18
reps <- 500
out2<- matrix(NA, length(ngroups), reps)

for (ii in 1:length(ngroups)) {
    for(j in 1:reps) {
        sls <- sleepstudy[sleepstudy$Subject %in% sample(unique(sleepstudy$Subject), ngroups[i], replace = FALSE),]
        out2[i,j] <- attr(VarCorr(update(m0, data = sls))$Subject, 'stddev')
    }
}
out2 <- data.frame(t(out2))
names(out2) <- paste0('s',3:18)
out2 <- gather(out2, run, value)

ggplot(out2, aes(x = value, fill = run)) + 
    geom_histogram(bins = 60) +
    geom_vline(xintercept = 37.12, linetype =  'longdash', 
               aes(colour = 'original')) +
    facet_wrap(~run, scales = 'free_y') +
    scale_x_continuous(breaks = seq(0, 100, by = 20)) + 
    theme_bw() + 
    guides(fill = FALSE)

Parece (para este ejemplo, al menos) que la varianza no se estabiliza realmente hasta que hay al menos 14 sujetos, si no más tarde.

Answer 4

También podría utilizar un modelo mixto bayesiano, en cuyo caso la incertidumbre en la estimación de los efectos aleatorios se tiene en cuenta en el cálculo de los intervalos de credibilidad de la predicción del 95%. El nuevo paquete R brms y la función brm por ejemplo, permiten una transición muy fácil de un lme4 modelo mixto frecuentista a uno bayesiano, ya que tiene una sintaxis casi idéntica.

Answer 5

La "econometría más inofensiva" de Angrist y Pischke tiene una sección titulada "Menos de 42 racimos", en la que dicen en tono de broma,

Por lo tanto, siguiendo el... dictamen de que la respuesta a la vida, el universo y todo lo demás es 42, creemos que la pregunta es: ¿Cuántos conglomerados son suficientes para una inferencia fiable utilizando el ajuste de conglomerados estándar [similar al estimador de la varianza en GEE]?

La forma en que mi instructor de econometría solía responder a preguntas como la tuya es: "Estados Unidos es un país libre, puedes hacer lo que quieras. Pero si quieres que te publiquen tu trabajo, tienes que ser capaz de defender lo que has hecho". En otras palabras, es probable que puedas ejecutar el código de R o Stata o HLM o Mplus o SAS PROC GLIMMIX con 6 sujetos (y cambiar a estos paquetes alternativos si el de tu elección no lo ejecuta), pero probablemente te será muy difícil defender este enfoque y justificar las pruebas asintóticas.

Creo que, por defecto, incluir una variable como pendiente aleatoria implica incluirla también como efecto fijo, y hay que pasar por un montón de aros de sintaxis si sólo se quiere tener esto como efecto aleatorio con la media de cero. Esta es una elección sensata que los desarrolladores del software han hecho por usted.