R语言因子分析

安装加载包

library(mvstats)
library(factoextra)
library(ggplot2)

准备数据

> mydata=read.csv("E:/R/data.csv")#读取车险数据，命名成mydata

任务一：描述统计

> summary(mydata)

    行驶里程         行驶时长          平均时速       时速标准差    
 Min.   :  0.04   Min.   :0.05056   Min.   : 0.70   Min.   : 2.088  
 1st Qu.:  1.00   1st Qu.:0.15326   1st Qu.:14.03   1st Qu.:11.987  
 Median :  3.00   Median :0.29264   Median :20.40   Median :16.496  
 Mean   : 12.17   Mean   :0.48330   Mean   :24.01   Mean   :18.058  
 3rd Qu.: 14.00   3rd Qu.:0.58257   3rd Qu.:31.83   3rd Qu.:23.438  
 Max.   :301.00   Max.   :3.66944   Max.   :90.22   Max.   :47.444  
 
   时速最大值         平稳性          转速平均值       转速标准差    
 Min.   :  6.00   Min.   :0.02432   Min.   : 746.4   Min.   : 84.16  
 1st Qu.: 43.00   1st Qu.:0.94901   1st Qu.:1057.5   1st Qu.:285.19  
 Median : 59.50   Median :1.05037   Median :1183.3   Median :367.32  
 Mean   : 63.96   Mean   :1.02643   Mean   :1224.7   Mean   :375.59  
 3rd Qu.: 82.00   3rd Qu.:1.13851   3rd Qu.:1356.6   3rd Qu.:457.03  
 Max.   :177.00   Max.   :1.47343   Max.   :2462.3   Max.   :890.78  
 
   转速最大值      早晚高峰         深夜出行          极度拥堵      
 Min.   : 982   Min.   :0.0000   Min.   :0.00000   Min.   :0.02439  
 1st Qu.:1790   1st Qu.:0.0000   1st Qu.:0.00000   1st Qu.:0.26104  
 Median :2117   Median :0.0000   Median :0.00000   Median :0.37893  
 Mean   :2176   Mean   :0.3451   Mean   :0.05001   Mean   :0.40231  
 3rd Qu.:2441   3rd Qu.:1.0000   3rd Qu.:0.00000   3rd Qu.:0.53182  
 Max.   :5411   Max.   :1.0000   Max.   :1.00000   Max.   :1.00000  
 
    高速行驶          疲劳驾驶            车号      
 Min.   :0.00000   Min.   :0.00000   Min.   :1.000  
 1st Qu.:0.00000   1st Qu.:0.00000   1st Qu.:3.000  
 Median :0.00000   Median :0.00000   Median :5.000  
 Mean   :0.01709   Mean   :0.09008   Mean   :4.676  
 3rd Qu.:0.00000   3rd Qu.:0.00000   3rd Qu.:6.000  
 Max.   :0.61943   Max.   :3.66944   Max.   :8.000  

由数据信息得出，变量包含：
    时间的维度有：行驶时长、疲劳驾驶、早晚高峰、深夜出现等4个变量；
    速度维度的有：平均时速、时速最大值、时速标准差、极度拥堵、高速行驶等5个变量；
    里程维度有：行驶里程这1个变量；
    机械性能的维度有转速平均值、转速最大值、转速标准差等3个变量；
    驾驶的平稳性指标变量为平稳性1个变量。
    其中车辆行程最大行驶里程为301千米，行驶时长最长为3.7小时，平均时速最大达到90.22千米每小时，时速标准的均值为18.05小时，时速最大值为177千米每小时，其中平稳性最大为1.47，转速的平均值的平均水平为1224.7千米每转，转速标准差的最大值最890.78，转速最大值的平均值为2176转，早晚高峰的时间平均在0.3451，深夜出行的时段平均在0.05，极度拥堵的平均水平为0.40231，高速行驶的平均水平为0.01709，疲劳驾驶的平均水平为0.09008，车号的水平共有8个，共8辆车。

任务二：对数行驶里程和平均时速的箱线图

对数行驶里程=log(mydata$行驶里程) #对行驶里程取对数
> par(mfrow=c(1,2))#画1*2图
> boxplot(对数行驶里程~mydata$车号,xlab="车号",ylab="对数行驶里程",col=c("lightgrey"))#画对数行驶里程箱线图
> boxplot(mydata$平均时速~mydata$车号,xlab="车号",ylab="平均时速(千米/每小时)",col=c("lightgrey"))#平均时速箱线图

