смешанные модели с двумя случайными эффектами - statsmodels

Question

import pandas as pd
import statsmodels.formula.api as smf

df = pd.read_csv('http://www.bodowinter.com/tutorial/politeness_data.csv')
df = df.drop(38)

В R Я бы сделал:

lmer(frequency ~ attitude + (1|subject) + (1|scenario), data=df)

который в R дает мне:

Random effects:
 Groups   Name        Variance Std.Dev.
 scenario (Intercept)  219     14.80   
 subject  (Intercept) 4015     63.36   
 Residual              646     25.42   
Fixed effects:
            Estimate Std. Error t value
(Intercept)  202.588     26.754   7.572
attitudepol  -19.695      5.585  -3.527

Я пытался сделать то же самое с statsmodels:

model = smf.mixedlm("frequency ~ attitude", data=df, groups=df[["subject","scenario"]]).fit()

Но model.summary() дает мне другой вывод:

      Mixed Linear Model Regression Results
=======================================================
Model:            MixedLM Dependent Variable: frequency
No. Observations: 83      Method:             REML     
No. Groups:       2       Scale:              0.0000   
Min. group size:  1       Likelihood:         inf      
Max. group size:  1       Converged:          Yes      
Mean group size:  1.0                                  
-------------------------------------------------------
                  Coef.  Std.Err. z P>|z| [0.025 0.975]
-------------------------------------------------------
Intercept        204.500                               
attitude[T.pol]    8.800                               
groups RE          0.000                               
=======================================================

Answer 1

Единственный способ, которым я мог бы придумать, как воспроизвести это, - просто объединить ваши группы.

df["grp"] = df["subject"].astype(str) + df["scenario"].astype(str)
model = smf.mixedlm("frequency ~ attitude", data=df, groups=df["grp"]).fit()

model.summary()
Out[87]: 
<class 'statsmodels.iolib.summary2.Summary'>
"""
            Mixed Linear Model Regression Results
==============================================================
Model:               MixedLM   Dependent Variable:   frequency
No. Observations:    83        Method:               REML     
No. Groups:          42        Scale:                615.6961 
Min. group size:     1         Likelihood:           -430.8261
Max. group size:     2         Converged:            Yes      
Mean group size:     2.0                                      
--------------------------------------------------------------
                 Coef.   Std.Err.   z    P>|z|  [0.025  0.975]
--------------------------------------------------------------
Intercept        202.588   10.078 20.102 0.000 182.836 222.340
attitude[T.pol]  -19.618    5.476 -3.582 0.000 -30.350  -8.885
groups RE       3650.021   50.224                             
==============================================================

"""

Answer 2

Код ниже воспроизводит результаты R. Поскольку это скрещенная модель без независимых групп, вам нужно поместить всех в одну группу и указать случайные эффекты, используя компоненты дисперсии.

import pandas as pd                                                                                                        
import statsmodels.api as sm                                                                                               

df = pd.read_csv('http://www.bodowinter.com/tutorial/politeness_data.csv')                                                 
df = df.dropna()                                                                                                           
df["group"] = 1                                                                                                            

vcf = {"scenario": "0 + C(scenario)", "subject": "0 + C(subject)"}                                                         
model = sm.MixedLM.from_formula("frequency ~ attitude", groups="group",                                                    
                                vc_formula=vcf, re_formula="0", data=df)                                                   
result = model.fit()  

Вот результаты:

            Mixed Linear Model Regression Results
==============================================================
Model:               MixedLM   Dependent Variable:   frequency
No. Observations:    83        Method:               REML     
No. Groups:          1         Scale:                646.0163 
Min. group size:     83        Likelihood:           -396.7268
Max. group size:     83        Converged:            Yes      
Mean group size:     83.0                                     
--------------------------------------------------------------
                 Coef.   Std.Err.   z    P>|z|  [0.025  0.975]
--------------------------------------------------------------
Intercept        202.588   26.754  7.572 0.000 150.152 255.025
attitude[T.pol]  -19.695    5.585 -3.526 0.000 -30.641  -8.748
scenario Var     218.991    6.476                             
subject Var     4014.616  104.614                             
==============================================================

Answer 3

lmer эквивалент вашей smf.mixedlm модели будет выглядеть примерно так:

lmer(frequency ~ attitude + (1 + attitude|subject) + (1 + attitude|scenario), data = df)

Объяснение терминов:

1. Глобальный перехват (вы можете отключить глобальный перехват с помощью frequency ~ 0 + attitude + ...)
2. Глобальный наклон для фиксированного эффекта attitude.
3. Случайный перехват vor subject (то есть для каждого уровня subject вы получаете отклонение от глобального перехвата), а отклонение от фиксированного наклона эффекта для attitude в пределах каждого уровня subject, что позволяет для корреляции между случайным перехватом и наклоном.
4. Эквивалентные термины случайного перехвата и наклона для scenario.

Обратите внимание: если вы хотите, чтобы произвольный перехват и наклон могли свободно изменяться (т. Е. Для обеспечения нулевой корреляции между перехватом и наклоном), вам придется заменить (1 + attitude|subject) на (1|subject) + (0 + attitude|subject), и аналогично для scenario.