一、实验背景

因子分析法是一种寻找公共因子的模型分析方法，其目的是用少数几个因子去描述许多指标或因素之间的联系，将联系比较密切的几个因子变量归为同一类，每一类变量即为一个因子，用少数几个因子反映大部分的信息。运用这种模型方法，我们可以很方便的找出影响原有变量的主要因素有哪些。各省会城市通常是各省的经济、政治、文化中心，带动周边经济发展，是该省份其他地区经济和社会发展的“引路者”，由此吸引了很多人口到省会城市工作、定居。各个城市的常住人口的收入、生活便利情况受到很多因素的影响，如平均工资、房价、储蓄、医院情况等。通过因子分析模型，我们可以将这些指标进行归类，从而将影响该城市常住人口生活水平的指标进行简化。

二、实验内容及数据

2.1 概述

本项目选择各省份中心城市2018年的相关数据进行分析，数据来自于中国统计年鉴。总共选取了7个因子指标，分别是国内生产总值(亿元)、在岗职工平均工资(元)、住宅商品房平均销售价格(元/平方米)、城乡居民储蓄年末余额(亿元)、医院数(个)、社会商品零售总额(亿元)、年末总人口(万人)。

2.2 变量介绍

下面对原始数据中的变量进行相关说明，变量说明如表所示。

国内生产总值(亿元)、在岗职工平均工资(元)、住宅商品房平均销售价格(元/平方米)、城乡居民储蓄年末余额(亿元)、医院数(个)、社会商品零售总额(亿元)、年末总人口(万人)。

变量名	说明
GDP	国内生产总值（亿元）

AW	在岗职工平均工资(元)
APH	住宅商品房平均销售价格(元/平方米)
HSB	城乡居民储蓄年末余额(亿元)
NH	医院数(个)

TRS	社会商品零售总额(亿元)
TP	年末总人口(万人)

三、实验步骤

3.1 导入模块和数据

先导入我们所需要的模块：

import pandas as pd    
import numpy as np
from csv import reader 
from sklearn import preprocessing
from sklearn.decomposition import PCA
from factor_analyzer import *
from scipy.stats import bartlett
import matplotlib.pyplot as plt
import matplotlib
import math as math
from matplotlib import cm
import numpy.linalg as nlg

然后，对要使用到的算法步骤进行定义。为了使结果更为直观，我们可以将最后对样本的评分结果用条形图来表示，那么我们需要对这一过程进行定义，如下所示：

#画条形图
def draw_hist(myname,mydata):
    fig = plt.figure(figsize=(10,5))
    plt.rcParams['font.sans-serif'] = ['SimHei'] 
    plt.rcParams['axes.unicode_minus'] = False
    ax = fig.add_subplot(111)
    index=np.arange(0,len(myname),1)
    ax.bar(index,mydata,width=[0.5])
    xticks = range(0,len(myname), 1)
    xlabels = [el for el in myname]
    ax.set_xticks(xticks)
    ax.set_xticklabels(myname,rotation=45)
    ax.set_xlabel("城市")
    ax.set_ylabel("综合得分")
    plt.title("常住人口生活水平")
    plt.show()

其次，是定义KMO检验法，这一检验法可以帮助判断我们所选择的数据是否适合做因子分析。通常来说，KMO在0.9以上，非常合适做因子分析；在0.8-0.9之间，很适合；在0.7-0.8之间，适合；在0.6-0.7之间，尚可；在0.5-0.6之间，表示很差；在0.5以下应该放弃。

# KMO测度
def kmo(dataset_corr):
    corr_inv = np.linalg.inv(dataset_corr)
    nrow_inv_corr, ncol_inv_corr = dataset_corr.shape
    A = np.ones((nrow_inv_corr, ncol_inv_corr))
    for i in range(0, nrow_inv_corr, 1):
        for j in range(i, ncol_inv_corr, 1):
            A[i, j] = -(corr_inv[i, j]) / (math.sqrt(corr_inv[i, i] * corr_inv[j, j]))
            A[j, i] = A[i, j]
    dataset_corr = np.asarray(dataset_corr)
    kmo_num = np.sum(np.square(dataset_corr)) - np.sum(np.square(np.diagonal(A)))
    kmo_denom = kmo_num + np.sum(np.square(A)) - np.sum(np.square(np.diagonal(A)))
    kmo_value = kmo_num / kmo_denom
    return kmo_value

另外，我们在相关系数矩阵的特征根和特征向量之后，需要对其进行排序，所以，还需要定义特征根和特征向量排序的步骤。

