Binomial glmer () Singular trotz "vieler" Daten
Ich habe verschiedene Fragen darüber gelesen, wie einfach es ist, beim Anpassen von Modellen mit gemischten Effekten auf Singularitäten zu stoßen glmer(). Im Allgemeinen besteht die Idee darin, dass Singularitäten aus sehr komplexen Zufallsstrukturen entstehen können. Wenn die Zufallsstruktur einfach ist, kann es auch vorkommen, dass die Daten nicht ausreichen, um die Varianz-Kovarianz-Matrix zu berechnen. Siehe zum Beispiel diese Seite von Ben Bolker , die Antwort von Robert Long auf diesen Beitrag oder die Hilfeseite vonisSingular() .
Das Modell, das ich anpassen möchte, ist jedoch sehr einfach:
mod.detection_rand <- glmer(reaction ~ Pedra + (1|Channel), family="binomial", data = garotes)
boundary (singular) fit: see ?isSingular
... und anscheinend habe ich genug Daten für die verschiedenen (festen und zufälligen) Prädiktorvariablenkombinationen:
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
Was denken Sie?
EDIT: Hier ist die Zusammenfassung des Modells, 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: Nach Billys Kommentar:
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: Nach Isabellas Antwort:
Ich habe die Struktur der Ergebnisvariablen ( reaction) überprüft . Hier ist die Ergebnistabelle:
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
Anscheinend gibt es beide Arten von Ergebnissen für alle Channelsund alle PedraBehandlungen ... es ist also nicht wie das Beispiel, das Isabella vorgestellt hat ... außerdem habe ich versucht, dieses GLMM mit dem zu modellieren, library(GLMMadaptive)und es konvergierte auch nicht.
EDIT4: Der Datensatz, den ich verwende, falls jemand neugierig ist.
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
Vielen Dank für all Ihre Antworten auf jeden Fall! Viel von ihnen lernen!
Antworten
Ein weiterer Kommentar: Ich habe mir Ihre Daten angesehen und es ist wieder klar, dass es keine Hinweise auf systematische Abweichungen zwischen den verschiedenen Kanälen gibt. Aus diesem Grund schätzt das gemischte Modell die Varianz zwischen den Kanälen$0$und macht das Modell einzigartig.
Sie können dies in der folgenden Abbildung sehen, in der sich die Standardfehler für fast jeden Kanal überlappen ...
... und kann dies durch ANOVA-Zerlegung eines GLM mit festen Effekten bestätigen, was zeigt, dass es keinen signifikanten Haupteffekt von Channel gibt (p = 0,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
Da es sich um ein binäres logistisches Regressionsmodell mit gemischten Effekten handelt, wird davon ausgegangen, dass Ihre Ergebnisvariable binär ist und die Werte entweder als 0 oder 1 codiert sind.
Sie müssen untersuchen, ob in Ihrer Antwortvariablen genügend Einsen für eine ausreichende Anzahl von "Probanden" vorhanden sind. (In Ihrem Fall steht Betreff für Kanal.)
Hier ist ein erfundenes Beispiel, das dieselbe Warnung ausgibt wie das, was Sie erhalten haben:
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)
In diesem Beispiel gibt es 5 Probanden, so dass 4 von ihnen nur 0 Ergebniswerte haben und einer von ihnen Ergebniswerte hat, die einen einzelnen Wert von 1 enthalten. (Jeder Proband hat insgesamt 3 Ergebniswerte.)
Selbst wenn Sie jedem Thema in diesem zusammengesetzten Beispiel den Wert 1 für seinen ersten Ergebniswert geben, wird beim Anpassen des Modells immer noch dieselbe Fehlermeldung angezeigt:
Outcome <- rep(0, 15)
Outcome[c(1, 4, 7, 10, 13)] <- 1
Wenn jedoch alle 4 Probanden, die ursprünglich nur 0 Werte hatten, diese Werte beibehalten dürfen und der erste Proband zwei Werte von 1 erhält, verschwindet die Fehlermeldung:
Outcome <- rep(0, 15)
Outcome[c(1,2)] <- 1
Sobald Sie das Muster der 0- und 1-Werte für die Ergebnisvariable unter Ihren Studienteilnehmern besser verstanden haben, können Sie versuchen, Ihr Modell mit der Funktion mix_model () aus dem GLMMadaptive- Paket in R anzupassen .
Für das hier bereitgestellte kleine Beispiel würde diese Funktion folgendermaßen verwendet:
library(GLMMadaptive)
m <- mixed_model(fixed = Outcome ~ 1,
random = ~ 1 | SubjectID,
data = Data,
family = binomial())
summary(m)
Isabella machte einige hervorragende Punkte. Dies kann auch passieren, wenn auf Kanalebene nur sehr geringe Abweichungen auftreten. Möglicherweise sind die Kanäle einander nur sehr ähnlich, sodass ihre Varianz tatsächlich nahe Null liegt und daher im Modell nicht benötigt wird. Sie können dies beurteilen, indem Sie ein glm anpassen und prüfen, ob die Schlussfolgerungen ähnlich sind.
Tritt der Fehler aus Neugier auf, wenn Sie einen alternativen Schätzer verwenden? Es könnte sein, dass der Schätzer aus irgendeinem Grund an einer Singularität hängen bleibt. Sie können einfach Folgendes versuchen : mod.alt_est <- allFit(mod.detection_rand). Alternativ benötigen Sie möglicherweise eine Bayes'sche Lösung, um die Schätzung zu regulieren und von einer Singularität zu entfernen (versuchen Sie das blmePaket, wenn die allFitFunktion keinen funktionierenden Schätzer erzeugt).