前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >使用Stata完成广西碳酸钙企业的主成分分析和因子分析

使用Stata完成广西碳酸钙企业的主成分分析和因子分析

作者头像
润森
发布2020-03-12 18:14:03
1.8K0
发布2020-03-12 18:14:03
举报
文章被收录于专栏:毛利学Python

原文首发:https://maoli.blog.csdn.net/article/details/104787308

我们选取2018年的广西碳酸钙企业的数据,来判断那间企业在2018年更具有竞争力。我们来做主成分分析和因子分析。

下面是数据来源:

企业

净利润(万元)

营业总收入(万元)

期间费用(万元)

总资产周转率(次)

成本总额(万元)

流动资产(万元)

每股收益(元)

应收账款周转天数(天)

存货周转天数(天)

资产负债率(%)

八菱科技

721

71000

7303

0.29

72500

86100

0.0300

44.11

61.56

22.36

南宁化工股份有限公司

4488

27500

2714

0.48

29600

39000

0.2300

32.07

71.54

23.28

河池化工

-27400

23100

9244

0.38

51300

5321

-0.9311

10.92

35.35

164.52

柳州化工

37600

201000

30500

0.68

289000

114000

0.9000

12.88

49.82

20.83

想到主成分分析和因子分析,都会使用SPSS做。但是由于,我的SPPS上个月删掉了,占用1.5g内存,而且没有破解。这次,我用最不怎么熟悉的Stata来做主成分分析和因子分析。

主成分分析

在实际生活工作中,往往会出现所搜集的变量之间存在较强相关关系的情况。如果直接利用数据进行分析,不仅会使模型变得复杂,而且会带来多重线性的问题。主成分分析方法提供了解决这一问题的办法。

我们创建上面数据为2018年碳酸钙企业,通过Stata导入xlsx,注意:必须选择:将第一行作为变量名,不然你无法选择列名,一开始我以为列名不能有中文和括号,结果浪费我好多时间。

在这里插入图片描述

主成分在stata中的命令就是 pca ,其实了解sklearn就知道PCA(Principal Component Analysis),就是降维抽取维度。

用Stata真的不需要记住什么命令,明明可以找得到,不会一个help(pca)命令就欧克。

在这里插入图片描述

在这里插入图片描述

在Results界面中给出了分析结果

代码语言:javascript
复制
. pca 净利润万元 营业总收入万元 期间费用万元 总资产周转率次 成本总额万元 流动资产万元 每股收益
> 元 应收账款周转天数天 存货周转天数天 资产负债率

Principal components/correlation                 Number of obs    =          4
                                                 Number of comp.  =          3
                                                 Trace            =         10
    Rotation: (unrotated = principal)            Rho              =     1.0000

# 第一个表格给出了相关系数矩阵的特征值和比例,结果界面和SPSS类似,只是多了一个Difference列

    --------------------------------------------------------------------------
       Component |   Eigenvalue   Difference         Proportion   Cumulative
    -------------+------------------------------------------------------------
           Comp1 |      6.27881      3.19118             0.6279       0.6279
           Comp2 |      3.08763      2.45408             0.3088       0.9366
           Comp3 |      .633554      .633554             0.0634       1.0000
           Comp4 |            0            0             0.0000       1.0000
           Comp5 |            0            0             0.0000       1.0000
           Comp6 |            0            0             0.0000       1.0000
           Comp7 |            0            0             0.0000       1.0000
           Comp8 |            0            0             0.0000       1.0000
           Comp9 |            0            0             0.0000       1.0000
          Comp10 |            0            .             0.0000       1.0000
    --------------------------------------------------------------------------


