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

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

Mi pregunta es: ¿puedo usar la identificación del sujeto como un efecto aleatorio?

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

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

• ¿El hecho de que tengo bastantes mediciones repetidas para cada sujeto ayuda en este sentido (no veo cómo importa)?

• Finalmente, si no puedo usar la identificación del sujeto como un efecto aleatorio, ¿incluirlo como un efecto fijo me permitirá controlar el hecho de que tengo medidas repetidas?

Editar: Me gustaría aclarar que cuando digo "puedo" usar la identificación del sujeto como un efecto aleatorio, quiero decir "es una buena idea". Sé que puedo ajustar el modelo con un factor con solo 2 niveles, pero ¿seguramente esto no sería defendible? Me pregunto en qué punto se vuelve razonable 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.

Respuesta corta: Sí, puedes usar ID como efecto aleatorio con 6 niveles.

Respuesta un poco más larga: las preguntas frecuentes de GLMM de @BenBolker dicen (entre otras cosas) lo siguiente bajo el título " ¿Debería tratar el factor xxx como fijo o aleatorio? ":

Un punto de particular relevancia para la estimación del modelo mixto 'moderno' (en lugar de la estimación del método de momentos 'clásico') es que, 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.

Entonces estás en el límite inferior, pero en el lado derecho.

En un intento por determinar el número mínimo de grupos para un modelo multinivel, miré el libro Análisis de datos usando regresión y modelos multinivel / jerárquico de Gelman y Hill (2007).

Parecen abordar este tema en el Capítulo 11, Sección 5 (página 247) donde escriben que cuando hay <5 grupos, los modelos multinivel generalmente agregan poco sobre los modelos clásicos. Sin embargo, parecen escribir que hay poco riesgo de aplicar un modelo multinivel.

Los mismos autores parecen volver a este tema en el Capítulo 12, Sección 9 (páginas 275-276). Allí escriben que el consejo sobre el número mínimo de grupos para un modelo multinivel es erróneo. Allí nuevamente dicen que los modelos multinivel a menudo agregan poco sobre 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 agrupación (donde la no agrupación 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 (p. Ej., Hombre versus mujer). Aquí escriben que típicamente expresan el modelo en forma clásica. Sin embargo, afirman que el modelado multinivel puede ser útil incluso con solo uno o dos grupos. Escriben que con uno o dos grupos el modelado multinivel se reduce a la regresión clásica.

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

Con base en lo anterior, mi impresión es que la regresión clásica y el modelado multinivel arrojarán estimaciones casi idénticas cuando solo hay dos grupos y que usar modelos multinivel con solo uno, dos, tres, cuatro, cinco o seis grupos está bien.

Intentaré modificar esta respuesta en el futuro con un Rcódigo y un pequeño conjunto de datos que comparen las estimaciones obtenidas con ambos enfoques al usar dos grupos.

Para lo que vale, hice un poco de estudio de simulación para observar la estabilidad de la estimación de la varianza para un LMM relativamente simple (utilizando el sleepstudyconjunto de datos disponible a través de lme4). La primera forma genera todas las combinaciones de sujetos posibles para el ngroupsnúmero de sujetos y ajusta el modelo para cada combinación posible. El segundo 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 ( s3grupos de tres sujetos, s4cuatro, 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, de todos modos) que la variación realmente no se estabiliza hasta que hay al menos 14 sujetos, si no más tarde.

La "Econometría en su mayoría inofensiva" de Angrist y Pischke tiene una sección titulada, "Menos de 42 grupos", 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 es 42, creemos que la pregunta es: ¿Cuántos grupos son suficientes para una inferencia confiable usando el ajuste de grupo estándar [similar al estimador de varianza en GEE]?

La forma en que mi instructor de econometría solía responder preguntas como la suya es: "Estados Unidos es un país libre, puedes hacer lo que quieras. Pero si quieres que se publique tu trabajo, debes ser capaz de defender lo que has hecho". " En otras palabras, es probable que pueda ejecutar R o Stata o HLM o Mplus o código SAS PROC GLIMMIX con 6 temas (y cambiar a estos paquetes alternativos si el que elija no ejecuta esto), pero es probable que tenga Es muy difícil defender este enfoque y justificar las pruebas asintóticas.

Creo que, por defecto, incluir una variable como una pendiente aleatoria implica incluir eso también como un efecto fijo, y debe saltar a través de muchos aros de sintaxis si solo desea tener esto como un efecto aleatorio con la media de cero. Esa es una elección sensata que los desarrolladores de software han hecho para usted.

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

No usaría un modelo de efectos aleatorios con solo 6 niveles. Los modelos que usan un efecto aleatorio de 6 niveles pueden ejecutarse en algún momento utilizando muchos programas estadísticos y, a veces, dan estimaciones imparciales, pero:

1. Creo que hay un consenso arbitrario en la comunidad estadística de que 10-20 es el número mínimo. Si desea que se publique su investigación, se le recomendará que busque una revista sin revisión estadística (o que pueda justificar su decisión con un lenguaje bastante sofisticado).
2. Con tan pocos grupos, es probable que la varianza entre grupos sea poco estimada. La mala estimación de la varianza entre grupos generalmente se traduce en una mala estimación del error estándar de los coeficientes de interés. (los modelos de efectos aleatorios se basan en el número de grupos que teóricamente van al infinito).
3. A menudo, los modelos simplemente no convergen. ¿Has intentado ejecutar tu modelo? Me sorprendería con solo 12-16 medidas por tema que los modelos converjan. Cuando logré que este tipo de modelo convergiera, tuve cientos de mediciones por grupo.

Este problema se aborda en la mayoría de los libros de texto estándar en el campo y usted los ha abordado en su pregunta. No creo que te esté dando ninguna información nueva.

Ha pasado mucho tiempo desde la pregunta original, pero pensé que podría agregar algunos puntos pertinentes a la selección del modelo.

1 - Siempre y cuando se identifique el modelo (es decir, tiene grados de libertad en el espacio de parámetros), debería ser capaz de INTENTAR para adaptarse al modelo. Dependiendo del método de optimización, el modelo puede o no converger. En cualquier caso, no trataría de incluir más de 1 o 2 efectos aleatorios y definitivamente no más de 1 interacción de nivel cruzado. En el caso específico del problema presentado aquí si sospechamos que una interacción entre las características específicas de la lagartija (por ejemplo, edad, tamaño, etc.) y las características de tratamiento / medida, el tamaño del grupo 6 puede no ser suficiente para hacer estimaciones lo suficientemente precisas.

2 - Como mencionan algunas respuestas, la convergencia puede ser un problema. Sin embargo, mi experiencia es que, si bien los datos de las ciencias sociales tienen un gran problema de convergencia debido a problemas de medición, las ciencias biológicas y especialmente las medidas bioquímicas repetidas tienen errores estándar mucho más pequeños. Todo depende del proceso de generación de datos. En los datos sociales y económicos tenemos que trabajar en varios niveles de abstracción. En biológicos y químicos, y ciertamente el error de medición de datos astronómicos es un problema menor.