R语言多元分析系列

R语言多元分析系列之一：主成分分析

主成分分析（principal components analysis， PCA）是一种分析、简化数据集的技术。它把原始数据变换到一个新的坐标系统中，使得任何数据投影的第一大方差在第一个坐标（称为第一主成分）上，第二大方差在第二个坐标（第二主成分）上，依次类推。主成分分析经常用减少数据集的维数，同时保持数据集的对方差贡献最大的特征。这是通过保留低阶主成分，忽略高阶主成分做到的。这样低阶成分往往能够保留住数据的最重要方面。但是在处理观测数目小于变量数目时无法发挥作用，例如基因数据。

R语言中进行主成分分析可以采用基本的princomp函数，将结果输入到summary和plot函数中可分别得到分析结果和碎石图。但psych扩展包更具灵活性。

1 选择主成分个数

选择主成分个数通常有如下几种评判标准：

根据经验与理论进行选择

根据累积方差贡献率 ，例如选择使累积方差贡献率达到80%的主成分个数。

根据相关系数矩阵的特征值，选择特征值大于1的主成分。

另一种较为先进的方法是平行分析（parallel analysis）。该方法首先生成若干组与原始数据结构相同的随机矩阵，求出其特征值并进行平均，然后和真实数据的特征值进行比对，根据交叉点的位置来选择主成分个数。我们选择USJudgeRatings数据集举例，首先加载psych包，然后使用fa.parallel函数绘制下图，从图中可见第一主成分位于红线上方，第二主成分位于红线下方，因此主成分数目选择1。

fa.parallel(USJudgeRatings[,-1], fa="pc",n.iter=100, show.legend=FALSE) 

2 提取主成分

pc=principal(USJudgeRatings[,-1],nfactors=1) 
 
  
    PC1   h2     u2 
1  0.92 0.84 0.1565 
2  0.91 0.83 0.1663 
3  0.97 0.94 0.0613 
4  0.96 0.93 0.0720 
5  0.96 0.92 0.0763 
6  0.98 0.97 0.0299 
7  0.98 0.95 0.0469 
8  1.00 0.99 0.0091 
9  0.99 0.98 0.0196 
10 0.89 0.80 0.2013 
11 0.99 0.97 0.0275 
 
                 PC1 
SS loadings    10.13 
Proportion Var  0.92 

从上面的结果观察到，PC1即观测变量与主成分之间的相关系数，h2是变量能被主成分解释的比例，u2则是不能解释的比例。主成分解释了92%的总方差。注意此结果与princomp函数结果不同，princomp函数返回的是主成分的线性组合系数，而principal函数返回原始变量与主成分之间的相关系数，这样就和因子分析的结果意义相一致。

3 旋转主成分

旋转是在保持累积方差贡献率不变条件下，将主成分负荷进行变换，以方便解释。成分旋转这后各成分的方差贡献率将重新分配，此时就不可再称之为“主成分”而仅仅是“成分”。旋转又可分为正交旋转和斜交旋转。正交旋转的流行方法是方差最大化，需要在principal中增加rotate='varimax'参数加以实现。也有观点认为主成分分析一般不需要进行旋转。

4 计算主成分得分

主成分得分是各变量的线性组合，在计算出主成分得分之后，还可以将其进行回归等做进一步分析处理。但注意如果输入数据不是原始数据时，则无法计算主成分得分。我们需要在principal中增加score=T的参数设置，结果将存放在结果的score元素中。

R语言多元分析系列之二：探索性因子分析

探索性因子分析（Exploratory Factor Analysis，EFA）是一项用来找出多元观测变量的本质结构、并进行处理降维的技术。 因而EFA能够将具有错综复杂关系的变量综合为少数几个核心因子。EFA和PCA的区别在于：PCA中的主成分是原始变量的线性组合，而EFA中的原始变量是公共因子的线性组合，因子是影响变量的潜在变量，变量中不能被因子所解释的部分称为误差，因子和误差均不能直接观察到。进行EFA需要大量的样本，一般经验认为如何估计因子的数目为N，则需要有5N到10N的样本数目。

