Binomial glmer () singular a pesar de "muchos" datos
He estado leyendo diferentes preguntas sobre lo fácil que es toparse con singularidades al ajustar modelos de efectos mixtos con glmer(). En general, la idea es que las singularidades pueden surgir de estructuras aleatorias muy complejas. Si la estructura aleatoria es simple, también podría suceder cuando los datos no sean suficientes para calcular la matriz de varianza-covarianza ... vea, por ejemplo, esta página de Ben Bolker , la respuesta de Robert Long a esta publicación o la página de ayuda deisSingular() .
Sin embargo, el modelo que intento ajustar es muy simple:
mod.detection_rand <- glmer(reaction ~ Pedra + (1|Channel), family="binomial", data = garotes)
boundary (singular) fit: see ?isSingular
... y aparentemente tengo suficientes datos para las diferentes combinaciones de variables predictoras (fijas y aleatorias):
library(tidyverse)
garotes %>%
group_by(Channel, Pedra) %>%
summarise(n = n())
# A tibble: 16 x 3
# Groups: Channel [8]
Channel Pedra n
<int> <fct> <int>
1 1 No 13
2 1 Yes 13
3 2 No 14
4 2 Yes 12
5 3 No 12
6 3 Yes 14
7 4 No 13
8 4 Yes 13
9 5 No 13
10 5 Yes 13
11 6 No 14
12 6 Yes 12
13 7 No 13
14 7 Yes 13
15 8 No 14
16 8 Yes 12
¿Qué piensas?
EDITAR: Aquí está el resumen del modelo, summary(mod.detection_rand)
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: binomial ( logit )
Formula: reaction ~ Pedra + (1 | Channel)
Data: garotes
AIC BIC logLik deviance df.resid
261.5 271.5 -127.7 255.5 205
Scaled residuals:
Min 1Q Median 3Q Max
-1.8533 -0.9449 0.5396 0.5396 1.0583
Random effects:
Groups Name Variance Std.Dev.
Channel (Intercept) 0 0
Number of obs: 208, groups: Channel, 8
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.1133 0.1946 -0.582 0.56
PedraYes 1.3473 0.3066 4.394 1.11e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr)
PedraYes -0.635
convergence code: 0
boundary (singular) fit: see ?isSingular
EDIT2: Siguiendo el comentario de Billy:
bobyqa : boundary (singular) fit: see ?isSingular
[OK]
Nelder_Mead : boundary (singular) fit: see ?isSingular
[OK]
nlminbwrap : boundary (singular) fit: see ?isSingular
[OK]
nmkbw : boundary (singular) fit: see ?isSingular
[OK]
optimx.L-BFGS-B : boundary (singular) fit: see ?isSingular
[OK]
nloptwrap.NLOPT_LN_NELDERMEAD : boundary (singular) fit: see ?isSingular
[OK]
nloptwrap.NLOPT_LN_BOBYQA : boundary (singular) fit: see ?isSingular
[OK]
EDIT3: Siguiendo la respuesta de Isabella:
Verifiqué la estructura de la variable de resultado ( reaction). Aquí está la tabla de resultados:
library(tidyverse)
garotes %>%
group_by(Channel, Pedra, reaction) %>%
summarise(n = n()) %>%
print(n = Inf)
# A tibble: 32 x 4
# Groups: Channel, Pedra [16]
Channel Pedra reaction n
<int> <fct> <int> <int>
1 1 No 0 6
2 1 No 1 7
3 1 Yes 0 3
4 1 Yes 1 10
5 2 No 0 7
6 2 No 1 7
7 2 Yes 0 2
8 2 Yes 1 10
9 3 No 0 8
10 3 No 1 4
11 3 Yes 0 6
12 3 Yes 1 8
13 4 No 0 7
14 4 No 1 6
15 4 Yes 0 3
16 4 Yes 1 10
17 5 No 0 8
18 5 No 1 5
19 5 Yes 0 1
20 5 Yes 1 12
21 6 No 0 6
22 6 No 1 8
23 6 Yes 0 2
24 6 Yes 1 10
25 7 No 0 6
26 7 No 1 7
27 7 Yes 0 2
28 7 Yes 1 11
29 8 No 0 8
30 8 No 1 6
31 8 Yes 0 4
32 8 Yes 1 8
Aparentemente, existen ambos tipos de resultados para todos Channelsy todos los Pedratratamientos ... así que no es como el ejemplo que presentó Isabella ... además, intenté modelar este GLMM con el library(GLMMadaptive)y tampoco convergió.
EDIT4: El conjunto de datos que estoy usando, en caso de que alguien tenga curiosidad.
Channel Pedra reaction
1 No 1
2 No 0
3 No 0
4 No 0
5 No 0
6 No 1
7 No 0
8 No 0
1 No 1
2 No 1
3 No 1
4 No 1
5 No 0
6 No 0
7 No 0
8 No 0
1 No 0
2 No 1
3 No 0
4 No 0
5 No 0
6 No 0
7 No 0
8 No 1
1 No 0
2 No 1
3 Yes 0
4 Yes 1
5 Yes 1
6 Yes 1
7 Yes 1
8 Yes 0
1 Yes 1
2 Yes 1
3 Yes 0
4 Yes 0
5 No 0
6 No 1
7 Yes 1
8 Yes 1
1 Yes 0
2 Yes 1
3 Yes 1
4 Yes 1
5 Yes 1
6 Yes 0
7 No 1
8 No 1
1 Yes 1
2 Yes 1
3 Yes 1
4 Yes 1
5 Yes 1
6 Yes 1
7 Yes 1
8 Yes 1
1 Yes 1
2 Yes 1
3 Yes 1
4 Yes 1
5 Yes 0
6 Yes 1
7 Yes 1
8 Yes 1
1 Yes 1
2 Yes 1
3 Yes 0
4 Yes 1
5 Yes 1
6 Yes 1
7 Yes 0
8 Yes 0
1 Yes 1
2 Yes 1
3 Yes 0
4 Yes 0
5 Yes 1
6 Yes 1
7 Yes 1
8 Yes 0
1 Yes 1
2 Yes 1
3 Yes 0
4 Yes 1
5 Yes 1
6 Yes 1
7 Yes 0
8 Yes 0
1 Yes 1
2 Yes 0
3 Yes 1
4 Yes 0
5 Yes 1
6 Yes 1
7 Yes 1
8 Yes 1
1 Yes 1
2 Yes 1
3 Yes 0
4 Yes 1
5 Yes 1
6 Yes 0
7 Yes 1
8 Yes 1
1 Yes 1
2 Yes 1
3 Yes 1
4 Yes 1
5 Yes 1
6 Yes 1
7 Yes 1
8 Yes 1
1 Yes 0
2 Yes 0
3 Yes 1
4 Yes 1
5 Yes 1
6 Yes 1
7 Yes 1
8 Yes 1
1 Yes 1
2 No 0
3 Yes 1
4 No 1
5 Yes 1
6 No 1
7 Yes 1
8 No 1
1 No 0
2 Yes 1
3 No 0
4 Yes 1
5 No 1
6 Yes 1
7 No 1
8 Yes 1
1 Yes 0
2 No 1
3 Yes 1
4 No 0
5 Yes 1
6 No 1
7 Yes 1
8 No 0
1 No 0
2 No 1
3 No 1
4 No 0
5 No 1
6 No 0
7 No 0
8 No 0
1 No 1
5 No 0
3 No 1
4 No 1
2 No 1
6 No 0
7 No 1
8 No 0
1 No 0
5 No 0
3 No 0
4 No 0
2 No 1
6 No 0
7 No 0
8 No 0
1 No 1
5 No 1
3 No 1
4 No 0
2 No 0
6 No 1
7 No 1
8 No 0
1 No 1
5 No 0
3 No 0
4 No 1
2 No 0
6 No 1
7 No 1
8 No 1
1 No 1
5 No 1
3 No 0
4 No 1
2 No 0
6 No 1
7 No 1
8 No 1
1 No 1
5 No 1
3 No 0
4 No 0
2 No 0
6 No 1
7 No 0
8 No 0
1 No 0
5 No 0
3 No 0
4 No 1
2 No 0
6 No 0
7 No 1
8 No 1
¡Muchas gracias por todas sus respuestas, en cualquier caso! ¡Aprendiendo mucho de ellos!
Respuestas
Un comentario adicional: eché un vistazo a sus datos y está claro, nuevamente, que no hay evidencia de una variación sistemática entre los diferentes canales. Es por esto que el modelo mixto estima que la varianza entre canales es$0$, haciendo que el modelo sea singular.
Puede ver esto en la figura siguiente, donde los errores estándar para casi todos los canales se superponen ...
... y puede confirmarlo mediante la descomposición ANOVA de un GLM de efectos fijos, mostrando que no hay un efecto principal significativo de Canal (p = .986).
m_fixed_effects = glm(cbind(n, total) ~ Pedra + factor(Channel),
data=positive, family=binomial)
car::Anova(m_fixed_effects)
# Analysis of Deviance Table (Type II tests)
#
# Response: cbind(n, total)
# LR Chisq Df Pr(>Chisq)
# Pedra 4.9148 1 0.02663 *
# factor(Channel) 1.3859 7 0.98600
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Código
library(tidyverse)
df = read.csv('/path/to/reaction.csv')
head(df)
# Channel Pedra reaction n
# 1 1 No 0 6
# 2 1 No 1 7
# 3 1 Yes 0 3
# 4 1 Yes 1 10
# 5 2 No 0 7
# 6 2 No 1 7
df = df %>%
group_by(Channel, Pedra) %>%
mutate(total = sum(n),
prop = n / total,
se = sqrt((prop * (1-prop)) / n)) %>%
ungroup()
positive = filter(df, reaction==1)
ggplot(positive, aes(Pedra, prop, group=Channel, color=factor(Channel))) +
geom_path(position = position_dodge(width=.1)) +
geom_point(position = position_dodge(width=.1)) +
stat_summary(fun.data=mean_se, group=1, color='black',
position = position_nudge(x=c(-.2, .2))) +
geom_linerange(mapping=aes(ymin=prop-se, ymax=prop+se),
position = position_dodge(width=.1)) +
geom_hline(linetype='dashed', yintercept=.5) +
coord_cartesian(ylim=c(0, 1)) +
labs(color='Channel', y='Proportion positive reactions',
caption='Error bars show SEM')
m_fixed_effects = glm(cbind(n, total) ~ Pedra + factor(Channel),
data=positive, family=binomial)
car::Anova(m_fixed_effects)
# Analysis of Deviance Table (Type II tests)
#
# Response: cbind(n, total)
# LR Chisq Df Pr(>Chisq)
# Pedra 4.9148 1 0.02663 *
# factor(Channel) 1.3859 7 0.98600
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Debido a que este es un modelo de regresión logística binaria de efectos mixtos, se supone que su variable de resultado es binaria con valores codificados como 0 o 1.
Lo que necesita investigar es si tiene suficientes 1 presentes en su variable de respuesta para un número suficiente de 'sujetos'. (En su caso, sujeto significa canal).
Aquí hay un ejemplo inventado que produce la misma advertencia que lo que recibió:
SubjectID <- rep(1:5, each = 3)
SubjectID
Outcome <- rep(0, 15)
Outcome[1] <- 1
Data <- data.frame(Outcome, SubjectID)
str(Data)
Data
library(lme4)
glmer(Outcome ~ 1 + (1|SubjectID), family="binomial", data = Data)
En este ejemplo, hay 5 sujetos, de modo que 4 de ellos tienen solo 0 valores de resultado y uno de ellos tiene valores de resultado que incluyen un valor único de 1. (Cada sujeto tiene 3 valores de resultado en total).
Incluso si le da a cada uno de los sujetos en este ejemplo inventado un valor de 1 para su primer valor de resultado, seguirá recibiendo el mismo mensaje de error al ajustar el modelo:
Outcome <- rep(0, 15)
Outcome[c(1, 4, 7, 10, 13)] <- 1
Sin embargo, si los 4 sujetos que inicialmente tenían solo 0 valores pueden mantener estos valores y el primer sujeto recibe dos valores de 1, el mensaje de error desaparece:
Outcome <- rep(0, 15)
Outcome[c(1,2)] <- 1
Una vez que comprenda mejor el patrón de los valores 0 y 1 para la variable de resultado entre los sujetos de estudio, la otra cosa que puede intentar es ajustar su modelo con la función mixed_model () del paquete GLMMadaptive en R.
Para el pequeño ejemplo proporcionado aquí, esta función se usaría así:
library(GLMMadaptive)
m <- mixed_model(fixed = Outcome ~ 1,
random = ~ 1 | SubjectID,
data = Data,
family = binomial())
summary(m)
Isabella hizo algunos puntos excelentes. Esto también puede suceder cuando hay muy poca variación en el nivel del canal. Quizás los canales son muy similares entre sí, por lo que su varianza realmente es cercana a cero y, por lo tanto, no es necesaria en el modelo. Puede evaluar esto ajustando un glm y ver si las inferencias son similares.
Por curiosidad, ¿surge el error cuando utiliza un estimador alternativo? Podría ser que el estimador, por alguna razón, se atasque en una singularidad. Es posible que sólo intente lo siguiente: mod.alt_est <- allFit(mod.detection_rand). Alternativamente, es posible que necesite una solución bayesiana para ayudar a regularizar la estimación y alejarla de una singularidad (pruebe el blmepaquete si la allFitfunción no produce un estimador que funcione).