机器学习之主成分分析建模

        大家好，我是带我去滑雪！

        本期介绍一种常见的非监督学习方法，即主成分分析。对于非监督学习，其数据中只含有特征变量x，而没有响应变量y。因此非监督学习的目标并非用x预测y，而是探索特征变量x本身的规律和模式。主成分分析是统计学中进行降维的经典方法，该降维的思想是将多个高度相关的特征变量转换为几个不存在线性关系的特征变量，转置后的变量称为主成分，每个主成分都是原始特征变量的线性组合，主成分可以反映原始数据的大部分信息，从而达到简化系统结构，抓住问题实质。一般特征变量个数太多或存在严重的多重共线性时，可以使用主成分分析对自变量进行处理，在科研中一般将主成分分析作为研究中的一个中间环节。

目录

1、主成分分析建模

（1）导入相关模块与数据

（2）计算特征变量之间的pearson相关系数，并绘制相关系数热力图

（3）数据标准化

（4）主成分拟合

（5）计算主成分载荷矩阵

（6）画图展示前4个主成分的载荷向量

（7）计算主成分得分

（8）利用K均值聚类对三个主成分聚类并进行可视化

2、主成分回归建模

（1）导入相关库与数据

（2）计算中国香港经济增长率与其他国家或地区经济增长率的相关系数

（3）划分训练集与测试集，并进行标准化

（4）使用留一交叉验证选择误差最小的时候的主成分个数

（5）主成分回归

---

1、主成分分析建模

        使用一个关于听觉测试的数据集audiometric.csv，该数据集包含100个观测值与8个特征变量，分别为100位39岁男子左右耳在四种不同频率上的听力度量。数据集可在文末获取。

（1）导入相关模块与数据

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from sklearn.model_selection import cross_val_score
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import LeaveOneOut
from mpl_toolkits import mplot3d
 
audiometric = pd.read_csv(r'E:\工作\硕士\博客\博客22-主成分分析与主成分回归\audiometric.csv')#读取数据
audiometric.shape

输出结果：

(100, 8)

audiometric.head()

输出结果：

	lft500	lft1000	lft2000	lft4000	rght500	rght1000	rght2000	rght4000
0	0	5	10	15	0	5	5	15
1	-5	0	-10	0	0	5	5	15
2	-5	0	15	15	0	0	10	25
3	-5	0	-10	-10	-10	-5	-10	10
4	-5	-5	-10	10	0	-10	-10	50

        其中lft500表示左耳在500Hz频率上的听力，而变量rht500表示右耳在500Hz频率上的听力，其他特征变量以次类推。

（2）计算特征变量之间的pearson相关系数，并绘制相关系数热力图

pd.options.display.max_columns = 10
round(audiometric.corr(), 2)#使用round函数将小数点保留两位

输出结果：

	lft500	lft1000	lft2000	lft4000	rght500	rght1000	rght2000	rght4000
lft500	1.00	0.78	0.40	0.26	0.70	0.64	0.24	0.20
lft1000	0.78	1.00	0.54	0.27	0.55	0.71	0.36	0.22
lft2000	0.40	0.54	1.00	0.43	0.24	0.45	0.70	0.33
lft4000	0.26	0.27	0.43	1.00	0.18	0.26	0.32	0.71
rght500	0.70	0.55	0.24	0.18	1.00	0.66	0.16	0.13
rght1000	0.64	0.71	0.45	0.26	0.66	1.00	0.41	0.22
rght2000	0.24	0.36	0.70	0.32	0.16	0.41	1.00	0.37
rght4000	0.20	0.22	0.33	0.71	0.13	0.22	0.37	1.00

sns.heatmap(round(audiometric.corr(), 2),cmap='Blues',annot=True)#annot=True表示热力图上展示数值
plt.savefig("squares3.png",
            bbox_inches ="tight",
            pad_inches = 1,
            transparent = True,
            facecolor ="w",
            edgecolor ='w',
            dpi=300,
            orientation ='landscape')

