小编YTD*_*YTD的帖子

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

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

1. 我发现该statsmodels函数的文档并不完全有用或清晰，所以我不完全确定如何正确使用该函数。
2. 到目前为止，我的尝试还没有产生可以复制我glmer在 R 中拟合模型时得到的结果。
3. 我希望该BinomialBayesMixedGLM函数不计算 p 值，因为它是贝叶斯，但我似乎无法弄清楚如何访问参数的完整后验分布。

作为测试用例，我使用了可用的泰坦尼克号数据集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 …

我写了一个脚本，相信应该在Python和R中产生相同的结果，但是它们产生的答案却截然不同。每种方法都尝试通过使用Nelder-Mead使偏差最小化来使模型适合模拟数据。总体而言，R的乐观表现要好得多。难道我做错了什么？R和SciPy中实现的算法是否不同？

Python结果：

>>> res = minimize(choiceProbDev, sparams, (stim, dflt, dat, N), method='Nelder-Mead')

 final_simplex: (array([[-0.21483287, -1.        , -0.4645897 , -4.65108495],
       [-0.21483909, -1.        , -0.4645915 , -4.65114839],
       [-0.21485426, -1.        , -0.46457789, -4.65107337],
       [-0.21483727, -1.        , -0.46459331, -4.65115965],
       [-0.21484398, -1.        , -0.46457725, -4.65099805]]), array([107.46037865, 107.46037868, 107.4603787 , 107.46037875,
       107.46037875]))
           fun: 107.4603786452194
       message: 'Optimization terminated successfully.'
          nfev: 349
           nit: 197
        status: 0
       success: True
             x: array([-0.21483287, -1.        , -0.4645897 , -4.65108495])

R结果：

> res <- optim(sparams, choiceProbDev, stim=stim, dflt=dflt, dat=dat, N=N,
             method="Nelder-Mead")

$par
[1] …