# 特征根与特征向量排序
def sort_valvector(Teig_value, Teigvector):
    Temp=[]
    Temp.append(Teig_value)
    TempNames=["a"]
    for i in range(Teigvector.shape[0]):
          TempNames.append("L"+str(i+1))
    TempMat=np.append(Temp,Teigvector,axis=0).T
    PDTempMat=pd.DataFrame(TempMat)
    PDTempMat.columns=TempNames
    PDTempMat.sort_values('a', ascending=False,inplace=True)
    Teig_value=np.array(PDTempMat['a']).T
    Teigvector=np.array(PDTempMat.iloc[0:,1:]).T
    return Teig_value,Teigvector

最后，导入数据集，删除数据的第一列“城市”，只保留需要分析的数据。

#导入数据
data1 = pd.read_csv('data.csv', sep=',')
myID=data1["城市"]
data=data1
del data["城市"]
ColName=data.columns
print(ColName)
data=data.astype(float)

3.2 数据分析

3.2.1 数据标准化

首先需要对数据进行标准化处理，计算出每列数据的平均值和标准差，利用标准化公式将每个数据标准化，计算得到标准化之后的数据，用dataz表示。

#数据标准化
data=np.array(data)
print(data.shape[0])
tempavg=data.mean(axis=0)  #计算每列的平均值
tempdev=data.var(axis=0)*data.shape[0]/(data.shape[0]-1)   #计算每列的标准差 
tempdev=np.sqrt(tempdev)
dataz = [[0 for i in range(data.shape[1])] for i in range(data.shape[0])]  #定义标准化矩阵
for i in range(data.shape[0]):
    for j in range(data.shape[1]):
        dataz[i][j]=(data[i][j]-tempavg[j])/tempdev[j]
print(np.round(dataz,2))

标准化后结果如下：

3.2.2 相关系数矩阵

对标准化后的数据求相关矩阵：

#求相关矩阵
pddataz=pd.DataFrame(dataz)
print(np.round(pddataz.corr(),2))
DF_corr=pddataz.corr()

最后得到一个7×7的矩阵：

将各个因子变量用热力图来表示：

#热力图
cmap = cm.Blues
fig=plt.figure()
ax=fig.add_subplot(111)
map = ax.imshow(DF_corr, interpolation='nearest', cmap=cmap, vmin=0, vmax=1)
plt.title('correlation coefficient--headmap')
ax.set_yticks(range(len(DF_corr.columns)))
ax.set_yticklabels(DF_corr.columns)
ax.set_xticks(range(len(DF_corr)))
ax.set_xticklabels(DF_corr.columns)
plt.colorbar(map)
plt.show()

图如下：

3.2.3 KMO方法检测

调用前面定义的KMO检测法，对数据进行检测，从而判断该数据是否适用因子分析模型进行分析。

#KMO测度
from factor_analyzer.factor_analyzer import calculate_bartlett_sphericity
chi_square_value,p_value=calculate_bartlett_sphericity(DF_corr)
print("p值:",np.round(p_value,3))
print("KMO测度:", kmo(DF_corr))

最后得到的结果如下图

可知p值<0.05，KMO测度值大于0.7，该数据很适合做因子分析，因此，我们可以用因子分析法对这些因素进行分析。

3.2.4 特征值和特征向量

求相关系数矩阵的特征值及其对应的特征向量

#求特征根
R= np.matrix(pddataz.corr())  
eig_value, eigvector = np.linalg.eig(R)
print(np.round(eig_value,3))
print(np.round(eigvector,3))

求得结果如下：

之后，按照由大到小的顺序对特征值进行排序，对应特征向量相应改变位置：

# 特征根与特征向量排序
eig_value, eigvector=sort_valvector(eig_value, eigvector)
print(np.round(eig_value,3))
print(np.round(eigvector,3))

结果如下：

3.2.5 求因子载荷矩阵

首先，我们要确定公共因子的个数。由上一步骤可以看到，前三个特征值的比重大于85%的标准，所以，公共因子的个数为3。之后，再求因子载荷矩阵。因子载荷矩阵表示的是9个因子变量分别在3个公共因子上的相对重要性。

#计算主成分载荷矩阵 sqrt(a[i])*L[i][j]
for i in range(eigvector.shape[1]):
    for j in range(eigvector.shape[0]):
        eigvector[j][i]=np.sqrt(eig_value[i])*eigvector[j][i]
#以0.85为基准，求解特征根的个数
k=0
sa=0
for i in range(len(eig_value)):
    sa=sa+eig_value[i]
    k=k+1
    if sa/sum(eig_value)>=0.85:
        break
Keigvector=eigvector[0:,0:k]    #取前K个特征向量
a = pd.DataFrame(Keigvector)
tempcolname=[]
for i in range(k):
    tempcolname.append("factor" + str(i+1))
a.columns = tempcolname
a.index = ColName
print("\n因子载荷阵\n", np.round(a,3))

求得的因子载荷矩阵结果如下：

3.2.6 因子分析

求特殊因子的方差及公共因子解释的总方差（贡献率）。