虽然EFA和PCA有本质上的区别，但在分析流程上有相似之处。下面我们用ability.cov这个心理测量数据举例，其变量是对人的六种能力，例如阅读和拼写能力进行了测验，其数据是一个协方差矩阵而非原始数据。R语言中stats包中的factanal函数可以完成这项工作，但这里我们使用更为灵活的psych包。

一、选择因子个数

一般选择因子个数可以根据相关系数矩阵的特征值，特征值大于0则可选择做为因子。我们仍使用平行分析法（parallel analysis）。该方法首先生成若干组与原始数据结构相同的随机矩阵，求出其特征值并进行平均，然后和真实数据的特征值进行比对，根据交叉点的位置来选择因子个数。根据下图我们可以观察到特征值与红线的关系，有两个因子都位于红线上方，显然应该选择两个因子。

library(psych) 
covariances = ability.cov$cov 
correlations = cov2cor(covariances) 
fa.parallel(correlations, n.obs=112, fa="fa", n.iter=100,show.legend=FALSE) 

二、提取因子

psych包中是使用fa函数来提取因子，将nfactors参数设定因子数为2，rotate参数设定了最大化方差的因子旋转方法，最后的fm表示分析方法，由于极大似然方法有时不能收敛，所以此处设为迭代主轴方法。从下面的结果中可以观察到两个因子解释了60%的总方差。Reading和vocabulary这两个变量于第一项因子有关，而picture、blocks和maze变量与第二项因子有关，general变量于两个因子都有关系。

fafa = fa(correlations,nfactors=2,rotate="varimax",fm="pa" ) 
         PA1  PA2   h2    u2 
general 0.49 0.57 0.57 0.432 
picture 0.16 0.59 0.38 0.623 
blocks  0.18 0.89 0.83 0.166 
maze    0.13 0.43 0.20 0.798 
reading 0.93 0.20 0.91 0.089 
vocab   0.80 0.23 0.69 0.313 
 
                PA1  PA2 
SS loadings    1.83 1.75 
Proportion Var 0.30 0.29 
Cumulative Var 0.30 0.60 

如果采用基本函数factanal进行因子分析，那么函数形式应该是factanal(covmat=correlations,factors=2,rottion='varimax')，这会得到相同的结果。此外，我们还可以用图形来表示因子和变量之间的关系

factor.plot(fa,labels=rownames(fa$loadings)) 

三、因子得分

得到公共因子后，我们可以象主成分分析那样反过来考察每个样本的因子得分。如果输入的是原始数据，则可以在fa函数中设置score=T参数来获得因子得分。如果象上面例子那样输入的是相关矩阵，则需要根据因子得分系数来回归估计。

fa$weights 
 
                 PA1         PA2 
general  0.017702900  0.21504415 
picture -0.007986044  0.09687725 
blocks  -0.198309764  0.79392660 
maze     0.019155930  0.03027495 
reading  0.841777373 -0.22404221 
vocab    0.190592536 -0.02040749 

参考资料：R in Action

R语言多元分析系列之三：多维标度分析

多维标度分析(MDS)是一种将多维空间的研究对象简化到低维空间进行定位、分析和归类，同时又保留对象间原始关系的数据分析方法。

设想一下如果我们在欧氏空间中已知一些点的座标，由此可以求出欧氏距离。那么反过来，已知距离应该也能得到这些点之间的关系。这种距离可以是古典的欧氏距离，也可以是广义上的“距离”。MDS就是在尽量保持这种高维度“距离”的同时，将数据在低维度上展现出来。从这种意义上来讲，主成分分析也是多维标度分析的一个特例。

一、距离的度量

多元分析中常用有以下几种距离，即绝对值距离、欧氏距离(euclidean)、马氏距离(manhattan)、 两项距离(binary)、明氏距离(minkowski)。在R中通常使用disk函数得到样本之间的距离。MDS就是对距离矩阵进行分析，以展现并解释数据的内在结构。

