python - 混合效应逻辑回归

Question

我正在尝试在 python 中实现混合效果逻辑回归。作为比较点，我使用的是 glmer来自 lme4 的函数R中的包。

我发现 statsmodels模块有一个 BinomialBayesMixedGLM那应该能够适合这样的模型。但是，我遇到了一些问题:

我找到了 statsmodels 的文档功能不完全有用或清晰，所以我不完全确定如何正确使用该功能。

到目前为止，我的尝试还没有产生可以复制我在使用 glmer 拟合模型时得到的结果。在河

我期待 BinomialBayesMixedGLM函数不计算 p 值，因为它是贝叶斯，但我似乎无法弄清楚如何访问参数的完整后验分布。

作为测试用例，我使用了可用的 Titanic 数据集 here .

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 值？最后，他们的混合效应逻辑回归的任何 python 实现是否更类似于 R 中的实现？

Answer 1

只需要对 Python 做类似的事情，正如评论中所建议的 Pymer4似乎提供了一种合适的方法(尤其是如果您熟悉 R 无论如何)。
使用问题中提到的示例数据集“titanic”:

from pymer4.models import Lmer

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

print(model.fit())

出去:

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 文档 here 复制类似模型时出现相同错误):

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

引用/pymer4/models/Lmer.py 包文件中的这一行:

--> 444         if design_matrix:

我通过将其更改为来解决此问题(不确定这是否是最优雅或最安全的方法，很高兴在这里得到纠正):

if design_matrix.any():

在我测试的少数情况下，这似乎使程序包运行并提供了与 R 等效的结果。
以下评论中建议的更好方法可能是:

if design_matrix is not None:

感谢 Giacomo Petrillo 指出这一点，但是我没有对此进行测试。