输出结果：

       结果显示，这8个特征变量之间相关性较高。

（3）数据标准化

scaler = StandardScaler()#创建StandardScaler类的一个实例
scaler.fit(audiometric)#使用fit方法进行估计
X = scaler.transform(audiometric)#使用transform（）方法完成标准化

（4）主成分拟合

model = PCA()#创建一个主成分类的实例
model.fit(X)#模型拟合
model.explained_variance_#每个主成分能解释的方差

输出结果：

array([3.96869223, 1.63466846, 0.98517655, 0.47149715, 0.34352524,
       0.31908202, 0.20213246, 0.15603398])

plt.plot(model.explained_variance_ratio_, 'o-')#绘制线图
plt.xlabel('Principal Component')
plt.ylabel('Proportion of Variance Explained')
plt.title('PVE')
plt.savefig("squares4.png",
            bbox_inches ="tight",
            pad_inches = 1,
            transparent = True,
            facecolor ="w",
            edgecolor ='w',
            dpi=300,
            orientation ='landscape')

输出结果：

         绘制累计方差解释百分比图，直观判断选几个主成分合适：

plt.plot(model.explained_variance_ratio_.cumsum(), 'o-')
plt.xlabel('Principal Component')
plt.ylabel('Cumulative Proportion of Variance Explained')
plt.axhline(0.9, color='k', linestyle='--', linewidth=1)
plt.title('Cumulative PVE')
plt.savefig("squares5.png",
            bbox_inches ="tight",
            pad_inches = 1,
            transparent = True,
            facecolor ="w",
            edgecolor ='w',
            dpi=300,
            orientation ='landscape')

输出结果：

         通过累积方差解释百分比图可以看出，选取前4个主成分，所能解释的累积方差解释比重占86%以上。

（5）计算主成分载荷矩阵

model.components_#计算主成分核载矩阵
columns = ['PC' + str(i) for i in range(1, 9)]#为更好展示主成分载荷矩阵，生成一个带有主成分名称的列表
pca_loadings = pd.DataFrame(model.components_.T, index=audiometric.columns,columns=columns)
round(pca_loadings, 2)#保留两位小数

输出结果：

	PC1	PC2	PC3	PC4	PC5	PC6	PC7	PC8
lft500	0.40	0.32	-0.16	-0.33	0.02	-0.45	-0.33	-0.55
lft1000	0.42	0.23	0.05	-0.48	-0.38	0.07	0.03	0.62
lft2000	0.37	-0.24	0.47	-0.28	0.44	0.06	0.53	-0.19
lft4000	0.28	-0.47	-0.43	-0.16	0.35	0.42	-0.43	0.08
rght500	0.34	0.39	-0.26	0.49	0.50	-0.19	0.16	0.34
rght1000	0.41	0.23	0.03	0.37	-0.35	0.61	0.08	-0.36
rght2000	0.31	-0.32	0.56	0.39	-0.11	-0.27	-0.48	0.15
rght4000	0.25	-0.51	-0.43	0.16	-0.40	-0.37	0.41	-0.05

        结果表明，第1个主成分PC1对于所有变量的载荷系数均为正，且大致相同，故第1个主成分可解释为耳朵的整体听力。第2个主成分PC2对于较低频率的听力（频率为500HZ和1000HZ）的载荷矩阵则为正，而对于较高频率的听力（2000HZ和4000HZ）的载荷系数为负，故第2个主成分为较低频与较高频听力的对比。总之，对于主成分的含义解释，一般可通过观察主成分载荷系数来进行，带有一定的主观性。

（6）画图展示前4个主成分的载荷向量

