Binomial glmer () singulier malgré "beaucoup" de données
J'ai lu différentes questions sur la facilité avec laquelle il est possible de se heurter à des singularités lors de l'ajustement de modèles d'effets mixtes glmer(). En général, l'idée est que les singularités peuvent provenir de structures aléatoires très complexes. Si la structure aléatoire est simple, cela peut aussi arriver lorsque les données ne sont pas suffisantes pour calculer la matrice de variance-covariance ... voir par exemple cette page de Ben Bolker , la réponse de Robert Long à cet article ou la page d'aide deisSingular() .
Cependant, le modèle que j'essaie d'adapter est très simple:
mod.detection_rand <- glmer(reaction ~ Pedra + (1|Channel), family="binomial", data = garotes)
boundary (singular) fit: see ?isSingular
... et apparemment, j'ai suffisamment de données pour les différentes combinaisons de variables prédictives (fixes et aléatoires):
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'est-ce que tu penses?
EDIT: Voici le résumé du modèle, 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: Suite au commentaire 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: Suite à la réponse d'Isabella:
J'ai vérifié la structure de la variable de résultat ( reaction). Voici le tableau des résultats:
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
Apparemment, il y a les deux types de résultats pour tous Channelset tous les Pedratraitements ... donc ce n'est pas comme l'exemple présenté par Isabella ... de plus, j'ai essayé de modéliser ce GLMM avec le library(GLMMadaptive)et il n'a pas non plus convergé.
EDIT4: L'ensemble de données que j'utilise, au cas où quelqu'un serait curieux.
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
Merci beaucoup pour toutes vos réponses, en tout cas! Apprendre beaucoup d'eux!
Réponses
Un autre commentaire: j'ai jeté un coup d'œil à vos données, et il est clair, encore une fois, qu'il n'y a aucune preuve de variance systématique entre les différents canaux. C'est pourquoi le modèle mixte estime la variance inter-canaux à$0$, rendant le modèle singulier.
Vous pouvez le voir dans la figure ci-dessous, où les erreurs standard pour presque tous les canaux se chevauchent ...
... et peut le confirmer par décomposition ANOVA d'un GLM à effets fixes, montrant qu'il n'y a pas d'effet principal significatif de Channel (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
Code
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
Comme il s'agit d'un modèle de régression logistique binaire à effets mixtes, il suppose que votre variable de résultat est binaire avec des valeurs codées comme 0 ou 1.
Ce que vous devez rechercher, c'est si vous avez suffisamment de 1 dans votre variable de réponse pour un nombre suffisant de «sujets». (Dans votre cas, sujet signifie chaîne.)
Voici un exemple inventé qui produit le même avertissement que wnat que vous avez obtenu:
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)
Dans cet exemple, il y a 5 sujets tels que 4 d'entre eux n'ont que 0 valeur de résultat et l'un d'eux a des valeurs de résultat qui incluent une seule valeur de 1. (Chaque sujet a 3 valeurs de résultat au total.)
Même si vous donnez à chacun des sujets de cet exemple inventé une valeur de 1 pour leur première valeur de résultat, vous obtiendrez toujours le même message d'erreur lors de l'ajustement du modèle:
Outcome <- rep(0, 15)
Outcome[c(1, 4, 7, 10, 13)] <- 1
Cependant, si les 4 sujets qui n'avaient initialement que 0 valeurs sont autorisés à conserver ces valeurs et que le premier sujet reçoit deux valeurs de 1, le message d'erreur disparaît:
Outcome <- rep(0, 15)
Outcome[c(1,2)] <- 1
Une fois que vous comprenez mieux le modèle des valeurs 0 et 1 pour la variable de résultat parmi les sujets de votre étude, vous pouvez également essayer d'ajuster votre modèle avec la fonction mixed_model () du package GLMMadaptive dans R.
Pour le petit exemple fourni ici, cette fonction serait utilisée comme ceci:
library(GLMMadaptive)
m <- mixed_model(fixed = Outcome ~ 1,
random = ~ 1 | SubjectID,
data = Data,
family = binomial())
summary(m)
Isabella a fait d'excellents points. Cela peut également se produire lorsqu'il y a très peu de variation au niveau du canal. Les canaux sont peut-être très similaires les uns aux autres, de sorte que leur variance est vraiment proche de zéro et n'est donc pas nécessaire dans le modèle. Vous pouvez évaluer cela en ajustant un glm et voir si les inférences sont similaires.
Par curiosité, l'erreur survient-elle lorsque vous utilisez un estimateur alternatif? Il se peut que l'estimateur soit pour une raison quelconque bloqué sur une singularité. Vous pouvez juste essayer ce qui suit: mod.alt_est <- allFit(mod.detection_rand). Alternativement, vous pouvez avoir besoin d'une solution bayésienne pour aider à régulariser l'estimation et l'éloigner d'une singularité (essayez le blmepackage si la allFitfonction ne produit pas d'estimateur qui fonctionne).