机器算法验证 - 尽管有“大量”数据，但二项式 glmer() 单数 - 吾爱随笔录

尽管有“大量”数据，但二项式 glmer() 单数

机器算法验证 r 混合模式广义线性模型 lme4-nlme 随机效应模型

2022-03-24 09:56:37

我一直在阅读不同的问题，这些问题是关于在使用glmer(). 一般来说，这个想法是奇点可能来自非常复杂的随机结构。如果随机结构很简单，那么当数据不足以计算方差 - 协方差矩阵时也可能发生这种情况......例如参见Ben Bolker 的这个页面、Robert Long 对这篇文章的回答或的帮助页面isSingular()。

但是，我要拟合的模型非常简单：

mod.detection_rand <- glmer(reaction ~ Pedra + (1|Channel), family="binomial", data = garotes)
boundary (singular) fit: see ?isSingular

...显然我有足够的数据用于不同的（固定和随机）预测变量组合：

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

你怎么看？

编辑：这是模型的摘要，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：按照比利的评论：

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：按照伊莎贝拉的回答：

我检查了结果变量 ( reaction) 的结构。这是结果表：

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

Channels显然，所有和所有治疗都有两种类型的结果Pedra......所以它不像伊莎贝拉提出的例子......此外，我试图用来模拟这个 GLMM，但library(GLMMadaptive)它也没有收敛。

EDIT4：我正在使用的数据集，以防有人好奇。

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

无论如何，非常感谢您的所有回复！向他们学习很多！

4个回答

伊莎贝拉提出了一些很好的观点。当通道级别的变化很小时，也会发生这种情况。也许通道彼此非常相似，因此它们的方差实际上接近于零，因此模型中不需要。您可以通过拟合 glm 来评估这一点，看看推论是否相似。

因为这是一个混合效应二元逻辑回归模型，它假设您的结果变量是二元的，其值编码为 0 或 1。

您需要调查的是您的响应变量中是否有足够的 1 用于足够数量的“受试者”。（在您的情况下，主题代表频道。）

这是一个虚构的示例，它会产生与您得到的 wnat 相同的警告：

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)

在此示例中，有 5 个受试者，其中 4 个只有 0 个结果值，其中一个具有包含单个值 1 的结果值。（每个受试者总共有 3 个结果值。）

即使你在这个虚构的例子中给每个主题的第一个结果值赋值为 1，在拟合模型时你仍然会得到相同的错误消息：

Outcome <- rep(0, 15)

Outcome[c(1, 4, 7, 10, 13)] <- 1

但是，如果最初只有 0 个值的所有 4 个主题都被允许保留这些值并且第一个主题收到两个值 1，则错误消息将消失：

Outcome <- rep(0, 15)

Outcome[c(1,2)] <- 1

一旦您更好地理解了研究对象中结果变量的 0 和 1 值的模式，您可以尝试的另一件事是使用R 中GLMMadaptive包中的mixed_model()函数来拟合您的模型。

对于此处提供的小示例，此函数将像这样使用：

library(GLMMadaptive)

m <- mixed_model(fixed = Outcome ~ 1, 
                 random = ~ 1 | SubjectID, 
                 data = Data,
                 family = binomial())
summary(m)

进一步评论：我查看了您的数据，很明显，没有证据表明不同渠道之间存在系统差异。这就是为什么混合模型估计通道间方差为 $0$ ，使模型奇异。

您可以在下图中看到这一点，其中几乎每个通道的标准误差都重叠......

...并且可以通过固定效应 GLM 的 ANOVA 分解来确认它，表明通道没有显着的主效应 (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

代码

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

出于好奇，当您使用替代估计器时会出现错误吗？估计器可能由于某种原因卡在奇点上。您可以尝试以下方法：mod.alt_est <- allFit(mod.detection_rand). 或者，您可能需要一个贝叶斯解决方案来帮助规范化估计并使其远离奇点（如果该函数不能产生有效的估计器，请尝试blme包）。allFit

其它你可能感兴趣的问题

上一篇回归是否适用于非正态分布的数据？下一篇机器学习VS统计学习VS统计