混合效应 logistic 回归

Question

我正试图在Python中实现混合效应逻辑回归。作为比较的一点，我使用了R中lme4包中的glmer函数。
我发现statsmodels模块有一个BinomialBayesMixedGLM, 可以拟合这样的模型。但是，我遇到了一些问题：

1. 我认为statsmodels函数文档并不完全有用或明确，因此我不完全确定如何适当地使用该函数。
2. 到目前为止，我的尝试没有产生与我在R中使用glmer拟合模型时获得的结果相同的结果。
3. 由于它是贝叶斯的，我预计BinomialBayesMixedGLM函数不会计算p值，但我似乎无法弄清如何访问参数的完整后验分布。

作为测试案例，我正在使用此处提供的泰坦尼克号数据集。

import os
import pandas as pd
import statsmodels.genmod.bayes_mixed_glm as smgb

titanic = pd.read_csv(os.path.join(os.getcwd(), 'titanic.csv'))

r = {"Pclass": '0 + Pclass'}
mod = smgb.BinomialBayesMixedGLM.from_formula('Survived ~ Age', r, titanic)
fit = mod.fit_map()
fit.summary()

#           Type    Post. Mean  Post. SD       SD SD (LB) SD (UB)
# Intercept M           3.1623    0.3616            
# Age       M          -0.0380    0.0061            
# Pclass    V           0.0754    0.5669    1.078   0.347   3.351

然而，除了年龄对结果的影响，这个句子似乎与我在 R 中使用 glmer(Survived ~ Age + (1 | Pclass), data = titanic, family = "binomial") 得到的结果不符：

Random effects:
 Groups Name        Variance Std.Dev.
 Pclass (Intercept) 0.8563   0.9254  
Number of obs: 887, groups:  Pclass, 3

Fixed effects:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)  0.961780   0.573402   1.677   0.0935 .  
Age         -0.038708   0.006243  -6.200 5.65e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

我在用Python创建模型的时候出了什么错误？一旦解决了这个问题，我该如何提取后验概率或p值？最后，是否有任何与R中实现更相似的混合效应 logistic 回归的 Python 实现？

- YTD

---

@petezurich -- 为了我的知识，为什么要编辑标题，以至于它不再清楚地表明我的问题特别是关于如何在Python中拟合这种类型的模型？新标题似乎对于其他搜索类似问题答案的人来说描述性少了很多。 - YTD

这就是标签的用处。请查看此帖子，其中提到：“在问题标题中包含标签完全是不必要的。” - petezurich

@petezurich 感谢您的澄清和建议。经过更多的工作，我已经更新了我的问题。 - YTD

也许是小鹿班比或 Pymer4。 - Blue482

我是Bambi开发者。现在你可以使用贝叶斯统计学和Python来处理混合模型，就像在lme4中所做的那样。请查看https://github.com/bambinos/bambi，并随时打开问题以寻求帮助。 - Tomas Capretto

Answer 1

我刚刚也需要用到类似的Python代码，正如评论中建议的那样，Pymer4似乎提供了一个合适的方法（特别是如果你已经熟悉R）。使用问题中提到的示例数据集“titanic”：

from pymer4.models import Lmer

model = Lmer("Survived  ~ Age  + (1|Pclass)",
             data=titanic, family = 'binomial')

print(model.fit())

OUT:

Formula: Survived~Age+(1|Pclass)

Family: binomial     Inference: parametric

Number of observations: 887  Groups: {'Pclass': 3.0}

Log-likelihood: -525.812     AIC: 1057.624

Random effects:

               Name    Var    Std
Pclass  (Intercept)  0.856  0.925

No random effect correlations specified

Fixed effects:

             Estimate  2.5_ci  97.5_ci     SE     OR  OR_2.5_ci  OR_97.5_ci  \
(Intercept)     0.962  -0.162    2.086  0.573  2.616       0.85       8.050   
Age            -0.039  -0.051   -0.026  0.006  0.962       0.95       0.974   

              Prob  Prob_2.5_ci  Prob_97.5_ci  Z-stat  P-val  Sig  
(Intercept)  0.723        0.460         0.889   1.677  0.093    .  
Age          0.490        0.487         0.493  -6.200  0.000  *** 

只是额外的一点评论（抱歉偏离主要问题），我在一个 Ubuntu 20.04 机器上使用了 Python 3.8.8。不确定是否有其他人遇到过这个问题，但是在使用 Pymer4 运行上面的模型时，该程序包会抛出一个错误（当我尝试从 Pymer4 文档这里复制类似模型时也遇到同样的错误）：

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

这可以通过更新pymer4版本来解决（请参见这个相关问题）：

conda install -c ejolly -c conda-forge -c defaults pymer4=0.7.8

- Hal Koons

---

使用 design_matrix.any() 很可能是错误的，if design_matrix: 的预期含义是检查它是否为 None（替换为 if design_matrix is not None:）或者检查它不是一个空数组（if design_matrix.size > 0:）。 - Giacomo Petrillo