PLS回归系数获取与验证：coef_调用失效问题及解决方法

问题背景 #

目标是验证Partial Least Squares（PLS）回归模型的系数，回归方程形式为：

y = b + b₁(X₁) + ... + bₙ(Xₙ)
其中y为因变量，b为截距，b₁_{bₙ为系数，X₁}Xₙ为自变量。

用户实现的PLS回归代码如下：

import numpy as np
from sklearn.cross_decomposition import PLSRegression
from sklearn.metrics import mean_squared_error


# Filtrar por municipios
municipios = ['Irapuato', 'Celaya', 'Salamanca']
data = df[df['municipio'].isin(municipios)]

columnas = ['t confirmados', 'imeca (pm2.5)', 'imeca so2', 'imeca co', 'imeca no2', 'imeca (pm 10)' ,'imeca o3']
data = data[columnas]

# Convertir a valores de numpy después de la selección y filtrado
data_values = data.values

# Variable dependiente y variables independientes
y = data_values[:, 0]  # 't confirmados'
X = data_values[:, 1:]  # 'imeca (pm2.5)', 'imeca so2', 'imeca co', 'imeca no2', 'imeca (pm 10)', 'imeca o3'

# Número de componentes para PLS
n_components = 6

# Almacenar RMSEP para cada componente
rmsep = []
best_model = None
best_rmsep = np.inf

for i in range(n_components):
    pls = PLSRegression(n_components=i+1)
    residuals = []
    
    # Validación cruzada dejando uno fuera
    for j in range(len(y)):
        X_train = np.delete(X, j, axis=0)
        y_train = np.delete(y, j, axis=0)
        X_test = X[j].reshape(1, -1)
        y_test = y[j]
        
        pls.fit(X_train, y_train)
        y_pred = pls.predict(X_test)
        
        residuals.append((y_pred - y_test)**2)
    
    # Calcular RMSEP
    rmsep_i = np.sqrt(np.mean(residuals))
    rmsep.append(rmsep_i)
    
    # Guardar el modelo con el RMSEP más pequeño
    if rmsep_i < best_rmsep:
        best_rmsep = rmsep_i
        best_model = pls

# Entrenar el mejor modelo con todos los datos
best_model.fit(X, y)
# Imprimir las predicciones
for i, pred in enumerate(y_pred):
    print(f'Datos de prueba {i+1}: t confirmados predichos = {pred[0]}')
print(best_model.coef_)

问题现象 #

• 调用best_model.predict(new_data)得到的预测结果准确
• 直接使用best_model.coef_输出的系数代入回归方程手动计算，结果与预测值不符
• 手动设置X值反推的系数却能匹配预测结果

原因分析 #

sklearn的PLSRegression默认会对输入的X和y做标准化处理（均值为0，标准差为1），具体逻辑如下：

1. 训练时先计算训练集X的均值x_mean_、标准差x_std_，以及y的均值y_mean_、标准差y_std_
2. 基于标准化后的X和y训练模型，coef_是标准化变量对应的系数
3. predict方法会自动对输入数据先做标准化，计算后再还原为原始y的尺度
4. 直接用原始X乘以coef_加截距，跳过了标准化步骤，所以结果和predict不一致

正确的系数提取与使用方法 #

要手动计算出和predict一致的结果，需要遵循模型内部的标准化逻辑，具体步骤如下：

步骤1：获取模型的标准化参数 #

PLSRegression对象存储了以下关键属性：

• x_mean_: 训练集X的特征均值
• x_std_: 训练集X的特征标准差
• y_mean_: 训练集y的均值
• y_std_: 训练集y的标准差
• intercept_: 模型的截距（基于标准化数据）

步骤2：手动计算匹配预测值 #

按以下流程计算，结果将和predict完全一致：

1. 对输入X做标准化：X_std = (X - x_mean_) / x_std_
2. 计算标准化后的预测值：y_std_pred = X_std @ coef_.T + intercept_
3. 还原为原始尺度的y：y_pred = y_std_pred * y_std_ + y_mean_

修正后的代码示例 #

# 假设已训练好best_model
# 提取标准化参数与模型系数
x_mean = best_model.x_mean_
x_std = best_model.x_std_
y_mean = best_model.y_mean_
y_std = best_model.y_std_
coef = best_model.coef_
intercept = best_model.intercept_

# 取一组测试数据（示例为原始X的第一行）
test_X = X[0].reshape(1, -1)

# 手动计算预测值
X_std = (test_X - x_mean) / x_std
y_std_pred = X_std @ coef.T + intercept
y_pred_manual = y_std_pred * y_std + y_mean

# 对比模型预测结果
y_pred_model = best_model.predict(test_X)

print(f"模型预测值: {y_pred_model[0][0]}")
print(f"手动计算值: {y_pred_manual[0][0]}")

直接获取原始尺度的系数（可选） #

如果想直接得到基于原始数据的回归系数（无需手动标准化），可以通过以下转换计算：

# 原始尺度的系数
coef_original = (coef / x_std) * y_std
# 原始尺度的截距
intercept_original = y_mean - np.sum(coef_original * x_mean) + intercept * y_std

此时直接用y = intercept_original + X @ coef_original.T计算，结果与predict一致。

内容的提问来源于stack exchange，提问作者Carlos Leonel Guerrero Rodrigu