从箱线图中可以看出，一共有8辆车，图中研究不同的车号和对数行驶里程、平均时速间的关系。其中对数行驶里程、平均时速和车号做的箱线图可以看出均有离开箱子外的点，说明数据仍存在不规整，需要进一步的数据处理和验证。
    从上图中左边的箱线图可以看出：
    （1）不同的车号行驶里程是不同的。
    （2）2号车、3号车、4号车的对数行驶里程平均水平远大于其它车号，可将这三个车号得行驶里程归为长途行驶类车辆，而其余车：1号、5号、6号、7号、8号对数行驶里程平均水平较长途行驶类车辆平均水平小，可归为短途类行驶车辆；
    （3）其中3号车的对数行驶里程在平均水平上远高于其它车辆，推测3号车可能是长途快运城；
    （4）而7号类车辆的对数行驶里程的平均水平远低于其它车辆，推测7号车可能是私家的短途车。
    （5）发现5号车的箱子长度最大，看出5号车的样本量最大，行驶里程的记录数据比其它车辆的数据更多，说明5号车在数据样本量上更加充足。
    从上图中右边的箱线图可以看出：
    （1）不同的车号产生的平均时速是会发生变化的。
    （2）其中3号车的平均时速的平均水平最大，而对应3号车的对数行驶里程平均水平也是最大的，推测3号车在相同里程下比其它车行驶花费的时间少；
    （3）1号车的平均时速的平均水平较其它车小，推测很可能在相同时间内，1号车行驶的里程可能较其它车的行驶距离小；
    （4）2号车、3号车、4号车的平均时速的平均水平远大于其它车号，可将这三个车号的平均时速归为快速类车辆，
    （5）而其余车：1号、5号、6号、7号、8号平均时速的平均水平较快速类车辆平均水平小，推测这些车可归为稳定行驶类行驶车辆。

任务三：因子分析

因子分析是关于降维简化数据的技术,往往因子的选取标准的方差贡献率大于80%，或者特征根（特征值）大于1，或者用碎石图的“”趋于平衡点点来判断。为此下面对公因子深入进行分析：

> mydata1=mydata[,-15] #删除第15列
> FA4=factpc(mydata1,4,rotation="varimax")#极大似然法进行因子旋转并且保留四个因子
 Factor Analysis for Princomp in Varimax: 
> FA4$Vars#旋转后方差贡献率
         Vars Vars.Prop Vars.Cum
Factor1 3.680    26.286    26.29
Factor2 3.301    23.582    49.87
Factor3 2.474    17.671    67.54
Factor4 1.236     8.828    76.37

1）因子分析是关于降维简化数据的技术。
    （2）去掉最后一列车号，对剩余变量做因子分析。保留4个因子。从上看出Vars是方差贡献、Vars.Prop是方差贡献率、Vars.Cum是累计方差贡献率；其中方差贡献和方差贡献率的值越大，说明相应因子变量越重要，是衡量因子变量重要性的重要指标。累计方差贡献率是因子方差贡献和得累计，因子个数设定的越多，累计方差贡献率越大。
    其中第一个因子的方差贡献是3.68，因子的波动是在3.68，方差贡献率是26.286%，贡献率的波动情况是26.286%.累计方差贡献率是26.29%。
    第二个因子的方差贡献是3.301，方差贡献率是23.582，累计方差贡献率是49.87，此时两个因子不足以解释整个变量；
    第三个因子的方差贡献是2.474，方差贡献率是17.671，累计方差贡献率是67.54,此时三个因子的解释所有变量的能力进一步提升；
    第四个因子的方差贡献是1.236，方差贡献率是8.282，累计方差贡献率是76.37，此时累计方差贡献率可以解释变量的能力达到76.37%，以不小于80%为累计方差贡献率解释的最佳值，达到76.37%的累计贡献率的解释效果说明较好，聚成4个因子用来解释变量的效果还是较好的，为此进一步对因子进行分析。

> z=round(x=FA4$loadings,digits=3)#旋转因子载荷矩阵
> 变量共同度=vector(length=14)
> for(i in 1:14){#计算因子变量个数为4时，原有变量的变量共同度
+ 变量共同度[i]=round(sum(z[i,1:4]^2),3)#因子载荷矩阵的平方和并保留3位数
+ }
> cbind(z[,1:4],变量共同度)#合并因子载荷矩阵和变量共同度
           Factor1 Factor2 Factor3 Factor4 变量共同度
