我正在尝试在 python 中实现混合效应逻辑回归。作为比较,我正在使用glmer
函数从lme4
R 中的包。
我发现statsmodels
模块有一个BinomialBayesMixedGLM
应该能够适合这样的模型。然而,我遇到了很多问题:
- 我找到了有关的文档
statsmodels
函数并不完全有帮助或不清楚,所以我不完全确定如何正确使用该函数。
- 到目前为止,我的尝试还没有产生复制我在拟合模型时得到的结果
glmer
in R.
- 我期望
BinomialBayesMixedGLM
函数不计算 p 值,因为它是贝叶斯函数,但我似乎无法弄清楚如何访问参数的完整后验分布。
作为测试用例,我使用可用的泰坦尼克号数据集here https://web.stanford.edu/class/archive/cs/cs109/cs109.1166/problem12.html.
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
然而,除了 Age 的斜率之外,这似乎与我在 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 中的实现?