如何在Python statsmodels线性混合效应模型中拥有多个组？

Question

我正在尝试使用 Python statsmodels 线性混合效应模型来拟合具有两个随机截距的模型，例如两组。我不知道如何初始化模型以便我可以做到这一点。

这是例子。我的数据如下所示（取自here）：

subject gender  scenario    attitude    frequency
F1  F   1   pol 213.3
F1  F   1   inf 204.5
F1  F   2   pol 285.1
F1  F   2   inf 259.7
F1  F   3   pol 203.9
F1  F   3   inf 286.9
F1  F   4   pol 250.8
F1  F   4   inf 276.8

我想制作一个具有两种随机效应的线性混合效应模型——一种用于主题组，一种用于场景组。我正在尝试这样做：

import statsmodels.api as sm
model = sm.MixedLM.from_formula("frequency ~ attitude + gender", data, groups=data[['subject', 'scenario']])
result = model.fit()
print result.summary()

我不断收到此错误：

LinAlgError: Singular matrix

它在 R 中工作得很好。当我在 R 中使用

lme4

进行基于公式的渲染时，它非常适合：

politeness.model = lmer(frequency ~ attitude + gender + 
        (1|subject)  + (1|scenario), data=politeness)

我不明白为什么会发生这种情况。当我使用任何一个随机效果/组时，它会起作用，例如

model = sm.MixedLM.from_formula("frequency ~ attitude + gender", data, groups=data['subject'])

然后我得到：

                 Mixed Linear Model Regression Results
===============================================================
Model:                MixedLM   Dependent Variable:   frequency
No. Observations:     83        Method:               REML     
No. Groups:           6         Scale:                850.9456 
Min. group size:      13        Likelihood:           -393.3720
Max. group size:      14        Converged:            Yes      
Mean group size:      13.8                                     
---------------------------------------------------------------
                 Coef.   Std.Err.   z    P>|z|  [0.025   0.975]
---------------------------------------------------------------
Intercept        256.785   15.226 16.864 0.000  226.942 286.629
attitude[T.pol]  -19.415    6.407 -3.030 0.002  -31.972  -6.858
gender[T.M]     -108.325   21.064 -5.143 0.000 -149.610 -67.041
Intercept RE     603.948   23.995                              
===============================================================

或者，如果我这样做：

model = sm.MixedLM.from_formula("frequency ~ attitude + gender", data, groups=data['scenario'])

这是我得到的结果：

              Mixed Linear Model Regression Results
================================================================
Model:               MixedLM    Dependent Variable:    frequency
No. Observations:    83         Method:                REML     
No. Groups:          7          Scale:                 1110.3788
Min. group size:     11         Likelihood:            -402.5003
Max. group size:     12         Converged:             Yes      
Mean group size:     11.9                                       
----------------------------------------------------------------
                 Coef.   Std.Err.    z    P>|z|  [0.025   0.975]
----------------------------------------------------------------
Intercept        256.892    8.120  31.637 0.000  240.977 272.807
attitude[T.pol]  -19.807    7.319  -2.706 0.007  -34.153  -5.462
gender[T.M]     -108.603    7.319 -14.838 0.000 -122.948 -94.257
Intercept RE     182.718    5.502                               
================================================================

我不知道发生了什么事。我觉得我在问题的统计数据中遗漏了一些基本的东西。

Answer 1

您正在尝试拟合具有“交叉随机效应”的模型，即，您希望允许跨场景的受试者之间存在一致的变化，以及跨场景的受试者之间存在一致的变化。您可以在统计模型中使用多个随机效应项，但它们必须是嵌套的。拟合交叉（而不是嵌套）随机效应需要更复杂的算法，事实上 statsmodels 文档 说（截至 2016 年 8 月 25 日，添加了强调）：

当前实现的一些限制是，它不支持残差误差上更复杂的结构（它们始终是同方差的），并且

不支持交叉随机效应

。我们希望在下一个版本中实现这些功能。

据我所知，您的选择是（1）退回到嵌套模型（即拟合模型，就像任一场景嵌套在主题中一样，或者

反之亦然

- 或者尝试两者并看看差异是否重要）； (2) 在 R 内或通过 rpy2 回退到

lme4

。 一如既往，您有权获得使用 statsmodels 所支付费用的全额退款......

Answer 2

我没有找到示例，并且文档似乎只更新了部分。

单元测试包含使用 MixedLM 公式接口的示例：

https://github.com/statsmodels/statsmodels/blob/master/statsmodels/regression/tests/test_lme.py#L284

Answer 3

https://eshinjolly.com/pymer4/

这是我使用它的代码示例：

from pymer4.models import Lmer
selected_columns = ['subid', 'accuracy', 'rt', 'response_switch', 'classwitch', 'Experiment']
df = df[selected_columns]
model = Lmer('accuracy ~ classwitch * response_switch + (1|subid) + (1|Experiment)', data=df, family='binomial')
result = model.fit()
result