R语言多元分析系列

R语言多元分析系列R语言多元分析系列之一：主成分分析主成分分析（principalcomponentsanalysis，PCA）是一种分析、简化数据集的技术。它把原始数据变换到一个新的坐标系统中，使得任何数据投影的第一大方差在第一个坐标（称为第一主成分）上，第二大方差在第二个坐标（第二主成分）上，依次类推。主成分分析经常用减少数据集的维数，同时保持数据集的对方差贡献最大的特征。这是通过保留低阶主成分，忽略高阶主成分做到的。这样低阶成分往往能够保留住数据的最重要方面。但是在处理观测数目小于变量数目时无法发挥作用，例如基因数据。R语言中进行主成分分析可以采用基本的princomp函数，将结果输入到summary和plot函数中可分别得到分析结果和碎石图。但psych扩展包更具灵活性。1选择主成分个数选择主成分个数通常有如下几种评判标准：根据经验与理论进行选择根据累积方差贡献率，例如选择使累积方差贡献率达到80%的主成分个数。根据相关系数矩阵的特征值，选择特征值大于1的主成分。另一种较为先进的方法是平行分析（parallelanalysis）。该方法首先生成若干组与原始数据结构相同的随机矩阵，求出其特征值并进行平均，然后和真实数据的特征值进行比对，根据交叉点的位置来选择主成分个数。我们选择USJudgeRatings数据集举例，首先加载psych包，然后使用fa.parallel函数绘制下图，从图中可见第一主成分位于红线上方，第二主成分位于红线下方，因此主成分数目选择1。fa.parallel(USJudgeRatings[,-1],fa=pc,n.iter=100,show.legend=FALSE)2提取主成分pc=principal(USJudgeRatings[,-1],nfactors=1)PC1h2u210.920.840.156520.910.830.166330.970.940.061340.960.930.072050.960.920.076360.980.970.029970.980.950.046981.000.990.009190.990.980.0196100.890.800.2013110.990.970.0275PC1SSloadings10.13ProportionVar0.92从上面的结果观察到，PC1即观测变量与主成分之间的相关系数，h2是变量能被主成分解释的比例，u2则是不能解释的比例。主成分解释了92%的总方差。注意此结果与princomp函数结果不同，princomp函数返回的是主成分的线性组合系数，而principal函数返回原始变量与主成分之间的相关系数，这样就和因子分析的结果意义相一致。3旋转主成分旋转是在保持累积方差贡献率不变条件下，将主成分负荷进行变换，以方便解释。成分旋转这后各成分的方差贡献率将重新分配，此时就不可再称之为“主成分”而仅仅是“成分”。旋转又可分为正交旋转和斜交旋转。正交旋转的流行方法是方差最大化，需要在principal中增加rotate='varimax'参数加以实现。也有观点认为主成分分析一般不需要进行旋转。4计算主成分得分主成分得分是各变量的线性组合，在计算出主成分得分之后，还可以将其进行回归等做进一步分析处理。但注意如果输入数据不是原始数据时，则无法计算主成分得分。我们需要在principal中增加score=T的参数设置，结果将存放在结果的score元素中。R语言多元分析系列之二：探索性因子分析探索性因子分析（ExploratoryFactorAnalysis，EFA）是一项用来找出多元观测变量的本质结构、并进行处理降维的技术。因而EFA能够将具有错综复杂关系的变量综合为少数几个核心因子。EFA和PCA的区别在于：PCA中的主成分是原始变量的线性组合，而EFA中的原始变量是公共因子的线性组合，因子是影响变量的潜在变量，变量中不能被因子所解释的部分称为误差，因子和误差均不能直接观察到。进行EFA需要大量的样本，一般经验认为如何估计因子的数目为N，则需要有5N到10N的样本数目。虽然EFA和PCA有本质上的区别，但在分析流程上有相似之处。下面我们用ability.cov这个心理测量数据举例，其变量是对人的六种能力，例如阅读和拼写能力进行了测验，其数据是一个协方差矩阵而非原始数据。R语言中stats包中的factanal函数可以完成这项工作，但这里我们使用更为灵活的psych包。一、选择因子个数一般选择因子个数可以根据相关系数矩阵的特征值，特征值大于0则可选择做为因子。我们仍使用平行分析法（parallelanalysis）。该方法首先生成若干组与原始数据结构相同的随机矩阵，求出其特征值并进行平均，然后和真实数据的特征值进行比对，根据交叉点的位置来选择因子个数。根据下图我们可以观察到特征值与红线的关系，有两个因子都位于红线上方，显然应该选择两个因子。library(psych)covariances=ability.cov$covcorrelations=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变量于两个因子都有关系。fa=fa(correlations,nfactors=2,rotate=varimax,fm=pa)PA1PA2h2u2general0.490.570.570.432picture0.160.590.380.623blocks0.180.890.830.166maze0.130.430.200.798reading0.930.200.910.089vocab0.800.230.690.313PA1PA2SSloadings1.831.75ProportionVar0.300.29CumulativeVar0.300.60如果采用基本函数factanal进行因子分析，那么函数形式应该是factanal(covmat=correlations,factors=2,rottion='varimax')，这会得到相同的结果。此外，我们还可以用图形来表示因子和变量之间的关系factor.plot(fa,labels=rownames(fa$loadings))三、因子得分得到公共因子后，我们可以象主成分分析那样反过来考察每个样本的因子得分。如果输入的是原始数据，则可以在fa函数中设置score=T参数来获得因子得分。如果象上面例子那样输入的是相关矩阵，则需要根据因子得分系数来回归估计。fa$weightsPA1PA2general0.0177029000.21504415picture-0.0079860440.09687725blocks-0.1983097640.79392660maze0.0191559300.03027495reading0.841777373-0.22404221vocab0.190592536-0.02040749参考资料：RinActionR语言多元分析系列之三：多维标度分析多维标度分析(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)三、非度量MDS第二例子中的数据是关于新泽西州议员投票行为的相似矩阵，这里我们用MASS包中的isoMDS函数进行分析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)参考资料：AHandbookofStatisticalAnalysesUsingR多元统计分析及R语言建模R语言多元分析系列之四：判别分析判别分析（discriminantanalysis）是一种分类技术。它通过一个已知类别的“训练样本”来建立判别准则，并通过预测变量来为未知类别的数据进行分类。判别分析的方法大体上有三类，即Fisher判别、Bayes判别和距离判别。Fisher判别思想是投影降维，使多维问题简化为一维问题来处理。选择一个适当的投影轴,使所有的样品点都投影到这个轴上得到一个投影值。对这个投影轴的方向的要求是：使每一组内的投影值所形成的组内离差尽可能小，而不同组间的投影值所形成的类间离差尽可能大。Bayes判别思想是根据先验概率求出后