python - 在 statsmodels 的混合线性模型中指定交叉效应的正确方法是什么？

Question

我有一个关于 statsmodels 中的交叉效应线性混合模型的问题。具体来说，我可以看到两种表示我的数据的方式，但我不确定哪种方式合适。任何建议表示赞赏!

我的数据如下。我想确定一本书的客观质量(“好”或“差”)是否能预测这本书的得分。头衔嵌套在质量变量中，但头衔和评分者是交叉的。 (这是假数据，所以我不担心模型会收敛。)

rater   title   quality score
john    book_1  good    0.600833333
frank   book_2  bad     0.683020833
emma    book_3  good    0.653645833
john    book_4  bad     0.6528125
frank   book_5  good    0.6040625
emma    book_1  good    0.600833333
john    book_2  bad     0.522
frank   book_3  good    0.600833333
emma    book_4  bad     0.619464286
john    book_5  good    0.600833333
frank   book_1  good    0.57125
emma    book_2  bad     0.6296875
john    book_3  good    0.607205882
frank   book_4  bad     0.61203125
emma    book_5  good    0.600833333

分析这些数据的一种方法是将质量作为我的自变量，将分数作为我的因变量，将评分者作为我的分组变量，并使用方差分量来捕捉对标题的交叉影响。这给出:

import statsmodels.api as sm
import statsmodels.fomula.api as smf

md = smf.mixedlm('score ~ quality', vc_formula = {"title":"0 + title"}, groups = data['rater'], data = data).fit().summary()

型号概要:

           Mixed Linear Model Regression Results
===========================================================
Model:              MixedLM   Dependent Variable:   score  
No. Observations:   15        Method:               REML   
No. Groups:         3         Scale:                0.0007 
Min. group size:    5         Log-Likelihood:       22.1997
Max. group size:    5         Converged:            Yes    
Mean group size:    5.0                                    
-----------------------------------------------------------
                Coef.  Std.Err.    z    P>|z| [0.025 0.975]
-----------------------------------------------------------
Intercept        0.620    0.001 841.098 0.000  0.618  0.621
quality[T.good] -0.015    0.013  -1.158 0.247 -0.041  0.011
title Var        0.001                                     
===========================================================

在我看来，这直觉上是正确的方法。它为我的 IV 提供了 p 值和系数，并说明了交叉效应。

但是，我在其他地方被告知，应该通过将数据集视为一组并完全使用方差分量来指定变化来指定这样的交叉效应。因此:

data['groups'] = 1

md = smf.mixedlm('score ~ 1', vc_formula = {"rater":"0 + rater", "title":"0 + title", "quality":"0 + quality"}, groups = data['groups'], data = data).fit().summary()

产量:

        Mixed Linear Model Regression Results
=====================================================
Model:            MixedLM Dependent Variable: score  
No. Observations: 15      Method:             REML   
No. Groups:       1       Scale:              0.0013 
Min. group size:  15      Log-Likelihood:     24.4023
Max. group size:  15      Converged:          No     
Mean group size:  15.0                               
-----------------------------------------------------
            Coef. Std.Err.   z    P>|z| [0.025 0.975]
-----------------------------------------------------
Intercept   0.612    0.017 35.011 0.000  0.577  0.646
quality Var 0.000                                    
rater Var   0.000    0.020                           
title Var   0.000                                    
=====================================================

这个模型没有提供 p 值、不同的系数和不同的模型测试统计数据。现在，我要么错误地使用了这两种模型，要么错误地使用了其中之一。谁能告诉我这是什么情况？谢谢。

Answer 1

您正在比较两个完全不同的模型，这给出了不同的解释。

如果您对质量如何影响学校感兴趣，同时将其他协变量作为随机截距，则第一个模型是正确的。我赶紧查了一下，发现this post提到交叉效应，如果您在方差分量中使用所有随机截距运行模型，您会得到大致相同的结果:

data['group'] = 1
md2 = smf.mixedlm('score ~ quality', 
                  vc_formula = {"title":"0 + title","rater":"0 + rater"}, 
                  groups = data['group'], data = data).fit().summary()

Model:  MixedLM Dependent Variable: score
No. Observations:   15  Method: REML
No. Groups: 1   Scale:  0.0014
Min. group size:    15  Log-Likelihood: 22.0951
Max. group size:    15  Converged:  No
Mean group size:    15.0        
Coef.   Std.Err.    z   P>|z|   [0.025  0.975]
Intercept   0.620   0.016   38.313  0.000   0.588   0.652
quality[T.good] -0.015  0.021   -0.736  0.462   -0.056  0.026
rater Var   0.000   0.030               
title Var   0.000           

问题在于“评估者”是否出现在此表中。 VC 也用于指定更复杂的模型，例如可变斜率，因此如果它是纯随机截距，我认为您可以使用 group 和 VC 的组合。

在第二个模型中，您只使用截距对分数进行建模，与质量无关，如果您对质量感兴趣，这没有意义。

最后，我们可以使用以下命令检查 R 中的结果:

df = structure(list(rater = structure(c(3L, 2L, 1L, 3L, 2L, 1L, 3L, 
2L, 1L, 3L, 2L, 1L, 3L, 2L, 1L), .Label = c("emma", "frank", 
"john"), class = "factor"), title = structure(c(1L, 2L, 3L, 4L, 
5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L), .Label = c("book_1", 
"book_2", "book_3", "book_4", "book_5"), class = "factor"), quality = structure(c(2L, 
1L, 2L, 1L, 2L, 2L, 1L, 2L, 1L, 2L, 2L, 1L, 2L, 1L, 2L), .Label = c("bad", 
"good"), class = "factor"), score = c(0.600833333, 0.683020833, 
0.653645833, 0.6528125, 0.6040625, 0.600833333, 0.522, 0.600833333, 
0.619464286, 0.600833333, 0.57125, 0.6296875, 0.607205882, 0.61203125, 
0.600833333)), class = "data.frame", row.names = c(NA, -15L))

library(lme4)
summary(lmer(score ~ quality + (1|rater) + (1|title),data=df))
boundary (singular) fit: see ?isSingular
Linear mixed model fit by REML ['lmerMod']
Formula: score ~ quality + (1 | rater) + (1 | title)
   Data: df

REML criterion at convergence: -44.4

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.60015 -0.09695 -0.09695  0.16712  1.67924 

Random effects:
 Groups   Name        Variance Std.Dev.
 title    (Intercept) 0.000000 0.00000 
 rater    (Intercept) 0.000000 0.00000 
 Residual             0.001416 0.03763 
Number of obs: 15, groups:  title, 5; rater, 3

Fixed effects:
            Estimate Std. Error t value
(Intercept)  0.61984    0.01536  40.351
qualitygood -0.01535    0.01983  -0.774

Correlation of Fixed Effects:
            (Intr)
qualitygood -0.775
convergence code: 0
boundary (singular) fit: see ?isSingular

或多或少与您的第一个模型相似，但由于您的数据，很难估计您的随机效应。