fig, ax = plt.subplots(2, 2)
plt.subplots_adjust(hspace=1, wspace=0.5)   
for i in range(1, 5):
    ax = plt.subplot(2, 2, i)
    ax.plot(pca_loadings['PC' + str(i)], 'o-')
    ax.axhline(0, color='k', linestyle='--', linewidth=1)
    ax.set_xticks(range(8))
    ax.set_xticklabels(audiometric.columns, rotation=30)
    ax.set_title('PCA Loadings for PC' + str(i))
    plt.savefig("squares6.png",
            bbox_inches ="tight",
            pad_inches = 1,
            transparent = True,
            facecolor ="w",
            edgecolor ='w',
            dpi=300,
            orientation ='landscape')

输出结果：

（7）计算主成分得分

pca_scores = model.transform(X)#计算主成分得分
pca_scores = pd.DataFrame(pca_scores, columns=columns)#将计算的主成分得分存入数据框中
pca_scores.shape#查看主成分得分数据框形状
pca_scores
输出结果：

	PC1	PC2	PC3	PC4	PC5	PC6	PC7	PC8
0	1.186389	0.684794	0.732349	-0.046970	-0.100041	0.266155	0.225687	-0.076763
1	-0.296519	1.089376	0.287918	1.172397	-0.942729	0.133382	-0.175504	0.216139
2	0.738238	-0.443551	1.098123	0.401944	0.345462	-0.215409	0.594275	0.189438
3	-2.142758	1.071569	0.025785	-0.694504	-0.989423	-0.244314	0.354743	0.032974
4	-1.444607	-0.269376	-1.723767	0.186417	-0.224337	-1.379977	0.915203	0.367099
...	...	...	...	...	...	...	...	...
95	0.756276	1.347916	0.482331	0.320696	0.325746	-1.097455	0.078947	-0.371867
96	-1.857851	-0.886236	-0.260644	-0.641353	1.101487	1.266948	0.091726	-0.393849
97	3.308567	-1.533835	-1.497721	-1.136666	-0.354971	-1.128547	0.323248	-0.308402
98	0.863854	-0.515402	0.159237	0.161329	-1.095811	0.354911	0.558742	0.151016
99	-1.014216	0.518773	1.353567	1.320252	-0.962398	-0.661237	-0.781416	-0.026140

100 rows × 8 columns

 

sns.scatterplot(x='PC1', y='PC2', data=pca_scores)#画出第1个主成分得分与第2个主成分得分的散点图
plt.title('Biplot')
plt.savefig("squares7.png",
            bbox_inches ="tight",
            pad_inches = 1,
            transparent = True,
            facecolor ="w",
            edgecolor ='w',
            dpi=300,
            orientation ='landscape')

输出结果：

 

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(pca_scores['PC1'], pca_scores['PC2'], pca_scores['PC3'], c='b')#画前3个主成分得分的三维散点图
ax.set_xlabel('PC1')
ax.set_ylabel('PC2')
ax.set_zlabel('PC3')
plt.savefig("squares8.png",
            bbox_inches ="tight",
            pad_inches = 1,
            transparent = True,
            facecolor ="w",
            edgecolor ='w',
            dpi=300,
            orientation ='landscape')

输出结果：

 

（8）利用K均值聚类对三个主成分聚类并进行可视化

from sklearn.cluster import KMeans
model = KMeans(n_clusters=3, random_state=1, n_init=20)
model.fit(X)
model.labels_
 
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(pca_scores['PC1'], pca_scores['PC2'], pca_scores['PC3'],
           c=model.labels_, cmap='rainbow')
ax.set_xlabel('PC1')
ax.set_ylabel('PC2')
ax.set_zlabel('PC3')
plt.savefig("squares9.png",
            bbox_inches ="tight",
            pad_inches = 1,
            transparent = True,
            facecolor ="w",
            edgecolor ='w',
            dpi=300,
            orientation ='landscape')

输出结果：

         通过图像可见，三种不同类别的散点，分布于三维空间的不同区域。

