Python中statsmodels多元线性假设检验（系数相等）实现求助

背景 #

已通过Statsmodels的MANOVA完成多元全局检验，现需执行线性假设检验（例如检验预测变量pv在两个响应变量se、em上的效应系数相等），Statsmodels暂无直接封装的对应函数，可通过手动构建假设矩阵结合MultivariateOLS的底层统计量实现。

实现步骤 #

1. 提取/拟合多元OLS模型 #

从已有的MANOVA模型中提取多元OLS实例，或直接拟合MultivariateOLS：

import pandas as pd
import numpy as np
from statsmodels.multivariate.manova import MANOVA
from statsmodels.multivariate.multivariate_ols import MultivariateOLS

# 假设df是你的数据集
mod = MANOVA.from_formula('se + em ~ pv + ed + fs + hp', df)
mv_ols = mod.mv_ols  # 提取多元OLS模型

2. 构建假设矩阵 #

要检验β_se(pv) = β_em(pv)，等价于β_se(pv) - β_em(pv) = 0，据此构建假设矩阵H：

# 获取所有系数的名称列表
coef_names = mv_ols.params.index.tolist()
# 定位pv在两个响应变量上的系数索引
pv_se_idx = coef_names.index('pv')
pv_em_idx = coef_names.index('pv', pv_se_idx + 1)

# 构造假设矩阵：一行，对应"se的pv系数 - em的pv系数 = 0"
H = np.zeros((1, len(coef_names)))
H[0, pv_se_idx] = 1
H[0, pv_em_idx] = -1

3. 计算检验统计量 #

基于多元线性假设检验的原理，手动计算Wilks' Lambda、F统计量和p值：

from scipy.stats import f

# 获取模型残差的平方和与交叉乘积矩阵（SSCP）
sscp_resid = mv_ols.mse_resid * (mv_ols.nobs - mv_ols.rank)
# 计算假设对应的SSCP
sscp_hyp = H @ mv_ols.cov_params() @ H.T * mv_ols.nobs

# 计算Wilks' Lambda统计量
wilks_lambda = np.linalg.det(sscp_resid) / np.linalg.det(sscp_resid + sscp_hyp)

# 转换为F统计量并计算p值
p = H.shape[0]  # 假设的维度数
q = mv_ols.endog.shape[1]  # 响应变量数量
n = mv_ols.nobs - mv_ols.rank  # 残差自由度

num_df = p * q
den_df = (n - q + p) * q if (n - q + p) > 0 else n * q - p*(q - p +1)/2

f_stat = (1 - wilks_lambda) / wilks_lambda * (den_df / num_df)
p_value = 1 - f.cdf(f_stat, num_df, den_df)

# 输出结果
print(f"Wilks' Lambda: {wilks_lambda:.6f}")
print(f"F统计量: {f_stat:.6f}")
print(f"p值: {p_value:.6f}")

4. 封装为可复用函数（可选） #

如果需要多次执行不同的假设检验，可封装成函数：

def multivariate_linear_hypothesis(mv_ols, target_vars, response_pairs):
    """
    检验指定预测变量在响应变量对中的系数相等
    
    参数:
        mv_ols: 拟合完成的MultivariateOLS模型
        target_vars: 需要检验的预测变量列表（如['pv']）
        response_pairs: 响应变量对列表（如[('se', 'em')]）
    """
    coef_names = mv_ols.params.index.tolist()
    H_rows = []
    
    for var in target_vars:
        for resp1, resp2 in response_pairs:
            # 定位两个响应变量下对应预测变量的系数索引
            idx1 = coef_names.index(var) if resp1 == mv_ols.endog.columns[0] else coef_names.index(var, coef_names.index(var)+1)
            idx2 = coef_names.index(var) if resp2 == mv_ols.endog.columns[0] else coef_names.index(var, coef_names.index(var)+1)
            
            row = np.zeros(len(coef_names))
            row[idx1] = 1
            row[idx2] = -1
            H_rows.append(row)
    
    H = np.array(H_rows)
    sscp_resid = mv_ols.mse_resid * (mv_ols.nobs - mv_ols.rank)
    sscp_hyp = H @ mv_ols.cov_params() @ H.T * mv_ols.nobs
    
    wilks_lambda = np.linalg.det(sscp_resid) / np.linalg.det(sscp_resid + sscp_hyp)
    p = H.shape[0]
    q = mv_ols.endog.shape[1]
    n = mv_ols.nobs - mv_ols.rank
    
    num_df = p * q
    den_df = (n - q + p) * q if (n - q + p) > 0 else n * q - p*(q - p +1)/2
    
    f_stat = (1 - wilks_lambda) / wilks_lambda * (den_df / num_df)
    p_value = 1 - f.cdf(f_stat, num_df, den_df)
    
    return pd.DataFrame({
        'Wilks\' Lambda': [wilks_lambda],
        'F统计量': [f_stat],
        '分子自由度': [num_df],
        '分母自由度': [den_df],
        'p值': [p_value]
    })

# 使用示例
result = multivariate_linear_hypothesis(mv_ols, ['pv'], [('se', 'em')])
print(result)

关键说明 #

• 该方法基于MultivariateOLS的底层计算结果（协方差矩阵、残差SSCP），完全适配Statsmodels的输出逻辑。
• 对于更复杂的线性假设（如多个变量的组合检验），仅需扩展假设矩阵H的行数即可实现。

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