在 statsmodels 的混合线性模型中指定交叉效应的正确方法是什么？
Which is the correct way to specify crossed effects in a mixed linear model in statsmodels?
我对 statsmodels 中的交叉效应线性混合模型有疑问。具体来说，我可以看到两种表示数据的方式，但我不确定哪种方式合适。任何建议表示赞赏！
我的数据如下。我想确定一本书的 objective 质量（'good' 或 'bad'）是否可以预测这本书的得分。标题嵌套在质量变量中，但标题和评分者交叉。 （这是假数据，所以我不担心模型收敛。）
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-value 和系数，并说明了交叉效应。
但是，我在别处被告知像这样的交叉效应应该通过将数据集视为一组并完全使用方差分量来指定变化来指定。因此：
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-value、不同的系数和不同的模型测试统计数据。现在，要么我对这两种模型的使用都是错误的，要么我对其中一种模型的使用是错误的。谁能告诉我这是怎么回事？谢谢。
您正在比较两个完全不同的模型，这给出了不同的解释。
如果您对质量如何影响学校感兴趣，同时将其他协变量作为随机截距，则第一个模型是正确的。我很快检查并发现 提到了交叉效应，如果你 运行 模型在方差分量中包含所有随机截距，你会得到大致相同的结果：
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
这是'rater'是否出现在这个table中的问题。 VC 也用于指定更复杂的模型，例如可变斜率，所以如果它是纯随机截距，我认为你可以使用组和 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
与您的第一个模型或多或少相似，但由于您的数据，很难估计您的随机效应。
我对 statsmodels 中的交叉效应线性混合模型有疑问。具体来说，我可以看到两种表示数据的方式，但我不确定哪种方式合适。任何建议表示赞赏！
我的数据如下。我想确定一本书的 objective 质量（'good' 或 'bad'）是否可以预测这本书的得分。标题嵌套在质量变量中，但标题和评分者交叉。 （这是假数据，所以我不担心模型收敛。）
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-value 和系数，并说明了交叉效应。
但是，我在别处被告知像这样的交叉效应应该通过将数据集视为一组并完全使用方差分量来指定变化来指定。因此：
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-value、不同的系数和不同的模型测试统计数据。现在，要么我对这两种模型的使用都是错误的，要么我对其中一种模型的使用是错误的。谁能告诉我这是怎么回事？谢谢。
您正在比较两个完全不同的模型，这给出了不同的解释。
如果您对质量如何影响学校感兴趣，同时将其他协变量作为随机截距，则第一个模型是正确的。我很快检查并发现
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
这是'rater'是否出现在这个table中的问题。 VC 也用于指定更复杂的模型，例如可变斜率，所以如果它是纯随机截距，我认为你可以使用组和 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
与您的第一个模型或多或少相似，但由于您的数据，很难估计您的随机效应。