R语言 如何在Python statmodels线性混合效应模型中有多个组?

mtb9vblg  于 11个月前  发布在  Python
关注(0)|答案(4)|浏览(139)

我试图使用Python statmodels线性混合效应模型来拟合一个有两个随机截距的模型,例如两个组。我不知道如何初始化模型,这样我就可以做到这一点。
下面是一个例子。我有一个看起来像下面这样的数据(取自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                               
================================================================

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

rqqzpn5f

rqqzpn5f1#

您试图拟合具有 * 交叉随机效应 * 的模型,即您希望允许受试者在不同情景之间的一致变异,以及受试者在不同情景之间的一致变异。您可以 * 在统计模型中使用多个随机效应项,但它们必须嵌套。拟合交叉(与嵌套相反)随机效应需要更复杂的算法,实际上statsmodels documentation说(截至2016年8月25日,重点添加):
当前实现的一些限制是它不支持残差上更复杂的结构(它们总是同方差的),并且它不支持交叉随机效应。我们希望在下一个版本中实现这些功能。
据我所知,你的选择是(1)回退到嵌套模型(即,拟合模型,好像任何一个场景都嵌套在主题中,或者 * 反之亦然 * -或者尝试两者,看看差异是否重要);(2)回退到lme4,无论是在R中还是通过rpy2
像往常一样,你有权获得全额退款的钱,你支付使用statsmodels.

dbf7pr2w

dbf7pr2w2#

多重或交叉随机截距交叉效应可以使用方差分量拟合,其实现方式与单组混合效应不同。
我没有找到一个例子,文档似乎只是部分更新。
单元测试包含一个使用MixedLM公式接口的示例:
https://github.com/statsmodels/statsmodels/blob/master/statsmodels/regression/tests/test_lme.py#L284

jum4pzuy

jum4pzuy3#

您可以使用Pymer4在Python中对多个组执行此操作。
下面是一个使用它的代码示例:

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

字符串

nhjlsmyf

nhjlsmyf4#

您可以在statsmodel Mixedlm中使用下面的方法(vc_formula),得到与R代码类似的结果(politeness.model = lmer(frequency ~ attitude + gender +(1|主题)+(1| scenario),data=politeness))

import pandas as pd
import numpy as np
import statsmodels.api as sm
import random
import statsmodels.api as sm
import statsmodels.formula.api as smf
df = pd.DataFrame(politeness)
fixed_formula='frequency ~ attitude + gender'
vc_formula={'g1': "0 + C(subject)",
'g2': "0 + C(scenario)"
 }
oo=np.ones(df.shape[0])
model=sm.MixedLM.from_formula(fixed_formula, data=df, groups=oo, 
 vc_formula=vc_formula)
results=model.fit()
print(results.summary())
random_effects_by_channel = results.random_effects
print(random_effects_by_channel)

字符串

相关问题