sklearn使用主成分回归与偏最小二乘回归

Mike123
577 阅读6分钟

1、导包

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import statsmodels.formula.api as smf
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LinearRegression
import warnings
warnings.filterwarnings('ignore')

2、查看数据

image.png

2.1查看相关系数

corr_matrix = df.corr()
mask = np.triu(np.ones_like(corr_matrix,dtype=bool))
plt.figure(figsize=(15, 12))
sns.heatmap(corr_matrix, mask=mask, annot=True, cmap='Blues', fmt='.2f', linewidths=1, square=True, annot_kws={"size": 9} )
plt.title('Correlation Matrix', fontsize=15)
plt.show()

image.png

可以看出y 与 x1、x2成正相关，与x3、x4成负相关

3.1使用普通最小二乘法

result = smf.ols(formula='y ~x1+x2+x3+x4',data = df).fit()
result.summary()

image.png

  • 查看预测结果
ordinary = pd.DataFrame(pd.concat([pd.Series(y), pd.Series(result.predict(), dtype='float')], axis=1)).reset_index(drop = True)
ordinary.columns = ['y','original_pred']
ordinary
image.png

image.png

普通最小二乘得到的回归方程为y^=62.41+1.55x1+0.51x2+0.10x3-0.14x4，从系数上看可以发现明显不合理的地方，x3与y是负相关，但它的系数确实正的,从t值也可以看出，方程不显著

3.2、使用后退法得到回归方程

import pandas as pd
from statsmodels.formula.api import ols

def backward_select(data, target):
    variate = set(data.columns)
    variate.remove(target)
    selected = list(variate)  # 初始化选择的自变量为所有变量
    current_score, best_new_score = float('inf'), float('inf')
    
    while len(selected) > 1:
        aic_with_variate = []
        
        for candidate in selected:
            remaining = list(set(selected) - set([candidate]))
            formula = "{}~{}".format(target, "+".join(remaining))
            print(formula,candidate)
            aic = ols(formula=formula, data=data).fit().aic
            aic_with_variate.append((aic, candidate))
        
        aic_with_variate.sort(reverse=True)
        print(aic_with_variate)
        best_new_score, worst_candidate = aic_with_variate.pop()
        
        if current_score > best_new_score:
            selected.remove(worst_candidate)
            current_score = best_new_score
            print("AIC is {}, continuing!".format(current_score))
        else:
            print("Feature selection is over!")
            break
    
    formula = "{}~{}".format(target, "+".join(selected))
    print("Final formula is {}".format(formula))
    model = ols(formula=formula, data=data).fit()
    return model

result = backward_select(df,'y')
result.summary()

y~x4+x1+x3 x2
y~x2+x1+x3 x4
y~x2+x4+x3 x1
y~x2+x4+x1 x3
[(67.46828661093299, 'x1'), (62.619952232581625, 'x2'), (61.90359687059506, 'x4'), (61.86628547186257, 'x3')]
AIC is 61.86628547186257, continuing!
y~x4+x1 x2
y~x2+x1 x4
y~x2+x4 x1
[(97.52172751297775, 'x1'), (65.63410626724038, 'x2'), (62.3123927621906, 'x4')]
Feature selection is over!
Final formula is y~x2+x4+x1

image.png 由输出结果可见，逐步回归筛选的最优子集为x1,x2,x4，但在显著性水平为0.05时x4的回归系数不显著，故删除x4再做分析

results = ols(formula='y~x1+x2', data=df).fit()
print(results.summary())
backward_pred = results.predict()
backward = pd.DataFrame(pd.concat([pd.Series(y), pd.Series(backward_pred, dtype='float')], axis=1)).reset_index(drop = True)
backward.columns = ['y','backward_pred']
backward

image.png 从上述输出结果可知，后退法得到的回归方程为y^=52.577+1.468x1+0.662x2。

  • 预测结果 image.png

3.3、使用主成分分析

scaler = StandardScaler()
X_scaler = scaler.fit_transform(np.array(X))#标准化数据
y_scaler = scaler.fit_transform(np.array(y).reshape(-1,1))

model = PCA()
model.fit(X)
#每个主成分能解释的方差
model.explained_variance_
#每个主成分能解释的方差的百分比
model.explained_variance_ratio_

plt.style.use('ggplot')
plt.plot(model.explained_variance_ratio_.cumsum(),'o--')
plt.xlabel('主成分')
plt.ylabel('解释的累积方差比例')
plt.axhline(0.9,color = 'k',linestyle = '--',linewidth = 1)
plt.grid(True)
plt.title('累计PVE')
plt.show()

可以看到当主成分个数为2时就能解释原始数据的90%以上了（坐标轴上是1是因为从0 开始的..）

model = PCA(n_components= 2)
model.fit(X)
X_pca = model.transform(X)
X_pca


array([[ 36.821826  ,   6.87087815],
       [ 29.60727342,  -4.61088196],
       [-12.98177572,   4.20491318],
       [ 23.71472572,   6.63405255],
       [ -0.55319168,   4.46173212],
       [-10.81249083,   3.64657117],
       [-32.58816661,  -8.97984628],
       [ 22.6063955 , -10.72590646],
       [ -9.26258724,  -8.98537335],
       [ -3.28396933,  14.15727734],
       [  9.22003112, -12.38608079],
       [-25.58490852,   2.78169315],
       [-26.90316183,   2.93097117]])

#计算主成分核载矩阵，显示了原始变量和主成分之间的关系。
pca_loadings = pd.DataFrame(model.components_.T,columns = columns,index = X.columns)
pca_loadings

image.png