行驶里程     0.394   0.829   0.181   0.048      0.878
行驶时长     0.192   0.847   0.289   0.041      0.839
平均时速     0.873   0.362   0.226   0.118      0.958
时速标准差   0.677   0.392   0.482   0.063      0.848
时速最大值   0.653   0.420   0.504   0.047      0.859
平稳性       0.473   0.263   0.201  -0.097      0.343
转速平均值   0.711   0.178   0.509   0.185      0.831
转速标准差   0.273   0.153   0.864   0.092      0.853
转速最大值   0.240   0.215   0.858   0.077      0.846
早晚高峰    -0.014  -0.092  -0.045  -0.761      0.590
深夜出行     0.121  -0.041   0.084   0.746      0.580
极度拥堵    -0.886   0.008  -0.094  -0.132      0.811
高速行驶     0.424   0.712   0.037  -0.011      0.688
疲劳驾驶    -0.015   0.865   0.138   0.038      0.769

从因子运行的结果前三列数值知，车行驶过程中涉及指标的14个变量和4个因子的相关系数，从相关系数值可以得出数值的绝对值越大，因子与变量的相关密切程度越高，也反映出该变量对公共因子的相对重要性。变量共同度描述出全部公共因子对变量总方差作出的总贡献，也成为个性方差。
    从因子的相关性的程度知：行驶里程与第二个因子的相关度为0.829远大于其它因子的相关度，行驶时长可归为第二个因子；行驶时长与第二个因子的相关度为0.847远大于其它因子的相关度，可归为第二个因子中；同理可得，对比所有变量与因子的相关度的大小，可将所有因子对应进行划分，划分结果为：
    （1）与第一个因子相关度最高的变量有：平均时速、时速标准差、极速拥堵、时速最大值、转速平均值等5个变量，同时说明Factor1包含的变量是这5个变量.
    （2）与第二个因子相关度最高的变量有：行驶里程、行驶时长、高速行驶、疲劳驾驶等4个变量。同时说明Factor2包含的变量是这4个变量。
    （3）与第三个因子相关度高的变量有：转速标准差、转速最大值等2个变量。同时说明Factor3包含的变量是这2个变量。
    （4）与第四个因子的相关度高的变量有：深夜出行、早晚高峰等2个变量。同时说明Factor4包含的变量是这2个变量.
    （5）其中平稳性变量与四个因子的相关系数最大值为0.473，而相关系数值在0.3到0.5属于低度相关，也就说明平稳性变量越因子的紧密程度较弱。与四个因子对比中，发现平稳性与第一个因子的相关系数值为0.473，大于和其它变量的相关度，如果划分平稳性到因子中可以划分到第一个因子中。
    通过因子的紧密性变量的联系进行归纳，对四个因子分别命名：Factor1命名为车速状况因子、Factor2命名为行驶状态因子、Factor3命名为车速波动因子、Factor4命名为行驶习惯因子，共4个水平因子。

任务四：选择保留合适公因子

分别尝试2个、3个、4个、5个公因子对车行驶的驾驶中的变量进行降维，而因子分析是对变量进行提炼，因子个数往往通过一定规则进行确定，而且因子的形式不是唯一确定，因子分析又与主成分分析有许多相似之处，因子分析的本质就是减少变量，达到降维的结果。下面对因子保留的个数从以下三种方法综合分析进行确定：

方法一：碎石图确定因子个数

> par(mfrow=c(1,1))
> screeplot(princomp(mydata1),type="lines")#用princomp函数（使用特征分解办法求解主成分),绘制碎石图，判断出保留3个因子

从图中可以看出，第一个特征值较大，第二个特征值次之，减小幅度很大，形成一个“陡峭的山坡”，第四个特征值都很小，从第4个的特征值都很小（根据第3个和第4个特征值间连线粗糙痕迹判断第四个特征值小），图形趋于平缓，类似“高山脚下的碎石”，因此选择“山脚”下特征值的编号应为因子变量的个数，选择4个公因子作为变量的降维结果。

方法二：对因子个数的累计方差贡献率逐一计算观察、分别尝试

尝试2个、3个、4个、5个公因子进行累计方差贡献率进行计算观察：

尝试2个公因子**

> FA2=factpc(mydata1,2,rotation="varimax")
 Factor Analysis for Princomp in Varimax: 
> FA2$Vars
         Vars Vars.Prop Vars.Cum
Factor1 4.315     30.82    30.82
Factor2 4.293     30.66    61.48

从上看出：发现2个因子的累计方差贡献率为61.48%，说明可以解释变量的水平为61.48%。

尝试3个公因子

 > FA3=factpc(mydata1,3,rotation="varimax")
 Factor Analysis for Princomp in Varimax: 