在经典MDS中，距离是数值数据表示，将其看作是欧氏距离。在R中stats包的cmdscale函数实现了经典MDS。它是根据各点的欧氏距离，在低维空间中寻找各点座标，而尽量保持距离不变。

非度量MDS方法中，“距离"不再看作数值数据，而只是顺序数据。例如在心理学实验中，受试者只能回答非常同意、同意、不同意、非常不同意这几种答案。在这种情况下，经典MDS不再有效。Kruskal在1964年提出了一种算法来解决这个问题。在R中MASS包的isoMDS函数可以实现这种算法，另一种流行的算法是由sammon函数实现的。

二、经典MDS

下面我们以HSAUR2包中的watervoles数据来举例。该数据是一个相似矩阵，表示了不同地区水田鼠的相似程度。首先加载数据然后用cmdscales进行分析。

library(ggplot2) 
data(watervoles, package = "HSAUR2") 
data(watervoles) 
voles.mds=cmdscale(watervoles,k=13,eig=T) 

下面计算前两个特征值在所有特征值中的比例，这是为了检测能否用两个维度的距离来表示高维空间中距离，如果达到了0.8左右则表示是合适的。

sum(abs(voles.mds$eig[1:2]))/sum(abs(voles.mds$eig)) 
sum((voles.mds$eig[1:2])^2)/sum((voles.mds$eig)^2) 

然后从结果中提取前两个维度的座标，用ggplot包进行绘图。 

x = voles.mds$points[,1] 
y = voles.mds$points[,2] 
p=ggplot(data.frame(x,y),aes(x,y,label = colnames(watervoles))) 
p+geom_point(shape=16,size=3,colour='red')+ 
  geom_text(hjust=-0.1,vjust=0.5,alpha=0.5) 

library("MASS") 
data(voting, package = "HSAUR2") 
voting_mds = isoMDS(voting) 
x = voting_mds$points[,1] 
y = voting_mds$points[,2] 
g=ggplot(data.frame(x,y),aes(x,y,label = colnames(voting))) 
g+geom_point(shape=16,size=3,colour='red')+ 
  geom_text(hjust=-0.1,vjust=0.5,alpha=0.5) 

library("MASS") 
data(voting, package = "HSAUR2") 
voting_mds = isoMDS(voting) 
x = voting_mds$points[,1] 
y = voting_mds$points[,2] 
g=ggplot(data.frame(x,y),aes(x,y,label = colnames(voting))) 
g+geom_point(shape=16,size=3,colour='red')+ 
  geom_text(hjust=-0.1,vjust=0.5,alpha=0.5) 

library(MASS) 
model1=lda(Species~.,data=iris,prior=c(1,1,1)/3) 
table(Species,predict(model1)$class) 
Species      setosa versicolor virginica 
setosa         50          0         0 
versicolor      0         48         2 
virginica       0          1        49 
从以上结果可观察到判断错误的样本只有三个。在判别函数建立后，还可以类似主成分分析那样对判别得分进行绘图 
ld=predict(model1)$x 
p=ggplot(cbind(iris,as.data.frame(ld)) 
,aes(x=LD1,y=LD2)) 
p+geom_point(aes(colour=Species),alpha=0.8,size=3) 

data=iris[,-5] 
distdist.e=dist(data,method='euclidean') 
heatmap(as.matrix(dist.e),labRow = F, labCol = F) 

model1=hclust(dist.e,method='ward') 
result=cutree(model1,k=3) 

mds=cmdscale(dist.e,k=2,eig=T) 
x = mds$points[,1] 
y = mds$points[,2] 
library(ggplot2) 
p=ggplot(data.frame(x,y),aes(x,y)) 
p+geom_point(size=3,alpha=0.8, 
             aes(colour=factor(result), 
               shape=iris$Species))