#打印主成分回归方程
df_pca['y'] = y
all_columns = "+".join(df_pca.columns[:-1])
print("回归方程为:y ~ -0.5537PC1 + 0.9196PC2 + 95.4231")
result = smf.ols(formula='y ~PC1+PC2',data = df_pca).fit()
result.summary()

image.png

3.4以PC1和PC2做因变量，以四个原始自变量做线性回归

data_pca = pd.concat([df_pca,X],axis = 1)
data_pca
linear_model1 = LinearRegression()
linear_model1.fit(data_pca[['x1','x2','x3','x4']],data_pca.PC1)
coef = linear_model1.coef_.squeeze()
intercept_1 = linear_model1.intercept_.squeeze()
print(f'PC1主成分回归方程为 :y = {coef[0]:.2f}X1 + {coef[1]:.2f}X2 + {coef[2]:.2f}X3 + {coef[3]:.2f}X4 + {linear_model1.intercept_.squeeze():.2f}')
linear_model2 = LinearRegression()
linear_model2.fit(data_pca[['x1','x2','x3','x4']],data_pca.PC2)
coef1 = linear_model2.coef_.squeeze()
intercept_2 = linear_model2.intercept_.squeeze()
print(f'PC2主成分回归方程为 :y = {coef1[0]:.2f}X1 + {coef[1]:.2f}X2 + {coef1[2]:.2f}X3 + {coef1[3]:.2f}X4 + {linear_model2.intercept_.squeeze():.2f}')

PC1主成分回归方程为 :y = -0.07X1 + -0.68X2 + 0.03X3 + 0.73X4 + 10.91
PC2主成分回归方程为 :y = 0.65X1 + -0.68X2 + -0.76X3 + 0.11X4 + -0.15
  • 查看还原之后对回归方程
  • 主成分回归方程为:y^=−0.5537PC1+0.9196PC2+95.42307692主成分回归方程为 :\hat{y} = -0.5537PC1 + 0.9196PC2 + 95.42307692
coef_all = (-0.5537* coef +0.9196 * coef1)
intercept_all = -0.5537 * intercept_1 + 0.9196 * intercept_2 + 95.42307692
print(f'主成分回归方程为 :y = {coef_all[0]:.2f}X1 + {coef_all[1]:.2f}X2 + {coef_all[2]:.2f}X3 + {coef_all[3]:.2f}X4 + {intercept_all:.2f}')
主成分回归方程为 :y = 0.63X1 + 0.39X2 + -0.71X3 + -0.30X4 + 89.25

3.5偏最小二乘法回归

# 现用y对前两个主成分做偏最小二乘回归
from sklearn.cross_decomposition import PLSRegression
pls_model = PLSRegression()
pls_model.fit(X,y)
pls_pred = pls_model.predict(X).squeeze()
pls = pd.DataFrame(pd.concat([pd.Series(y), pd.Series(pls_pred, dtype='float')], axis=1)).reset_index(drop = True)
pls.columns = ['y','pls_pred']
pls

4查看每个方法预测值

#查看每一个方法的预测值
merge1 = pd.merge(ordinary, backward, how='left')
merge2 = pd.merge(merge1, pca, on='y', how='left')
merge3 = pd.merge(merge2, pls, on='y', how='left')
all_pred = pd.merge(merge3, pca_all, on='y', how='left')
all_pred

image.png

#查看每一个方法与真实值的差距
column_diffs = all_pred.subtract(all_pred['y'], axis=0)
diff = column_diffs.T
diff['MAE'] = (np.abs(column_diffs).sum().sort_values())
diff = diff[['MAE'] + [col for col in diff.columns if col != 'MAE']]  # 重新排序列
diff.sort_values(by = 'MAE')

image.png

总结:

  • 偏最小二乘回归的预测效果最好，其次是后退法，接着是还原后的主成分回归方程，最后是普通最小二乘回归和使用PCA1和PCA2的主层次回归

    • pls_pred的平均绝对误差为 22.043398。
    • backward_pred的平均绝对误差为 24.821288。
    • pca_all_pred的平均绝对误差为 33.824407。
    • original_pred和pca_pred的平均绝对误差相同，都是 33.828070。

  • 每个模型预测的表现：

    • 在样本0和样本2中，pls_pred的误差比其他模型小。
    • backward_pred在样本1和样本8上表现相对较好。
    • pca_all_pred在样本3上表现相对较好。
    • original_pred和pca_pred在所有样本上的表现相似。

  • 模型之间的差异：

    • pls_pred相对其他模型在整体上表现更为稳定，因为它的平均绝对误差较小。
    • backward_pred在某些样本上表现优异，但整体上略逊于pls_pred。
    • pca_all_pred在平均绝对误差方面显示最大的差异。
    • original_pred和pca_pred之间没有明显的差异。

  • 总体而言，pls_pred可能是最优选择，因为它在整体表现上具有较小的平均绝对误差，并且在某些样本上的表现优于其他模型。

  • 偏最小二乘回归（PLSR）： PLSR 是一种基于线性模型的多元统计方法，旨在通过同时考虑自变量和因变量的变化来找到它们之间的最大协方差方向。 PLSR 将自变量和因变量同时投影到低维空间，找到最能解释因变量方差的方向，并基于这些方向建立模型。 PLSR 的目标是找到称为潜在变量（latent variables）或者部分最小二乘（partial least squares）的新特征，以最大程度地解释两个变量之间的协方差。 PLSR 的预测变量是根据自变量和因变量之间的共同方差来选择的。

  • 主成分回归（PCR）：

    PCR 首先使用主成分分析（PCA）来降低自变量的维度，将它们投影到低维的主成分空间。 然后，使用这些主成分作为新的自变量来建立线性回归模型。 PCR 的目标是通过降低自变量的维度来减少多重共线性，并使用较少数量的主成分来建立模型。

avatar
文章
阅读
粉丝