Binomial glmer () singular apesar de “muitos” dados
Tenho lido diferentes perguntas sobre como é fácil esbarrar em singularidades ao ajustar modelos de efeitos mistos glmer(). Em geral, a ideia é que singularidades podem surgir de estruturas aleatórias muito complexas. Se a estrutura aleatória for simples, também pode acontecer quando os dados não são suficientes para calcular a matriz de variância-covariância ... veja por exemplo esta página de Ben Bolker , a resposta de Robert Long a este post ou a página de ajuda deisSingular() .
No entanto, o modelo que estou tentando ajustar é muito simples:
mod.detection_rand <- glmer(reaction ~ Pedra + (1|Channel), family="binomial", data = garotes)
boundary (singular) fit: see ?isSingular
... e, aparentemente, tenho dados suficientes para as diferentes combinações de variáveis preditoras (fixas e aleatórias):
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
O que você acha?
EDIT: Aqui está o resumo do 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: Seguindo o comentário 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: Seguindo a resposta de Isabella:
Verifiquei a estrutura da variável de resultado ( reaction). Aqui está a tabela 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, existem os dois tipos de resultados para todos Channelse todos os Pedratratamentos ... então não é como o exemplo que a Isabella apresentou ... além disso, tentei modelar este GLMM com o library(GLMMadaptive)e também não convergiu.
EDIT4: O conjunto de dados que estou usando, caso alguém esteja curioso.
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
Muito obrigado por todas as suas respostas, em qualquer caso! Aprendendo muito com eles!
Respostas
Um comentário adicional: eu dei uma olhada em seus dados e está claro, novamente, que não há evidência de variação sistemática entre os diferentes canais. É por isso que o modelo misto estima que a variância entre os canais seja$0$, tornando o modelo singular.
Você pode ver isso na figura abaixo, onde os erros padrão para quase todos os canais se sobrepõem ...
... e pode confirmar por ANOVA decomposição de um GLM de efeitos fixos, mostrando que não há efeito principal significativo de Canal (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
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
Como este é um modelo de regressão logística binária de efeitos mistos, ele assume que sua variável de resultado é binária com valores codificados como 0 ou 1.
O que você precisa investigar é se você tem 1s suficientes presentes em sua variável de resposta para um número suficiente de 'sujeitos'. (No seu caso, assunto significa canal.)
Aqui está um exemplo inventado que produz o mesmo aviso que você recebeu:
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)
Neste exemplo, existem 5 assuntos, de modo que 4 deles têm apenas 0 valores de resultado e um deles tem valores de resultado que incluem um único valor de 1. (Cada assunto tem 3 valores de resultado no total).
Mesmo se você der a cada um dos sujeitos neste exemplo inventado um valor de 1 para seu primeiro valor de resultado, você ainda receberá a mesma mensagem de erro ao ajustar o modelo:
Outcome <- rep(0, 15)
Outcome[c(1, 4, 7, 10, 13)] <- 1
No entanto, se todos os 4 sujeitos que inicialmente tinham apenas 0 valores puderem manter esses valores e o primeiro sujeito receber dois valores de 1, a mensagem de erro desaparecerá:
Outcome <- rep(0, 15)
Outcome[c(1,2)] <- 1
Depois de entender melhor o padrão de valores 0 e 1 para a variável de resultado entre os sujeitos de seu estudo, a outra coisa que você pode tentar é ajustar seu modelo com a função mixed_model () do pacote GLMMadaptive em R.
Para o pequeno exemplo fornecido aqui, esta função seria usada assim:
library(GLMMadaptive)
m <- mixed_model(fixed = Outcome ~ 1,
random = ~ 1 | SubjectID,
data = Data,
family = binomial())
summary(m)
Isabella fez alguns pontos excelentes. Isso também pode acontecer quando há muito pouca variação no nível do canal. Talvez os canais sejam muito semelhantes uns aos outros, então sua variância é realmente próxima de zero e, portanto, não é necessária no modelo. Você pode avaliar isso ajustando um glm e ver se as inferências são semelhantes.
Por curiosidade, o erro surge quando você usa um estimador alternativo? Pode ser que o estimador esteja, por algum motivo, travando em uma singularidade. Você pode apenas tentar o seguinte: mod.alt_est <- allFit(mod.detection_rand). Alternativamente, você pode precisar de uma solução bayesiana para ajudar a regularizar a estimativa e afastá-la de uma singularidade (tente blmepackage se a allFitfunção não produzir um estimador que funcione).