使用statsmodels估计OLS时出现‘Full rank’错误的排查求助

背景与目标 #

拥有某地区作物年产量、年平均气温与年降水量的历史数据，目标是估计线性模型：

$y = \beta_0 + \beta_1 t + \beta_2 t^2 + \beta_3 tmp + \beta_4 tmp^2 + \beta_5 p + \beta_6 p^2 + \beta_7 tmp \times p + \beta_8 tmp^2 \times p^2 + \epsilon$
（其中$y$为作物年产量，$t$代表年份，$tmp$为年平均气温，$p$为年降水量，平方项用于捕捉极端值的影响）

代码与报错 #

使用的Python代码如下：

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

df = pd.read_csv('https://raw.githubusercontent.com/kevinkuranyi/data/main/crop_yield.csv')

model = smf.ols(formula = 'y_banana ~ year+year2+tmp+tmp2+pre+pre2+tmp_pre+tmp2_pre2',
 data=df, missing='drop').fit(cov_type='HAC', cov_kwds={'maxlags': 2})
model.summary()

运行后触发警告：

/usr/local/lib/python3.10/dist-packages/statsmodels/base/model.py:1888: ValueWarning: covariance of constraints does not have full rank. The number of constraints is 8, but rank is 5
  warnings.warn('covariance of constraints does not have full '

用户怀疑是多重共线性问题，但无论剔除哪个变量，只要纳入超过4个变量（即使不含交互项或可能构成线性组合的平方项），仍会出现该错误，已尝试多种变量组合。

问题根源与解决思路 #

1. 核心问题：样本量严重不足 #

先检查数据集规模：运行print(df.shape)会发现，样本观测数远小于模型变量数。比如如果样本只有10-15个观测，却要拟合8个自变量+截距，自由度严重不足，直接导致设计矩阵不满秩，HAC协方差估计无法正常计算——HAC需要足够样本量来估计滞后项的协方差，样本太少时约束矩阵的秩必然不够。

2. HAC协方差估计的额外限制 #

HAC（异方差自相关一致）协方差估计要求样本量远大于变量数和滞后阶数。你设置了maxlags=2，加上8个自变量，至少需要变量数3-5倍以上的样本量，否则会出现协方差矩阵秩不足的问题。

3. 具体解决步骤 #

• 第一步：确认样本规模
  先执行print(df.shape)明确观测数。如果样本量小于20，必须简化模型：
  • 先移除交互项，只保留主效应和平方项；
  • 甚至先从简单的趋势+气温+降水线性模型开始，再逐步添加高阶项。
• 第二步：先放弃HAC，验证模型可行性
  暂时去掉cov_type='HAC', cov_kwds={'maxlags': 2}参数，用普通OLS协方差拟合模型，先确认模型是否能正常运行，再用VIF值（通过statsmodels.stats.outliers_influence.variance_inflation_factor计算）判断是否真的存在多重共线性。
• 第三步：变量中心化缓解共线性
  年份、气温、降水的原始值做平方或交互项时，容易因数值范围差异产生共线性。对变量做中心化处理后再构建高阶项，能有效缓解该问题：

  # 中心化变量
  df['year_c'] = df['year'] - df['year'].mean()
  df['year2_c'] = df['year_c'] ** 2
  df['tmp_c'] = df['tmp'] - df['tmp'].mean()
  df['tmp2_c'] = df['tmp_c'] ** 2
  df['pre_c'] = df['pre'] - df['pre'].mean()
  df['pre2_c'] = df['pre_c'] ** 2
  df['tmp_pre_c'] = df['tmp_c'] * df['pre_c']
  df['tmp2_pre2_c'] = df['tmp2_c'] * df['pre2_c']
  
  # 用中心化后的变量建模
  model = smf.ols(formula='y_banana ~ year_c+year2_c+tmp_c+tmp2_c+pre_c+pre2_c+tmp_pre_c+tmp2_pre2_c',
                  data=df, missing='drop').fit()
  print(model.summary())
  

  若模型能正常拟合，且样本量足够，再考虑添加cov_type='HAC'参数。

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