¿Cómo afectarán los efectos aleatorios con sólo 1 observación a un modelo lineal mixto generalizado?

Question

Tengo un conjunto de datos en el que la variable que me gustaría utilizar como efecto aleatorio sólo tiene una única observación para algunos niveles. Basándome en las respuestas a preguntas anteriores, he deducido que, en principio, esto puede estar bien.

¿Puedo ajustar un modelo mixto con sujetos que sólo tienen 1 observación?

Modelo de interceptos aleatorios - una medición por sujeto

Sin embargo, en el segundo enlace, la primera respuesta dice:

"...asumiendo que no estás usando un modelo mixto lineal generalizado GLMM donde en ese caso entran en juego cuestiones de dispersión excesiva"

Estoy considerando utilizar un GLMM, pero no entiendo muy bien cómo afectarán al modelo los niveles de efectos aleatorios con observaciones únicas.

---

He aquí un ejemplo de uno de los modelos que intento ajustar. Estoy estudiando las aves, y me gustaría modelar los efectos de la población y la temporada en el número de paradas durante la migración. Me gustaría utilizar el individuo como efecto aleatorio, porque para algunos individuos tengo hasta 5 años de datos.

library(dplyr)
library(lme4)
pop <- as.character(c("BF", "BF", "BF", "BF", "BF", "BF", "BF", "BF", "BF", "BF", "BF", "BF", "BF", "BF", "BF", "BF", "BF", "BF", "BF", "BF", "BF", "MA", "MA", "MA", "MA", "MA", "MA", "MA", "MA", "MA", "MA", "MA", "MA", "MA", "MA", "MA", "NU", "NU", "NU", "NU", "NU", "NU", "NU", "NU", "NU", "SA", "SA", "SA", "SA", "SA", "SA", "SA", "SA", "SA", "SA", "SA", "SA", "SA", "SA", "SA", "SA", "SA", "SA", "SA", "SA", "SA", "SA", "SA", "SA", "SA", "SA", "SA", "SA", "SA"))
id <- "2 2 4 4 7 7 9 9 10 10 84367 84367 84367 84368 84368 84368 84368 84368 84368 84369 84369 33073 33073 33073 33073 33073 33073 33073 33073 33073 80149 80149 80149 80150 80150 80150 57140 57141 126674 126677 126678 126680 137152 137152 137157 115925 115925 115925 115925 115925 115925 115925 115925 115926 115926 115926 115926 115926 115926 115927 115928 115929 115929 115929 115930 115930 115930 115930 115931 115931 115931 115932 115932 115932"
id <- strsplit(id, " ")
id <- as.numeric(unlist(id))
year <- "2014 2015 2014 2015 2014 2015 2014 2015 2014 2015 2009 2010 2010 2009 2010 2010 2011 2011 2012 2009 2010 2009 2009 2010 2010 2011 2011 2012 2012 2013 2008 2008 2009 2008 2008 2009 2008 2008 2013 2013 2013 2013 2014 2015 2014 2012 2013 2013 2014 2014 2015 2015 2016 2012 2013 2013 2014 2014 2015 2013 2012 2012 2013 2013 2012 2013 2013 2014 2013 2014 2014 2013 2014 2014"
year <- strsplit(year, " ")
year <- as.numeric(unlist(year))
season <- as.character(c("fall", "spring", "fall", "spring", "fall", "spring", "fall", "spring", "fall", "spring", "fall", "fall", "spring", "fall", "fall", "spring", "fall", "spring", "spring", "fall", "spring", "fall", "spring", "fall", "spring", "fall", "spring", "fall", "spring", "spring", "fall", "spring", "spring", "fall", "spring", "spring", "fall", "fall", "fall", "fall", "fall", "fall", "fall", "spring", "fall", "fall", "fall", "spring", "fall", "spring", "fall", "spring", "spring", "fall", "fall", "spring", "fall", "spring", "spring", "fall", "fall", "fall", "fall", "spring", "fall", "fall", "spring", "spring","fall", "fall", "spring", "fall", "fall", "spring"))
stops <- "0 0 0 0 0 0 1 0 2 1 1 0 0 3 2 0 1 1 0 1 1 2 0 1 0 2 0 4 0 0 2 1 1 2 5 2 1 0 9 6 2 3 4 7 2 0 0 0 0 0 2 0 0 1 0 0 0 0 0 0 1 1 0 0 1 1 0 0 1 1 0 0 0 0"
stops <- strsplit(stops, " ")
stops <- as.numeric(unlist(stops))

stopdata <- data.frame(pop = pop, id = id, year = year, season = season, stops = stops, stringsAsFactors = FALSE)

stopdata <- group_by(stopdata, pop, id)
summary1 <- summarise(stopdata, n.years = length(year))
table(summary1$n.years)

