加国头条 i4cc.ca ｜基于“洛杉矶街区”数据的因子分析和主成分分析

一.数据描述

读取数据，并观察数据结构，可以看到共有110个观测值，15个变量。

> data<-read.csv("LA.Neighborhoods.csv")
> str(data)~q~data.frame~q~:   110 obs. of  15 variables:
 $ LA.Nbhd   : Factor w/ 110 levels "Adams_Normandie",..: 1 2 3 4 5 6 7 8 9 10 ...
 $ Income    : int  29606 65649 31423 53872 37948 208861 168104 63039 105253 33235 ...
 $ Schools   : int  691 719 687 762 656 924 0 791 872 689 ...
 $ Diversity : num  0.6 0.4 0.8 0.9 0.4 0.2 0.1 0.2 0.2 0.1 ...
 $ Age       : int  26 29 31 34 36 46 45 38 39 25 ...
 $ Homes     : num  0.26 0.29 0.31 0.34 0.36 0.46 0.45 0.38 0.39 0.25 ...
 $ Vets      : num  0.05 0.07 0.05 0.06 0.1 0.13 0.1 0.05 0.08 0.03 ...
 $ Asian     : num  0.05 0.11 0.13 0.2 0.05 0.08 0.04 0.05 0.07 0.02 ...
 $ Black     : num  0.25 0.02 0.25 0.01 0.71 0.01 0.02 0.02 0.04 0.01 ...
 $ Latino    : num  0.62 0.72 0.57 0.51 0.17 0.05 0.03 0.06 0.06 0.94 ...
 $ White     : num  0.06 0.13 0.05 0.22 0.03 0.83 0.88 0.82 0.8 0.02 ...
 $ Population: int  31068 31068 22106 14888 30123 7928 10610 21417 6080 92785 ...
 $ Area      : num  0.8 3.1 1 1.8 3 6.6 7.9 1.7 0.8 6.5 ...
 $ Longitude : num  -118 -118 -118 -118 -118 ...
 $ Latitude  : num  34 34.2 34 34.1 34 ...

各个变量的含义如下，包含了各个街区中收入、种族、人口、面积、位置等信息，我们进一步想要分析出各个街区在这些不同方面的表现，由于变量数目较多，可以采用降维的方法。

LA.Nbhd:街区名字
Income：收入中位数
Schools：公立学校API中位数
Diversity：种族多样性（0-1）打分
Age：年龄中位数
Homes：有房家庭比例
Vets：复员军人比例 
Asian：亚裔人口比例
Black：非裔比例 
Latino：拉美裔比例
White：欧裔比例
Population：人口
Area：面积 
Longitude：经度
Latitude：纬度

二.数据处理

观察变量，对于人口和面积变量，我们可以用人口/面积将其变为一个人口密度变量，而经纬度只是表示各个街区的地理位置，不能表现街区的社会特征，因此，我们暂且将这两个变量剔除，并将数据进行标准化。

> data$density<-data$Population/data$Area
> data1<-data[,-c(12:15)]###去掉人口、面积、经纬度变量> head(data1)
                 LA.Nbhd Income Schools Diversity Age Homes Vets Asian Black1        Adams_Normandie  29606     691       0.6  26  0.26 0.05  0.05  0.252                 Arleta  65649     719       0.4  29  0.29 0.07  0.11  0.023      Arlington_Heights  31423     687       0.8  31  0.31 0.05  0.13  0.254        Atwater_Village  53872     762       0.9  34  0.34 0.06  0.20  0.015 Baldwin_Hills/Crenshaw  37948     656       0.4  36  0.36 0.10  0.05  0.716                Bel-Air 208861     924       0.2  46  0.46 0.13  0.08  0.01
  Latino White   density1   0.62  0.06 38835.0002   0.72  0.13 10021.9353   0.57  0.05 22106.0004   0.51  0.22  8271.1115   0.17  0.03 10041.0006   0.05  0.83  1201.212> scale_data<-scale(data1[,-1])###将数据标准化