> FA3$Vars
         Vars Vars.Prop Vars.Cum
Factor1 5.067    36.190    36.19
Factor2 3.372    24.085    60.28
Factor3 1.260     9.001    69.28

从上看出：3个因子累计方差贡献率为69.28%，说明可以解释变量的水平为69.28%。

尝试4个公因子

> FA4=factpc(mydata1,4,rotation="varimax")
Factor Analysis for Princomp in Varimax: 
> FA4$Vars
         Vars Vars.Prop Vars.Cum
Factor1 3.680    26.286    26.29
Factor2 3.301    23.582    49.87
Factor3 2.474    17.671    67.54
Factor4 1.236     8.828    76.37

从上看出：4个因子的方差贡献率为76.37%，说明可以解释变量的水平为76.37%。

尝试5个公因子

> FA5=factpc(mydata1,5,rotation="varimax")
Factor Analysis for Princomp in Varimax: 
> FA5$Vars
         Vars Vars.Prop Vars.Cum
Factor1 3.489    24.925    24.92
Factor2 3.346    23.899    48.82
Factor3 2.620    18.711    67.53
Factor4 1.089     7.779    75.31
Factor5 1.026     7.328    82.64

从上看出：5个因子的方差贡献率为82.64%，说明可以解释变量的水平为82.64%。
    通常认为：当累计方差贡献率大于80%时的因子个数即为因子变量的选取个数，而因子分析的目的是减少原有变量的个数，为了达到降维的效果。如果选择4个因子相比较80%水平下，此时76.73%非常接近于80%，与80%求差后只相差3.4%，变量从14个降维变为4个因子变量；同理，如果选择5个因子相比较1，此时82.64%超过80%，与80%求差只增加2.64%，变量从14个降维为5个因子变量。围绕浓缩原有变量的本质为出发点，而且在不同的行业中要求方差贡献率的水准也不相同，有的60%即可，有的70%水平即可，为此从因子分析的降维思想角度考虑，选择4个因子作为最终的因子选择个数。

方法三：利用特征值和累计方差贡献率确定

> RMatrix=cor(mydata1)#计算相关系数矩阵
> Result=eigen(RMatrix)#求相关系数矩阵的特征值和特征向量
> A=NULL
> for(i in 1:14){   #因子载荷矩阵，因子个数等于原有变量个数
+  k=sqrt(Result$values[i])*Result$vectors[,i]# 计算平方根
+  A=cbind(A,k)#合并数据
+ }
> colnames(A)=paste("k",1:14,sep="") #重命名
> Result$values
 [1] 6.94933 1.65805 1.09122 0.99274 0.87847 0.75452 0.51511 0.41693 0.27876
[10] 0.19185 0.11581 0.09828 0.03806 0.02087
> cumsum(Result$values/sum(Result$values))
 [1] 0.4964 0.6148 0.6928 0.7637 0.8264 0.8803 0.9171 0.9469 0.9668 0.9805
[11] 0.9888 0.9958 0.9985 1.0000

确定因子变量个数往往选取的标准为：一般选取大于1的特征值，根据因子变量的方差贡献判断因子的变量重要性，若因子变量是不应略去的重要因子，则它至少应该能够解释14个原有变量中的1个；通常，累计方差贡献率大于0.8时的因子个数即为因子变量确定的最终个数。 从特征值可以看出：大于1的特征值因子为前3个，第4个因子的特征值为0.99274非常接近于1；第5个因子的方差贡献率为82.64%，而第四个因子的方差贡献率为76.37%，非常接近于80%。为此从特征值和方差贡献率看，因子的选取个数可能为3个、4个、5个，取一个平均水平，在此选择因子的个数为4个公因子。

综上，通过3种方法分析，保留4个公因子作为变量的最终的降维结果。

任务五：聚类分析

> wss=(nrow(mydata1)-1)*sum(apply(mydata1,2,var))
> for (i in 1:14)
+ wss[i]=sum(kmeans(mydata1,centers=i)$withinss)
> plot(1:14, wss, type="b", xlab="Number of Clusters",ylab="Within groups sum of squares")#每个聚类内部的距离平方和

通过任务三的因子的得分，对聚类指标做聚类碎石图，碎石图中的横轴为各类的距离，纵轴为类数目，随着类的不断凝聚，类数目的不断减少。类间距离在不断增大。从第4类开始。形成极为“平坦的碎石路。”根据类间距离小形成类的相似性大，类间距离大形成类的相似性小的原则，可以找到“山脚”3作为确定分类数目的参考，从第5类起，“山脚的路”趋于平稳，为此选择4类以作为确定分类数目的最终聚类标准。