2、主成分回归建模

         使用中国香港经济增长率的季度数据集growth.csv，该数据集最初用于研究中国香港与内地经济整合之后的处理效应。2003年9月，中国内地与香港签署《内地与香港关于建立更紧密贸易关系的安排》协议，并于2004年1月1日生效。Hsiao et al.（2012）利用与中国香港相邻或有密切贸易关系的24个国家或地区的经济增长率，预测中国香港如果未与内地经济整合的“反事实结果”。其中Quarter为时间，HongKong_CN为响应变量，即中国香港的经济增长率，数据可在文末下载。

（1）导入相关库与数据

 import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from sklearn.model_selection import cross_val_score
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import LeaveOneOut
from mpl_toolkits import mplot3d
growth=pd.read_csv(r'E:\工作\硕士\博客\博客22-主成分分析与主成分回归\growth.csv')#读取数据
growth

输出结果：

	Quarter	HongKong_CN	Australia	Austria	Canada	...	Indonesia	Malaysia	Thailand	Taiwan_CN	China
0	1993Q1	0.062	0.040489	-0.013084	0.010064	...	0.064024	0.085938	0.080000	0.064902	0.1430
1	1993Q2	0.059	0.037857	-0.007581	0.021264	...	0.066068	0.131189	0.080000	0.065123	0.1410
2	1993Q3	0.058	0.022509	0.000543	0.018919	...	0.057959	0.109666	0.080000	0.067379	0.1350
3	1993Q4	0.062	0.028747	0.001181	0.025317	...	0.062365	0.075801	0.080000	0.069164	0.1350
4	1994Q1	0.079	0.033990	0.025511	0.043567	...	0.049743	0.049147	0.112509	0.069451	0.1250
...	...	...	...	...	...	...	...	...	...	...	...
56	2007Q1	0.055	0.058013	0.036198	0.030712	...	0.095297	0.039884	0.049846	0.041019	0.1110
57	2007Q2	0.062	0.059519	0.032570	0.039827	...	0.110900	0.080276	0.051197	0.051073	0.1167
58	2007Q3	0.068	0.056649	0.031558	0.034742	...	0.110100	0.093361	0.050126	0.066369	0.1002
59	2007Q4	0.069	0.045825	0.019095	0.038128	...	0.110700	0.151739	0.079587	0.062929	0.1017
60	2008Q1	0.073	0.027523	0.017431	0.029217	...	0.133400	0.167034	0.049144	0.058841	0.1238

61 rows × 26 columns

（2）计算中国香港经济增长率与其他国家或地区经济增长率的相关系数

growth.index = growth['Quarter']#设置时间索引
growth = growth.drop(columns=['Quarter'])#去掉原来的Quarter变量
growth.corr().iloc[:, 0]#计算香港和其他地区的相关系数

输出结果：

HongKong_CN      1.000000
Australia        0.135181
Austria          0.119270
Canada           0.490202
Denmark          0.396235
Finland         -0.219702
France          -0.121903
Germany         -0.036767
Italy           -0.063108
Japan            0.551349
Korea            0.487735
Mexico           0.450608
Netherlands      0.064113
NewZealand       0.340908
Norway           0.709428
Switzerland      0.214031
UnitedKingdom   -0.119820
UnitedStates     0.086041
Singapore        0.742841
Philippines      0.577118
Indonesia        0.609668
Malaysia         0.756037
Thailand         0.655752
Taiwan_CN        0.414107
China            0.479530
Name: HongKong_CN, dtype: float64

         结果显示，有不少国家或地区的经济增长率与中国香港有关。

（3）划分训练集与测试集，并进行标准化

        由于中国香港与内地经济整合政策从2004Q1开始生效，故使用1993Q1-2003Q4，共44个数据作为训练集，样本的特征变量为24个，即24个国家或地区的经济增长率。由于特征变量个数相对于样本数量过多，如果直接使用线性回归容易导致过拟合。接下来，我们使用主成分回归来预测中国香港经济增长率。定义前44个样本为训练集X_train(变量HongKong_CN)。