三.主成分分析

从相关阵出发，计算相关阵的特征值、特征向量，根据累计贡献率选取主成分个数。

cor_matrix<-cor(scale_data)##相关矩阵eigen<-eigen(cor_matrix)###特征值、特征向量gongxian_rate<-eigen$values/sum(eigen$values)##方差贡献率cum_rate<-cumsum(eigen$values)/sum(eigen$values)###累计方差贡献率

结果如下：

> eigen
$values
 [1] 4.6930331513 1.8176421789 1.1687855390 1.0269743864 0.8062945744
 [6] 0.5577210494 0.4099288177 0.2437073056 0.1870581845 0.0882870577[11] 0.0005677551$vectors
             [,1]        [,2]         [,3]         [,4]         [,5]
 [1,] -0.40863471  0.04870608  0.204649345 -0.029705695 -0.026111021
 [2,] -0.03056206 -0.25663203  0.097169733  0.787452013  0.514858847
 [3,]  0.04582099 -0.50221591 -0.556053141 -0.005890927  0.005689376
 [4,] -0.42685697 -0.03623545 -0.119166660 -0.175554395  0.090259490
 [5,] -0.31637860 -0.04822810 -0.116671747  0.321273568 -0.319983287
 [6,] -0.36732010  0.20840539 -0.333535403  0.133826157 -0.047212706
 [7,] -0.02429262 -0.63666603 -0.159244974 -0.191916789 -0.160931505
 [8,]  0.14568670  0.46744106 -0.601051017  0.085945094  0.137411291
 [9,]  0.38404445 -0.03971545  0.229242391  0.223172109 -0.395285603[10,] -0.40399716 -0.06681935  0.247309411 -0.184073248  0.279824848[11,]  0.28326513 -0.05087439 -0.009350548 -0.314879884  0.584899674
             [,6]        [,7]        [,8]         [,9]       [,10]
 [1,]  0.11520700  0.02414449 -0.40720337  0.779078572 -0.02101538
 [2,] -0.05830369 -0.16370841 -0.06754548 -0.001444467 -0.06181177
 [3,] -0.23198685  0.56181935  0.01537205  0.223076194 -0.12960575
 [4,]  0.06102746 -0.16261473 -0.11938577 -0.281797151 -0.79960993
 [5,]  0.69514540  0.17211476  0.40286718 -0.027098252  0.04916595
 [6,] -0.03391422  0.14448585 -0.58384102 -0.408311715  0.40078456
 [7,]  0.14445665 -0.60618509 -0.13365058 -0.006753609  0.24867552
 [8,]  0.04955672 -0.30625699  0.11599645  0.280302229 -0.05727508
 [9,]  0.12290085  0.18728813 -0.38361188 -0.084606342 -0.26533256[10,] -0.15965016  0.20472602  0.28604268 -0.088487009  0.19728755[11,]  0.61440385  0.20442388 -0.22730006 -0.052971337  0.04554907
             [,11]
 [1,]  0.011160715
 [2,] -0.004331797
 [3,] -0.012713363
 [4,] -0.010332579
 [5,]  0.007493128
 [6,] -0.003923512
 [7,] -0.194680705
 [8,] -0.425663391> gongxian_rate
 [1] 0.4266393774 0.1652401981 0.1062532308 0.0933613079 0.0732995068
 [6] 0.0507019136 0.0372662562 0.0221552096 0.0170052895 0.0080260962[11] 0.0000516141> cum_rate
 [1] 0.4266394 0.5918796 0.6981328 0.7914941 0.8647936 0.9154955 0.9527618
 [8] 0.9749170 0.9919223 0.9999484 1.0000000

画碎石图：

plot(1:11,eigen$values,type="o",pch=17,col=4,main="Scree Plot",
xlab="Component Numbers",ylab="Eigen Values")

