在Python的Statsmodels中实现交叉随机效应模型的疑问

一、R端数据生成与模型拟合 #

数据生成与模型代码 #

set.seed(123)

# 被试和项目数量
num_subjects <- 30
num_items <- 20

# 生成被试和项目ID
subject <- factor(rep(1:num_subjects, each=num_items))
item <- factor(rep(1:num_items, num_subjects))

# 生成固定效应预测变量
x <- rnorm(num_subjects * num_items, mean=0, sd=1)

# 生成被试和项目的随机截距
subject_intercepts <- rnorm(num_subjects, mean=0, sd=2)
item_intercepts <- rnorm(num_items, mean=0, sd=1)

# 创建数据框
dt <- data.frame(subject, item, x)

# 添加随机截距到数据框
dt$subject_intercept <- subject_intercepts[dt$subject]
dt$item_intercept <- item_intercepts[dt$item]

# 生成残差
residual_error <- rnorm(num_subjects * num_items, mean=0, sd=1)

# 生成响应变量y
dt$y <- 5 + 3 * dt$x + dt$subject_intercept + dt$item_intercept + residual_error

# 拟合线性混合效应模型
model <- lmer(y ~ x + (1|subject) + (1|item), data=dt)

# 模型摘要
summary(model)

# 保存数据（移除随机截距列）
dt$subject_intercept <- NULL
dt$item_intercept <- NULL
write.csv(dt, "crossed_re_01.csv")

R模型方差成分输出 #

Random effects:
 Groups   Name        Variance Std.Dev.
 subject  (Intercept) 4.818    2.195   
 item     (Intercept) 0.778    0.882   
 Residual             1.030    1.015   

二、Python端模型拟合与结果 #

Python拟合代码 #

import numpy as np
import pandas as pd
import statsmodels.api as sm
from statsmodels.regression.mixed_linear_model import MixedLM

data = pd.read_csv("crossed_re_01.csv")

# 拟合包含被试和项目随机截距的线性混合效应模型
model = MixedLM.from_formula('y ~ x', data, re_formula='1', groups=data['subject'],
                             vc_formula={'item': '0 + C(item)'})
result = model.fit(reml=True, method='lbfgs')

# 打印模型摘要
print(result.summary())

# 提取方差成分
subject_var = result.cov_re.iloc[0, 0]  # 被试随机效应方差
item_var = result.vcomp[0]              # 项目随机效应方差
residual_var = result.scale             # 残差方差

print("\nRandom Effects Variance Components:")
print(f"Subject Variance: {subject_var}")
print(f"Item Variance: {item_var}")

print("\nResidual Variance:")
print(f"Residual Variance: {residual_var}")

Python模型关键输出 #

---------------------------------------------------------
Intercept  4.526     0.403  11.225  0.000   3.736   5.316
x          2.943     0.058  50.820  0.000   2.829   3.056
Group Var  4.786                                         
item Var   0.443                                         
========================================================


Random Effects Variance Components:
Subject Variance: 4.786
Item Variance: 0.442

Residual Variance:
Residual Variance: 1.363

三、问题：方差成分结果差异 #

被试随机截距方差（R:4.818 vs Python:4.786）吻合较好，但项目随机截距方差（R:0.778 vs Python:0.442）和残差方差（R:1.030 vs Python:1.363）差异显著，需确认Python代码是否正确。

四、代码修正方案 #

原Python代码中vc_formula={'item': '0 + C(item)'}写法错误，这会将项目作为固定效应处理而非随机截距。正确写法应使用'1'表示项目的随机截距，修正后代码如下：

import numpy as np
import pandas as pd
import statsmodels.api as sm
from statsmodels.regression.mixed_linear_model import MixedLM

data = pd.read_csv("crossed_re_01.csv")

# 修正后的模型拟合：指定项目为随机截距
model = MixedLM.from_formula(
    'y ~ x', 
    data=data, 
    re_formula='1', 
    groups=data['subject'],
    vc_formula={'item': '1'}  # 此处修改为'1'，表示项目的随机截距
)
result = model.fit(reml=True, method='lbfgs')

# 打印模型摘要和方差成分
print(result.summary())

subject_var = result.cov_re.iloc[0, 0]
item_var = result.vcomp[0]
residual_var = result.scale

print("\n修正后方差成分：")
print(f"被试方差: {subject_var:.3f}")
print(f"项目方差: {item_var:.3f}")
print(f"残差方差: {residual_var:.3f}")

修正后，项目随机截距方差和残差方差会更接近R的输出结果，符合模拟数据的真实参数（项目方差真实值为1，残差方差真实值为1）。

内容的提问来源于stack exchange，提问作者Lynchian