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

26

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.

Chris
fuente

Respuestas:

21

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.

Henrik
fuente
12

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.

Mark Miller
fuente
10

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.). ingrese la descripción de la imagen aquí

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)

ingrese la descripción de la imagen aquí

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.

alexforrence
fuente
1
+1. Por supuesto, cuanto menor es el número de sujetos, mayor es la varianza del estimador de varianza. Pero no creo que esto sea lo que importa aquí. La pregunta es, ¿qué número de sujetos permite obtener algunos resultados razonables? Si definimos el resultado "no sensible" como la obtención de la varianza cero, entonces en su simulación ocurre con bastante frecuencia con n = 5 o menos. A partir de n = 6 o n = 7, casi nunca se obtiene una estimación exacta de la varianza 0, es decir, el modelo está convergiendo a una solución no degenerada. Mi conclusión sería que n = 6 es límite aceptable.
ameba dice Reinstate Monica
1
Por cierto, esto está en línea con rpubs.com/bbolker/4187 .
ameba dice Reinstate Monica
8

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.

StasK
fuente
1
Supongo que la respuesta a la pregunta es, en cierta medida, "cuánto dura un trozo de cuerda". Pero, no confiaría mucho en estimar una media o varianza de una muestra de menos de 15-20, por lo que no se aplicaría la misma regla general a los niveles de un efecto aleatorio. Nunca he visto a nadie incluir la identificación del sujeto como un efecto fijo y aleatorio en los estudios longitudinales. ¿Es esta una práctica común?
Chris
Además de un pequeño número de sujetos en el modelo mixto, sus efectos aleatorios no son observados, por lo que debe eliminarlos de los datos, y posiblemente necesite relativamente más datos para hacerlo de manera confiable que solo estimar la media y el varianza cuando todo se observa. Así 42 vs. 15-20 :). Creo que me refería a pendientes aleatorias, ya que estás en lo correcto en las identificaciones de los sujetos que aparecen solo como efectos aleatorios, de lo contrario no se identificarán. Los economistas no creen en los efectos aleatorios, por cierto, y publican casi exclusivamente lo que llaman "efectos fijos", es decir, estimaciones dentro del tema.
StasK
2
+1 @StasK para obtener una muy buena respuesta a una pregunta que es muy difícil de tratar. Sin embargo, creo que hay un tinte de sarcasmo innecesario y podrías considerar editar tu respuesta para ser un poco más respetuoso con el OP.
Michael R. Chernick
@Michael, probablemente tengas razón en que esta es una respuesta malhumorada, y probablemente innecesariamente. Sin embargo, el OP aceptó la respuesta que querían escuchar, por lo que obtuvo una resolución al respecto. Una respuesta más seria apuntaría a una buena evidencia de simulación o a un análisis asintótico de orden superior, pero desafortunadamente no conozco tales referencias.
StasK
3
Por lo que vale, creo que el número mágico "42" no se trata de cuándo se justifican los efectos aleatorios, sino de cuando uno puede escapar sin preocuparse por las correcciones de tamaño finito (por ejemplo, pensar en grados de libertad de denominador efectivo / correcciones de Kenward-Roger / otros enfoques similares).
Ben Bolker
7

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.

Tom Wenseleers
fuente
4

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.

Charles
fuente
¿Se rechazó esta votación por una razón relacionada con su contenido técnico?
N Brouwer
¿Con qué tipo de datos estás trabajando? No estoy seguro de por qué te sorprende escuchar que el modelo convergerá con 12-16 medidas por persona. No puedo comentar sobre el sesgo en los modelos resultantes, pero nunca he tenido problemas con la convergencia en lme4modelos mixtos y a menudo los ejecuto en tamaños de muestra similares a los del OP (también estoy trabajando con conjuntos de datos de biología).
RTbecard
1

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.

m_e_s
fuente