#因子分析
fa = FactorAnalyzer(n_factors=k)
fa.loadings_ = a
print("\n公共因子方差:\n", fa.get_communalities())
var = fa.get_factor_variance()  # 给出贡献率
print("\n解释的总方差（即贡献率）:\n",var)

结果为：

由上可以看出，最后形成2个公共因子。从方差贡献率可以看出，其中第一个公因子解释了总体方差的67.54%，2个公共因子的方差贡献率为89.01%，可以较好的解释总体方差。

3.2.7 因子旋转

为了对提取的因子更好命名，且更好解释，我们需要进行因子旋转：

# 因子旋转
rotator = Rotator()
b = pd.DataFrame(rotator.fit_transform(fa.loadings_))
b.columns =tempcolname
b.index = ColName
print("\n因子旋转:\n", np.round(b,3))

因子旋转后的矩阵为：

我们可以看到，公共因子factor1与GDP、AW、APH、HSB的相关系数较大；公共因子factor2与NH、TRS、TP的相关系数较大。

再次进行因子分析：

#旋转后因子分析
fa = FactorAnalyzer(n_factors=k)
fa.loadings_ = b
print(fa.loadings_)
print("\n旋转后公共因子方差:\n", fa.get_communalities())
var = fa.get_factor_variance()  # 给出贡献率
print("\n旋转后解释的总方差（即贡献率）:\n",var)

其结果为：

3.2.8 旋转后的因子得分

计算旋转后的因子得分

# 因子得分
X1 = np.mat(DF_corr)
X1 = nlg.inv(X1)
b = np.mat(b)
factor_score = np.dot(X1, b)
factor_score = pd.DataFrame(factor_score)
factor_score.columns = tempcolname
factor_score.index = ColName
print("\n旋转后因子得分：\n", np.round(factor_score,3))

其结果为：

3.2.9 计算因子得分

首先，计算三个公共因子的权重：

#求各因子权重
AW=np.array(var)
AW=AW.T
TW=AW[:,1]
W=[]
for i in range(TW.shape[0]):
    W.append(TW[i]/sum(TW))
print("公共因子权重：\n",W)

其结果为：

之后计算各省份的因子得分：

#计算样本得分
fa_t_score = np.dot(np.mat(dataz), np.mat(factor_score))
print("样本的K个因子得分：\n",pd.DataFrame(fa_t_score))

计算结果为：

3.3 样本综合评分

根据各省份的三个因子得分，按照贡献率进行加权，可得得到各省份最后的综合得分，我们计算各省份的综合得分，并对其进行排名：

# 综合得分
wei = np.mat(W).T   #wei为权重向量
fa_t_score = np.dot(fa_t_score, wei) 
fa_t_score = pd.DataFrame(fa_t_score)
fa_t_score.columns = ['综合得分']
fa_t_score.insert(0,'ID',myID)
print("\n综合得分：\n", fa_t_score.sort_values(by='综合得分', ascending=False))

计算结果为：

最后的评分结果我们可以使用条形图来表示：

myname=np.array(fa_t_score["ID"])
mydata=np.array(fa_t_score["综合得分"])
draw_hist(myname,mydata)

依据我们所选择的指标进行评分，北京、上海、重庆、广州、深圳的常住居民的生活水平较高，位列前五名。从原始数据来看，这几个城市的平均工资、居民储蓄年末余额等，与其他市场相比都是比较高的，因此，该得分较为合理。另外，乌鲁木齐、呼和浩特、银川、海口、西宁五个城市的生活水平较低，排名靠后，这几个城市受限于地理位置、资源等因素，经济发展较为落后，总体评分较低。

四、实验总结

本文使用因子分析法来评估各个省会城市以及直辖市常住人口的生活水平，首先对7个因子进行分析，综合得到两个公共因子，计算公共因子得分及其权重，最后计算出各个城市的综合得分，依据得分高低对城市进行排序。选择的数据虽然适合用因子分析法来分析，但只选取了7个因子指标，指标量较少，最后得到的公共因子也只有两个。另外，评价一个城市居民的生活水平不仅要从经济方面来评价，还有考虑交通、环境等因素，而本文选取的多为经济方面的评价因子，所以最后的评分结果可能有失偏颇，希望以后在选取评价因子时能有所改进。

数据分析案例-基于因子分析探究各省份中心城市经济发展状况

一、实验背景

二、实验内容及数据

2.1 概述

2.2 变量介绍

三、实验步骤

3.1 导入模块和数据

3.2 数据分析

3.2.1 数据标准化

3.2.2 相关系数矩阵

3.2.3 KMO方法检测

3.2.4 特征值和特征向量

3.2.5 求因子载荷矩阵

3.2.6 因子分析

3.2.7 因子旋转

3.2.8 旋转后的因子得分

3.2.9 计算因子得分

3.3 样本综合评分

四、实验总结

热门文章

最新文章

相关课程

相关电子书

相关实验场景