按因子得分聚类

> FA4=factpc(mydata1,4,rotation="varimax")
 Factor Analysis for Princomp in Varimax: 
> fa4.score=FA4$scores#获得每段路程的因子得分，每行是一个行程，每列是该行程在相应因子上的得分
> colnames(fa4.score)=c("车速状况","行驶状态","车速波动","行驶习惯")
> set.seed(1234)#设定随机种子
> clur=kmeans(fa4.score,4,nstart=20)#按4类进行聚类
> clur$centers#查看4类的类中心
  车速状况 行驶状态 车速波动 行驶习惯
1   0.4826 -0.26480   0.6219 -0.32622
2  -0.5256 -0.07889  -0.6380 -0.03331
3   0.5727 -0.47528   0.3313  3.14438
4   0.3812  4.12063   0.5179  0.04307

类中心解读：第一类的中心向量为（0.4826，-0.26480，0.6219，-0.32622）；第二类的中心向量为（-0.5256，-0.07889，-0.6380，-0.03331,）；第三类的中心向量为（0.5727，-0.47528,0.3313,3.1448）；第四类中心向量为（0.3812,4.12063,0.5179,0.04307）

> clur$size#查看每类的样本量
[1] 324 361  37  32

类样本量情况：将数据聚类为4类，第一类样本量为324个，第二类样本量为361个，第三类样本量为37个，第四类样本量为32个。

每一类行程代表的含义和驾驶习惯： 第一类的驾驶习惯：该类车的转速标准差和转速最大值较大，可能是爱加油的车、司机猛踩油门，转速较高，常在空旷路段上行驶加速。 第二类的驾驶习惯：车速较快，平均时速快，司机的技巧娴熟，速度的值较为稳定，经常遇到红绿灯等拥堵情况，开车的效果很稳定，可能该车是属于室内的出租车内型。 第三类的驾驶习惯：行驶速度和行驶强度较大，开车的时间较长，夜间行驶，长途行驶开车的性能效果较好，司机的水平较好的老司机开车。 第四类的驾驶习惯：开车的行驶频率强度大，时间长，且开车效果较稳。

反思任务三中是需要去掉最后一列车号，对剩余的变量做因子分析。如果不去掉最后一列车号，聚类的结果又是如何的呢？同样分为四类，聚类的结果也分成4类，与驾驶习惯进行对应得出每辆车的聚类结果，一一对应后，按相同的标准可以进一步得出每类车的驾驶习惯。能得出每辆车的出行驾驶习惯，推测进而可以说明该车号司机的驾驶习惯。

任务六：描述驾驶行为

好的司机开车往往水平技术好，而评定驾驶习惯和驾驶行为需要从多方面来考虑，为此从以下两点进行分析：

驾驶习惯分析

> temp=table(mydata$车号,clur$cluster) #查看8位驾驶人行为分布
> temp.scale=apply(temp,1,function(x) x/sum(x)) #转化成比例数据，使结果更直观
> colors=c("lightgreen","lightblue","lightgrey","orange")
> classification=c("类别1","类别2","类别3","类别4")
> barplot(temp.scale,main = "",xlab = "比例",ylab = "车号",col = colors,horiz=T,
+ legend.text=paste("类",c(1:4),sep=""),args.legend=list(x="top",bty="n",inset=-0.24))  

从图中可以得出，第一类的车有：4号车、2号车、5号车；第二类的车有：6号车、7号车、8号车、3号车，第三类的车有：1号车。对每一类的车与上面的车进行对比，可以看出4号车和2号车的转速标准差和转速最大值较大，可能是爱加油的车、司机猛踩油门，转速较高，常在空旷路段上行驶加速；6号车、7号车、8号车、3号车的车速较快，平均时速快，司机的技巧娴熟，速度的值较为稳定，经常遇到红绿灯等拥堵情况，开车的效果很稳定，可能该车是属于室内的出租车内型；其中1号车行驶速度和行驶强度较大，开车的时间较长，夜间行驶，长途行驶开车的性能效果较好，司机的水平较好的老司机开车。

> boxplot(mydata$平稳性~mydata$车号,xlab="车号",ylab="平稳性",col=c("green"))

从箱线图可以看出，5号车的平稳性值最高，8号车的平稳性最差。推测出5号车的司机技术熟练度可能比8号车的司机技术好。那么区分好车的平稳性后，晕车的人或者寻求刺激速度快的乘客可以进行有效的选择车的行驶状况，从而能达到满意的就乘车辆。