使用statsmodels拟合MixedLM后执行ANOVA出现维度不匹配问题

问题场景 #

已成功使用statsmodels的MixedLM.from_formula拟合包含虚拟变量的GLIMMIX模型，但执行Type 2 ANOVA时出现维度不匹配错误。

模型代码 #

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

# 补充原代码缺失的K定义
K = 12
T_ctr_vec = np.random.rand(15,12)
T_pat_vec = np.random.rand(15,12)

mask_K = np.ones((K,K)) - np.eye(K,K)
mask_K_ind = np.where(mask_K != 0)

n_pat = T_pat_vec.shape[0]
n_ctr = T_ctr_vec.shape[0]

Y = np.concatenate((T_ctr_vec, T_pat_vec), axis=1).T.flatten()

TT = np.tile(mask_K_ind[0], n_ctr + n_pat)

G_pat = np.full(n_pat, np.nan)
G_pat[0:15] = 3
G = np.concatenate((np.ones(n_ctr), G_pat), axis=0)
G = np.repeat(G, len(mask_K_ind[0]))

S = np.tile(np.arange(1, n_ctr+n_pat+1), (1, len(mask_K_ind[0])))
S = S.flatten()

tbl_gmle_ = pd.DataFrame({'TT':TT,'S':S,'G':G}).astype('category')
tbl_gmle_t_Y = pd.DataFrame({'Y':Y})
tbl_gmle_t = tbl_gmle_.join(tbl_gmle_t_Y)

# 拟合模型
formula = 'Y ~ 1 + TT * G + C(S)'
model = sm.MixedLM.from_formula(formula, data=tbl_gmle_t, groups=tbl_gmle_t["S"])
model.family = sm.families.Gaussian()
result = model.fit()

# 执行Type 2 ANOVA时出错
anova_table = sm.stats.anova_lm(result,typ=2)

错误信息 #

Cell In [814], line 2
      1 # Perform ANOVA
----> 2 anova_table = sm.stats.anova_lm(result,typ=2)
      4 # Print the ANOVA table
      5 print(anova_table)

File /opt/homebrew/Caskroom/miniforge/base/envs/working°env/lib/python3.8/site-packages/statsmodels/stats/anova.py:349, in anova_lm(*args, **kwargs)
    347 if len(args) == 1:
    348     model = args[0]
---> 349     return anova_single(model, **kwargs)
    351 if typ not in [1, "I"]:
    352     raise ValueError("Multiple models only supported for type I. "
    353                      "Got type %s" % str(typ))

File /opt/homebrew/Caskroom/miniforge/base/envs/working°env/lib/python3.8/site-packages/statsmodels/stats/anova.py:80, in anova_single(model, **kwargs)
     77     return anova1_lm_single(model, endog, exog, nobs, design_info, table,
     78                             n_rows, test, pr_test, robust)
     79 elif typ in [2, "II"]:
---> 80     return anova2_lm_single(model, design_info, n_rows, test, pr_test,
     81                             robust)
     82 elif typ in [3, "III"]:
     83     return anova3_lm_single(model, design_info, n_rows, test, pr_test,
     84                             robust)

File /opt/homebrew/Caskroom/miniforge/base/envs/working°env/lib/python3.8/site-packages/statsmodels/stats/anova.py:201, in anova2_lm_single(model, design_info, n_rows, test, pr_test, robust)
    198 L2 = np.eye(model.model.exog.shape[1])[L2]
    200 if L2.size:
---> 201     LVL = np.dot(np.dot(L1,robust_cov),L2.T)
    202     from scipy import linalg
    203     orth_compl,_ = linalg.qr(LVL)

File <__array_function__ internals>:180, in dot(*args, **kwargs)

ValueError: shapes (6,37) and (38,38) not aligned: 37 (dim 1) != 38 (dim 0)

错误原因 #

1. 固定效应与随机效应的冗余冲突：
   公式中的C(S)将分组变量S转为虚拟变量作为固定效应，但同时又通过groups=tbl_gmle_t["S"]指定S为随机效应分组，这会导致设计矩阵的固定效应列数与模型参数的协方差矩阵维度不匹配。原错误中的(6,37)是Type 2 ANOVA计算时的对比矩阵维度，而(38,38)是完整参数的协方差矩阵维度，两者因重复定义S的效应出现维度差异。

2. statsmodels对混合模型ANOVA的支持局限：
   sm.stats.anova_lm最初是为普通线性模型（OLS）设计的，对混合效应模型（MixedLM）的Type 2/3 ANOVA支持不完善，其内部计算逻辑没有考虑混合模型的参数结构（固定效应+随机效应），导致矩阵运算时出现维度对齐失败。

解决方法 #

• 移除固定效应中的C(S)：
  由于groups=S已经为每个S分组添加了随机截距，无需再将C(S)放入固定效应公式，修改后的公式为：

  formula = 'Y ~ 1 + TT * G'
  

  这样可以消除冗余效应，确保设计矩阵与协方差矩阵维度匹配。

• 使用混合模型专用的ANOVA方法：
  若需要对混合模型做Type 2/3 ANOVA，可选择：

  • 手动计算平方和：基于模型的固定效应设计矩阵和参数协方差矩阵，自行实现Type 2平方和的计算逻辑；
  • 切换到更适配的库：比如pymer4（Python），该工具对混合效应模型的ANOVA支持更成熟。

• 验证维度一致性：
  修改模型后，可通过以下代码检查设计矩阵和协方差矩阵的维度是否匹配：

  print("设计矩阵维度:", model.model.exog.shape)
  print("协方差矩阵维度:", result.cov_params().shape)
  

  确保两者的列数/行数一致。

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