X_train = growth.iloc[:44, 1:]#定义特征变量的训练集
X_test = growth.iloc[44:, 1:]#定义特征变量的测试集
y_train = growth.iloc[:44, 0]#定义响应变量的训练集
y_test = growth.iloc[44:, 0]#定义响应变量测试集
 
scaler = StandardScaler()
scaler.fit(X_train)
X_train = scaler.transform(X_train)#特征变量训练集标准化
X_test = scaler.transform(X_test)#特征变量测试集标准化

（4）使用留一交叉验证选择误差最小的时候的主成分个数

scores_mse = []
for k in range(1, 24):
    model = PCA(n_components=k)
    model.fit(X_train)
    X_train_pca = model.transform(X_train)
    loo = LeaveOneOut()
    mse = -cross_val_score(LinearRegression(), X_train_pca, y_train, 
                           cv=loo, scoring='neg_mean_squared_error')
    scores_mse.append(np.mean(mse))
min(scores_mse)
 
index = np.argmin(scores_mse)
index
 
plt.plot(range(1, 24), scores_mse)
plt.axvline(index + 1, color='k', linestyle='--', linewidth=1)
plt.xlabel('Number of Components')
plt.ylabel('Mean Squared Error')
plt.title('Leave-one-out Cross-validation Error')
plt.tight_layout()
plt.savefig("squares10.png",
            bbox_inches ="tight",
            pad_inches = 1,
            transparent = True,
            facecolor ="w",
            edgecolor ='w',
            dpi=300,
            orientation ='landscape')

输出结果：

       通过图像可以发现主成分个数为6时最小，下面进行仅保留6个主成分的主成分分析。

（5）主成分回归

model = PCA(n_components = index + 1)
model.fit(X_train)
#得到主成分得分
X_train_pca = model.transform(X_train)
X_test_pca = model.transform(X_test)
X_train_pca#训练集的主成分得分
 
#进行线性回归拟合
reg = LinearRegression()
reg.fit(X_train_pca, y_train)
 
#全样本预测
X_pca = np.vstack((X_train_pca, X_test_pca))#将训练集与测试集的主成分得分以垂直方式合并
pred = reg.predict(X_pca)
 
y = growth.iloc[:, 0]#中国香港经济实际增长率
 
#将中国香港实际增长率与预测增长率进行可视化
plt.figure(figsize=(10, 5))
ax = plt.gca()
plt.plot(y, label='Actual', color='k')
plt.plot(pred, label='Predicted', color='k', linestyle='--')
plt.xticks(range(1, 62))
ax.set_xticklabels(growth.index, rotation=90)
plt.axvline(44, color='k', linestyle='--', linewidth=1)#添加辅助线
plt.xlabel('Quarter')
plt.ylabel('Growth Rate')
plt.title("Economic Growth of HongKong_CN")
plt.legend(loc='upper left')
plt.tight_layout()
plt.savefig("squares11.png",
            bbox_inches ="tight",
            pad_inches = 1,
            transparent = True,
            facecolor ="w",
            edgecolor ='w',
            dpi=300,
            orientation ='landscape')

输出结果：

        从图可见，在2004Q1政策实施之前，主成分回归的预测增长率较好地拟合了中国香港的实际增长率。而在2004Q1政策实施后，二者开始背离，中国香港实际增长率一直高于预测增长率，而二者的差距即为对该政策处理效应的估计。

---

需要数据集的家人们可以去百度网盘（永久有效）获取：

链接：https://pan.baidu.com/s/1qAu4ueR_ucJO_WbOKAYg2w?pwd=7woc 
提取码：7woc 
--来自百度网盘超级会员V5的分享

---

更多优质内容持续发布中，请移步主页查看。

   点赞+关注,下次不迷路！