Hay 27 individuos. 9 individuos tienen una sola observación. 18 individuos tienen de 2 a 9 observaciones.

¿De qué hay que preocuparse si 1/3 de los niveles de efectos aleatorios sólo tienen una observación?

---

Lo he estado considerando:

Opción 1: GLMM como se ha descrito anteriormente

stops.glmm <- glmer(stops ~ pop + season + (1|id), data=stopdata, family = poisson)

Opción 2: Modelo lineal generalizado ponderado GLM utilizando las medias de los individuos con observaciones múltiples

aggfun <- function(data, idvars=c("pop", "season", "id"), response){
#select id variables, response variable, and year
sub1 <- na.omit(data[,c(idvars, "year", response)])
#aggregate for mean response by year
agg1 <- aggregate(sub1[names(sub1) == response],by=sub1[idvars],FUN=mean)
#sample size for each aggregated group
aggn <- aggregate(sub1[response],by=sub1[idvars],FUN=length)
#rename sample size column
names(aggn)[4] <- "n"
agg2 <- merge(agg1, aggn)
agg2}

#Create weighted dataset
stops.weight <- aggfun(data = stopdata, response = "stops")
stops.weight$stops <- round(stops.weight$stops)

#Weighted GLM
stops.glm <- glm(stops~pop + season, data=stops.weight, family = poisson, weights = n)

Answer 1

En general, tiene un problema de identificabilidad. Los modelos lineales con un efecto aleatorio asignado a un parámetro con una sola medición no pueden distinguir entre el efecto aleatorio y el error residual.

Una ecuación lineal típica de efecto mixto tendrá el siguiente aspecto:

$E = \beta + \eta_i + \epsilon_j$

Dónde $\beta$ es el efecto fijo, $\eta_i$ es el efecto aleatorio para el nivel $i$ y $\epsilon_j$ es la variabilidad residual del $j$ ª medición. Cuando sólo se tiene una observación de un nivel con un efecto aleatorio, es difícil distinguir entre $\eta$ y $\epsilon$ . Normalmente, ajustará una varianza o una desviación típica a $\eta$ y $\epsilon$ Por lo tanto, con sólo una medida por individuo, no estará tan seguro de tener una estimación precisa para $SD(\eta)$ y $SD(\epsilon)$ pero la estimación de la suma de las varianzas ( $var(\eta) + var(\epsilon)$ ) debería ser relativamente robusto.

Paso a la respuesta práctica: Si aproximadamente 1/3 de sus observaciones tienen una única observación por individuo, probablemente esté bien en general. El resto de la población debería proporcionar una estimación razonable para $SD(\eta)$ y $SD(\epsilon)$ y estos individuos deberían ser contribuyentes menores en general. Por otro lado, si tiene a todos los individuos en un efecto fijo y un efecto aleatorio específicos con una única medida (por ejemplo, para su ejemplo, tal vez toda una población tal vez eso signifique especies para usted), entonces confiaría menos en el resultado.

Answer 2

En este ejemplo se utilizó la distribución poisson, que por defecto es link log, ¿no? Yo esperaría que el efecto aleatorio, incluso si TODOS los grupos tuvieran una observación, para estimar la varianza en el registro, mientras que el error residual sería entonces mediría el error en la predicción incluyendo el efecto aleatorio.

Ejemplo, 10 observaciones, la primera sin efecto aleatorio:

 summary(pp <- glmmTMB(nvn~offset(log(vnhr)),data=foc10,family="poisson"))
 Family: poisson  ( log )
 Formula:          nvn ~ offset(log(vnhr))
 Data: foc10
      AIC      BIC   logLik deviance df.resid
    111.5    111.8    -54.7    109.5        9
 Conditional model:
              Estimate Std. Error z value Pr(>|z|)   
  (Intercept) -1.14680    0.09129  -12.56   <2e-16 ***

y ahora con efecto aleatorio

> summary(pp <- glmmTMB(nvn~offset(log(vnhr))+(1|focal),data=foc10,family="poisson"))
 Family: poisson  ( log )
Formula:          nvn ~ offset(log(vnhr)) + (1 | focal)
Data: foc10
     AIC      BIC   logLik deviance df.resid 
    54.1     54.7    -25.0     50.1        8 

[ observe un AIC mucho más bajo ]

 Random effects:
 Conditional model:
 Groups Name        Variance Std.Dev.
 focal  (Intercept) 5.088    2.256   
Number of obs: 10, groups:  focal, 10
Conditional model:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)   -2.692      1.166  -2.309   0.0209 *

En cualquier caso, no hay quejas de convergencia.
Mientras que los hay si dejo de lado el poisson.