由累计方差贡献率和碎石图，前四个主成分的累计方差贡献率为79.15%，可提取原始变量中的大部分信息，因此只选取前4个主成分。主成分与变量的系数矩阵为前四个特征向量矩阵，可写出各主成分的表达式，结合系数的大小及实际意义，想找出各主成分代表的综合指标，但是并不是所有主成分的实际意义都比较好解释，这里的主成分就不太好解释。

下面计算出四个主成分与原始变量的相关系数（载荷）。

loading<-sweep(eigen$vectors,2,sqrt(eigen$values),"*")
> loading
             [,1]        [,2]        [,3]         [,4]         [,5]
 [1,] -0.88524211  0.06566552  0.22124720 -0.030103676 -0.023446105
 [2,] -0.06620785 -0.34599120  0.10505057  0.798001859  0.462311874
 [3,]  0.09926389 -0.67708730 -0.60115119 -0.005969851  0.005108713
 [4,] -0.92471774 -0.04885262 -0.12883153 -0.177906376  0.081047522
 [5,] -0.68538391 -0.06502110 -0.12613428  0.325577814 -0.287325496
 [6,] -0.79574058  0.28097207 -0.36058641  0.135619085 -0.042394133
 [7,] -0.05262609 -0.85835290 -0.17216035 -0.194487984 -0.144506686
 [8,]  0.31560707  0.63020386 -0.64979857  0.087096538  0.123386967
 [9,]  0.83197122 -0.05354435  0.24783483  0.226162046 -0.354942387[10,] -0.87519561 -0.09008582  0.26736716 -0.186539360  0.251265664[11,]  0.61364884 -0.06858884 -0.01010891 -0.319098470  0.525204269

绘制出载荷图：

plot(loading[,1:2],type="n",main="Loadings",xlab="Component1(42%)",
ylab="Component2(17%)",xlim=c(-1.2,1.2),ylim=c(-1.2,1.2))
name<-names(data1[,-1])
text(loading[,1],loading[,2],name)
abline(h=0)
abline(v=0)
plot(loading[,3:4],type="n",main="Loadings",xlab="Component3(11%)",
ylab="Component4(9%)",xlim=c(-1.2,1.2),ylim=c(-1.2,1.2))
text(loading[,3],loading[,4],name)
abline(h=0)
abline(v=0)

上面两幅图，变量在各主成分坐标轴上对应的数字表示该变量与该主成分的相关系数。可以看到，与第一主成分相关性大的变量有Income(负相关)、Age（负相关）、Homes（负相关）、Vets（负相关）、White（负相关）、density(正相关)、Lyatino（正相关），说明第一主成分小的对应高收入、高年龄、拥有住房的家庭多、复员军人多及欧裔多的社区，第一主成分大的对应拉美裔人口多、人口密度高的社区；同理，第二主成分大的对应非裔人口多的社区，第二主成分小的对应亚裔人较多、种族多样性较高的社区；第三主成分小的对应非裔较多、种族多样性较高的社区；第四主成分大的对应学校成绩较好的社区。

下面计算各个样本的主成分得分，即用标准化的原始数据乘特征向量矩阵。

> score<-as.matrix(scale_data)%*%eigen$vectors[,1:4]
> head(score)
            [,1]        [,2]       [,3]        [,4]
[1,]  3.14574473  0.36189223 -0.1621218 -0.74425898[2,]  0.92266517  0.01063399  0.9401550  0.35840978[3,]  1.98839943 -0.48717706 -0.8641400 -0.35372780[4,]  0.25666559 -1.70951387 -0.3518666  0.09088856[5,]  0.05765137  2.41953248 -2.4673846  0.24018071[6,] -5.79124555  0.77941752  1.2343587  0.65368578

根据主成分得分，将各个样本绘制出来：