# STATA给出了所有特征值的特征向量,各列的数值平方之和等于1,不信你求Comp1的平方和
# 系数越大,说明主成份对该变量的代表性越大。SPSS比Stata没有 Unexplained 这一列
Principal components (eigenvectors)
    ----------------------------------------------------------
        Variable |    Comp1     Comp2     Comp3 | Unexplained
    -------------+------------------------------+-------------
      净利润万元 |   0.3900    0.1078   -0.1188 |           0
     营业总收~元 |   0.3834   -0.1160    0.2365 |           0
    期间费用万元 |   0.3417   -0.2835    0.1727 |           0
     总资产周~次 |   0.3258   -0.1813   -0.6052 |           0
    成本总额万元 |   0.3699   -0.2012    0.1573 |           0
    流动资产万元 |   0.3564    0.1525    0.4544 |           0
      每股收益元 |   0.3698    0.1969   -0.1850 |           0
     应收账款~天 |  -0.0694    0.5357    0.3637 |           0
     存货周转~天 |   0.0564    0.5384   -0.3662 |           0
      资产负债率 |  -0.2634   -0.4271    0.0454 |           0
    ----------------------------------------------------------

SPSS比Stata没有 Unexplained 这一列。我们通过predict pc1 pc2 pc3 score命令就可以查看SPSS的界面如何

代码语言:javascript
复制
. predict pc1 pc2 pc3 score
(score assumed)
(extra variables dropped)

Scoring coefficients
    sum of squares(column-loading) = 1

    --------------------------------------------
        Variable |    Comp1     Comp2     Comp3
    -------------+------------------------------
      净利润万元 |   0.3900    0.1078   -0.1188
     营业总收~元 |   0.3834   -0.1160    0.2365
    期间费用万元 |   0.3417   -0.2835    0.1727
     总资产周~次 |   0.3258   -0.1813   -0.6052
    成本总额万元 |   0.3699   -0.2012    0.1573
    流动资产万元 |   0.3564    0.1525    0.4544
      每股收益元 |   0.3698    0.1969   -0.1850
     应收账款~天 |  -0.0694    0.5357    0.3637
     存货周转~天 |   0.0564    0.5384   -0.3662
      资产负债率 |  -0.2634   -0.4271    0.0454
    --------------------------------------------

下面我们来查看命令corr pc1 pc2 pc3 净利润万元 营业总收入万元 期间费用万元 总资产周转率次 成本总额万元 流动资产万元 每股收益元 应收账款周转天数天 存货周转天数天 资产负债率相关性矩阵

在这里插入图片描述

