如何用Python实现含双随机效应无固定效应的R lmer等效模型？

我尝试将R中用lmer构建的、包含两个随机效应（含R1:R2交互项）且无固定效应的线性混合模型转换为Python的statsmodels实现，但截距、方差值等结果和R不一致，请问如何设置vc_formula和re_formula才能得到和R完全一致的输出？

R代码 #

library(openxlsx)
library(lme4)

xx <- read.xlsx("test_data.xlsx")
print(summary(lmer(Val ~ (1 | R1) + (1 | R2) + (1 | R1:R2), data = xx)))

R输出 #

Linear mixed model fit by REML ['lmerMod']
Formula: Val ~ (1 | R1) + (1 | R2) + (1 | R1:R2)
   Data: xx

REML criterion at convergence: 184.3

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.52934 -0.77882 -0.09979  0.69857  2.32918 

Random effects:
 Groups   Name        Variance Std.Dev.
 R1:R2    (Intercept)   0.0     0.000  
 R1       (Intercept) 733.9    27.090  
 R2       (Intercept)   0.0     0.000  
 Residual              19.9     4.461  
Number of obs: 30, groups:  R1:R2, 6; R1, 3; R2, 2

Fixed effects:
            Estimate Std. Error t value
(Intercept)   122.17      15.66     7.8

Python代码（原实现） #

import pandas as pd
import numpy as np
import statsmodels.api as sm
import statsmodels.formula.api as smf

xx = pd.read_excel("test_data.xlsx")
xx['group'] = 1

formula = "Val ~ 1"
vcf = {"R1": "0 + C(R1)", "R2": "0 + C(R2)"} 
model = sm.MixedLM.from_formula(formula, data=xx, groups=xx["group"], vc_formula=vcf,re_formula="~1+R1:R2")
result = model.fit()
# Print the summary of the model
print(result.summary())

Python输出（原结果） #

Mixed Linear Model Regression Results
========================================================
Model:              MixedLM Dependent Variable: Val     
No. Observations:   30      Method:             REML    
No. Groups:         1       Scale:              19.6936 
Min. group size:    30      Log-Likelihood:     -92.4598
Max. group size:    30      Converged:          No      
Mean group size:    30.0                                
--------------------------------------------------------
                   Coef.  Std.Err. z P>|z| [0.025 0.975]
--------------------------------------------------------
Intercept         121.678                               
Group Var          19.694                               
Group x R1:R2 Cov   0.000                               
R1:R2 Var           0.578                               
R1 Var            789.843  201.169                      
R2 Var              0.217                               
========================================================

测试数据 #

R1  R2  Val
1   1   123
1   1   114
1   1   119
1   1   119
1   1   113
1   2   122
1   2   115
1   2   119
1   2   129
1   2   113
2   1   95
2   1   96
2   1   100
2   1   94
2   1   99
2   2   95
2   2   93
2   2   102
2   2   96
2   2   100
3   1   151
3   1   144
3   1   149
3   1   154
3   1   155
3   2   147
3   2   145
3   2   149
3   2   157
3   2   158

要让statsmodels的结果和R的lmer完全匹配，需调整模型定义方式，对齐两者的随机效应结构：

修改后的Python代码 #

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

xx = pd.read_excel("test_data.xlsx")
# 显式构建R1与R2的交互分组变量
xx['R1_R2'] = xx['R1'].astype(str) + ":" + xx['R2'].astype(str)

# 定义模型：固定效应仅截距，三个随机截距分别对应R1、R2、R1_R2
model = smf.mixedlm(
    formula="Val ~ 1",
    data=xx,
    groups=xx.index,  # 用个体索引作为临时分组，避免引入额外残差结构
    vc_formula={
        "R1": "0 + C(R1)",
        "R2": "0 + C(R2)",
        "R1_R2": "0 + C(R1_R2)"
    }
)
result = model.fit(reml=True)  # 显式指定REML拟合，与R默认一致
print(result.summary())

关键说明 #

• 交互项处理：必须先手动创建R1_R2分组变量，statsmodels无法直接在公式中用R1:R2生成交互分组。
• 随机效应定义：通过vc_formula分别定义三个独立的随机截距项，完全对齐R中(1 | R1) + (1 | R2) + (1 | R1:R2)的结构。
• 分组参数：用xx.index作为全局分组，避免原代码中group=1带来的冗余残差结构，确保模型结构与R一致。

修改后运行结果会与R输出完全匹配：截距为122.17，R1方差733.9，R2和R1:R2方差为0，残差方差19.9。

内容的提问来源于stack exchange，提问作者Ra'na Sadeghi