plot(score[,1],score[,2],main="Sample Principal Components",xlab="Component1",
ylab="Component2",type="n",xlim=c(-7,6),ylim=c(-6,6))
head(data1)
text(score[,1],score[,2],data1[,1],cex=0.7)
abline(h=0,col=2)
abline(v=0,col=2)
plot(score[,3],score[,4],main="Sample Principal Components",xlab="Component3",
ylab="Component4",type="n",xlim=c(-7,6),ylim=c(-6,6))
head(data1)
text(score[,3],score[,4],data1[,1],cex=0.6)
abline(h=0,col=2)
abline(v=0,col=2)

将上图与主成分的载荷图一起来看，可以看出各个街区的分布特点，如Pacific_Palisades州在第一主成分上得分非常低，说明该州具有整体收入较高、年龄较大，拥有住房家庭多等特点；Gramercy_Park州在第二主成分上得分较高，说明该州的非裔较多。

四.因子分析

因子分析也是从原始数据的协方差阵或者相关阵出发，求特征根、特征向量，因子载荷矩阵与主成分的成分矩阵相差倍。估计载荷矩阵可以用主成分法、主因子解法、极大似然法，因子得分的估计可以用加权最小二乘估计、回归法。值得注意的是，因子载荷矩阵并不唯一。为了更方便地对公共因子进行解释，通常要对因子载荷矩阵进行旋转，常用最大方差旋转法，使得各因子载荷矩阵的每一列元素的平方按列向0或1两极转化。

根据主成分法提取公因子，选取4个公因子进行方差分析，使用回归法估计因子得分，并对因子载荷矩阵进行方差最大旋转。

> factor<-factanal(factors=4,scale_data,scores="regression",rotation="varimax")
> factor$loadings

Loadings:
          Factor1 Factor2 Factor3 Factor4
Income     0.643   0.493   0.274  -0.140 Schools                    0.157         Diversity                          0.545 Age        0.661   0.651   0.109   0.117 Homes      0.690   0.162                 Vets       0.790   0.379  -0.205  -0.126 Asian                      0.324   0.943 Black                     -0.986         Latino    -0.379  -0.913   0.116         White      0.368   0.806   0.431  -0.161 density   -0.587  -0.141  -0.120         

               Factor1 Factor2 Factor3 Factor4
SS loadings      2.587   2.349   1.454   1.286Proportion Var   0.235   0.214   0.132   0.117Cumulative Var   0.235   0.449   0.581   0.698

与主成分一样，画出载荷图：

可以看到，第一公因子大的对应高收入、拥有住房的家庭多、年龄较大、复员军人比例多的社区；

第二公因子对应欧裔比例高的社区；第三公因子主要描述了变量Black；第四公因子大对应亚裔人口比例大。

同样，由因子得分画出样本的分布图：

> summary(factor$scores)
    Factor1            Factor2           Factor3            Factor4       
 Min.   :-2.37461   Min.   :-1.8506   Min.   :-3.93076   Min.   :-1.2748  
 1st Qu.:-0.64708   1st Qu.:-0.6653   1st Qu.:-0.04804   1st Qu.:-0.6820  
 Median :-0.05652   Median :-0.2235   Median : 0.41276   Median :-0.2193  
 Mean   : 0.00000   Mean   : 0.0000   Mean   : 0.00000   Mean   : 0.0000  
 3rd Qu.: 0.53148   3rd Qu.: 0.8871   3rd Qu.: 0.52365   3rd Qu.: 0.4803  
 Max.   : 2.81319   Max.   : 2.1912   Max.   : 0.98750   Max.   : 3.7930 plot(factor$scores[,1:2],type="n",main="Loadings",xlab="Factor1",
ylab="Factor2",xlim=c(-2.5,3),ylim=c(-2,2.5))
text(factor$scores[,1],factor$scores[,2],data1[,1],cex=0.7)
abline(h=0)
abline(v=0)
plot(factor$scores[,3:4],type="n",main="Loadings",xlab="Factor3",
ylab="Factor4",xlim=c(-4.2,1.5),ylim=c(-1.5,4))
text(factor$scores[,3],factor$scores[,4],data1[,1],cex=0.7)
abline(h=0)
abline(v=0)