我们如何去除呢?加个comp(2) covariance

  • com(#)指定最多只能有#个主成分
  • covariance利用协方差矩阵,默认利用相关系数矩阵

命令pca 净利润万元 营业总收入万元 期间费用万元 总资产周转率次 成本总额万元 流动资产万元 每股收益元 应收账款周转天数天 存货周转天数天 资产负债率, comp(2) covariance

在这里插入图片描述

我们用命令 screeplot作为碎石图,英语不好的我连eigenvalues什么意思都不知道,查了下原来是特征值

我们在使用命令loadingplot画载荷图,选择出最具有成分的两个成分的作为相关图,我们从相关图就完全看出是什么元素决定成分了。

在这里插入图片描述

因子分析

下面我们做因子分析,做前,我先吹下什么是因子分析:

因子分析(factor analysis)是用少数的不可观察的潜变量表示多数可观察的相关的变量 。

就是今天你突然间要上大街要饭,我要找出原因,你为什么今天跑去大街要饭,还用想吗,还不是因为你没钱没房

金钱因子。

在Python中主要使用from sklearn.decomposition import FactorAnalysispip install factor_analyzer

今天用Stata,点点就ok了 ,敲什么代码,再说了我又不会敲。

在这里插入图片描述

在这里插入图片描述

命令:factor 净利润万元 营业总收入万元 期间费用万元 总资产周转率次 成本总额万元 流动资产万元 每股收益元 应收账款周转天数天 存货周转天数天 资产负债率

在这里插入图片描述

因子1是利润营收决定的,说到底就是钱,因子2就是你的碳酸钙什么时候可以周转出去,最快时间把钱要到手里。负债的企业还有什么竞争力,赶紧关门跑路。

一般因子分析,就需要预测,命令很简单predict f1

在这里插入图片描述

排名第一的就是净利润,企业就是要赚钱才有竞争力。这不是常识了,分析和没分析差不多,我还是干Spring,Django吧。

看看上面内容,没有代码怎么行?虽然是化工专业,照样敲码两不误。虽说我专业都快挂完了,心痛一秒钟。

下面是百度百科给因子分析模型定义的,我抄了下。

的一种统计方法,是一种降维技术. 做因子分析的前提是自变量之间有相关关系. 这里的潜变量就是我们所求的因子,自变量是因子 的表征。

因子分析模型是把原观测变量分解成公共因子和特殊因子两部分

其中是原始变量标准化后的数据,是公共因子,是特殊因子。

因子分析的一般步骤

  • 将原始数据标准化处理
  • 计算相关矩阵
  • 计算相关矩阵的特征值和特征向量
  • 确定公共因子个数
  • 构造初始因子载荷矩阵 ,其中为的特征向量
  • 建立因子模型
  • 对初始因子载荷矩阵A进行旋转变换,旋转变换是使初始因子载荷矩阵结构简化,关系明确,使得因子变量更具有可解释性,如果初始因子不相关,可以用方差极大正交旋转,如果初始因子间相关,可以用斜交旋转,进过旋转后得到比较理想的新的因子载荷矩阵.
  • 将因子表示成变量的线性组合,其中的系数可以通过最小二乘法得到.
  • 计算因子得分

看看一般步骤,读取数据我就pass了

在这里插入图片描述

在这里插入图片描述

在这里插入图片描述

初始因子和Stata的结果一样

在这里插入图片描述

在Stata中我们没有旋转变换,

在这里插入图片描述

旋转变换的后的

答案是柳州化工,我听说柳州螺蛳粉,五菱。原来柳州的碳酸钙,用来干嘛,建房 遇到盐酸放出,是重要的碱性产品。

在这里插入图片描述

对了,如果想取出没用因子,只需要在Stata命令参数添加 pcf,

在这里插入图片描述

旋转更简单 rotate,搞定,如果报错 ssc install rotate安装下

在这里插入图片描述

Python代码如何设置,将确定公共因子个数解释度>=1 改为 >=0.8,这个就去除了没用因子

下面完整代码

代码语言:javascript
复制
import pandas as pd
import numpy as np
datafile = '2018年碳酸钙.xlsx'
data_2018 = pd.read_excel(datafile,index_col=0)
data_2018.head()
data_2018_mat=(data_2018-data_2018.mean())/data_2018.std()# 0均值规范化
data_2018_mat.corr()

import numpy.linalg as nlg #导入nlg函数 矩阵求逆
eig_value,eig_vector=nlg.eig(data_2018_mat.corr()) #计算特征值和特征向量
eig=pd.DataFrame() #利用变量名和特征值建立一个数据框
eig['names']=data_2018.columns#列名
eig['eig_value']= eig_value.real.astype('float64')#特征值 (复数的j都是0) 用real取出

for k in range(1,9): #确定公共因子个数
    if eig['eig_value'][:k].sum()/eig['eig_value'].sum()>=0.8: #如果解释度达到80%, 结束循环
        print(k)
        break
        
from math import sqrt
col0=list(sqrt(eig_value[0].real)*eig_vector[:,0].real) #因子载荷矩阵第1列
col1=list(sqrt(eig_value[1].real)*eig_vector[:,1].real) #因子载荷矩阵第2列

col0 =[-l for l in col0]
col1 =[-l for l in col1]
A2018=pd.DataFrame([col0,col1]).T #构造因子载荷矩阵A
A2018.columns=['factor1','factor2'] #因子载荷矩阵A的公共因子
A2018.index = data_2018_mat.corr().index

A2018

h2018=np.zeros(10) #变量共同度,反映变量对共同因子的依赖程度,越接近1,说明公共因子解释程度越高,因子分析效果越好
D2018=np.mat(np.eye(10))#特殊因子方差,因子的方差贡献度 ,反映公共因子对变量的贡献,衡量公共因子的相对重要性
A2018=np.mat(A2018) #将因子载荷阵A矩阵化
for i in range(10):
    a2018=A2018[i,:]*A2018[i,:].T #A的元的行平方和
    h2018[i]=a2018[0,0]  #计算变量X共同度,描述全部公共因子F对变量X_i的总方差所做的贡献,及变量X_i方差中能够被全体因子解释的部分
    D2018[i,i]=1-a2018[0,0] #因为自变量矩阵已经标准化后的方差为1,即Var(X_i)=第i个共同度h_i + 第i个特殊因子方差
    
from numpy import eye, asarray, dot, sum, diag #导入eye,asarray,dot,sum,diag 函数
from numpy.linalg import svd #导入奇异值分解函数
def varimax(Phi, gamma = 1.0, q =20, tol = 1e-6): #定义方差最大旋转函数
    p,k = Phi.shape #给出矩阵Phi的总行数,总列数
    R = eye(k) #给定一个k*k的单位矩阵
    d=0
    for i in range(q):
        d_old = d
        Lambda = dot(Phi, R)#矩阵乘法
        u,s,vh = svd(dot(Phi.T,asarray(Lambda)**3 - (gamma/p) * dot(Lambda, diag(diag(dot(Lambda.T,Lambda)))))) #奇异值分解svd
        R = dot(u,vh)#构造正交矩阵R
        d = sum(s)#奇异值求和
        if d_old!=0 and d/d_old < 1 + tol: break
    return dot(Phi, R)#返回旋转矩阵Phi*R
rotation_mat=varimax(A2018)#调用方差最大旋转函数
rotation_mat=pd.DataFrame(rotation_mat,columns=['factor1','factor2'],index=data_2018_mat.corr().index)#数据框化

rotation_mat

factor_score_2018=(data_2018_mat).dot(rotation_mat) #计算因子得分
factor_score_2018=pd.DataFrame(factor_score_2018)#数据框
factor_score_2018.columns=['利润因子','拿钱速度因子'] #对因子变量进行命名
factor_score_2018.index =  data_2018.index
factor_score_2018

在除去没用因子,旋转的变换和Stata完全一样。

在这里插入图片描述

这是计算各因子的得分

在这里插入图片描述

各因子的计算原理

在这里插入图片描述

然而Stata计算总因子得分没有命令,计算公式:因子得分*因子方差的贡献率/累计方差贡献率作为权重。然后计算

方差百分比

八菱科技总因子得分:(0.6131* (-1.637494) + 0.3236 * 2.396507 )/ 0.9367

柳州化工总因子得分: (0.6131* (8.935557) + 0.3236 * 0.262321 )/ 0.9367

虽然,数据很少,而且我们我很容易看出柳州化工是行业的大佬,

Stata相对于SPSS,就是没有那个各因子得分,计算总因子得分变得非常困难,因此计算的时候,反而想下载SPSS,但是SPSS免费试用14天已过。又不知道哪里下载盗版的,反而使用Python从原理计算出因子得分。

使用SPSS比Stata更适合主成分分析和因子分析,但是Stata是一款医学研究的软件,提供了大量的统计分析

相对的SPSS的更全,比如生存,时间序列,甚至有时连Python深度模型跑出来的,还不如用Stata点一点,Stata虽然命令多,但是完全不需要记忆,在窗口中完全可以找到,或者一个 help(命令)查看示例。

而SPSS两款工具,SPSS Modeler和SPSS Statistics是SPSS中的“哼哈二将”,一个负责统计分析,一个负责挖掘。

还有不要老是敲代码,有时候工具点几点就ok,excel也是,不要老是跑动不动就pandas读取文件,然后一无所知,比如做个时间序列,keras跑深度学习写几十行代码以为自己很牛,却不知道有的人使用Stata点一点就完成。我们在完成需求的同时,应该寻找最佳的方案,比如一个excel就可以完成的,干嘛要打开jupyter notebook?我们还是干回爬虫,Spring和Djiango吧

本文参与 腾讯云自媒体同步曝光计划,分享自微信公众号。
原始发表:2020-03-11,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 小刘IT教程 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体同步曝光计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
目录
  • 主成分分析
  • 因子分析
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档