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

主成分分析（principal components analysis， PCA）是一種分析、簡化數(shù)據(jù)集的技術(shù)。它把原始數(shù)據(jù)變換到一個新的坐標(biāo)系統(tǒng)中，使得任何數(shù)據(jù)投影的第一大方差在第一個坐標(biāo)（稱為第一主成分）上，第二大方差在第二個坐標(biāo)（第二主成分）上，依次類推。主成分分析經(jīng)常用減少數(shù)據(jù)集的維數(shù)，同時保持?jǐn)?shù)據(jù)集的對方差貢獻(xiàn)最大的特征。這是通過保留低階主成分，忽略高階主成分做到的。這樣低階成分往往能夠保留住數(shù)據(jù)的最重要方面。但是在處理觀測數(shù)目小于變量數(shù)目時無法發(fā)揮作用，例如基因數(shù)據(jù)。

R語言中進(jìn)行主成分分析可以采用基本的princomp函數(shù)，將結(jié)果輸入到summary和plot函數(shù)中可分別得到分析結(jié)果和碎石圖。但psych擴(kuò)展包更具靈活性。

1 選擇主成分個數(shù)

選擇主成分個數(shù)通常有如下幾種評判標(biāo)準(zhǔn)：

根據(jù)經(jīng)驗與理論進(jìn)行選擇

根據(jù)累積方差貢獻(xiàn)率 ，例如選擇使累積方差貢獻(xiàn)率達(dá)到80%的主成分個數(shù)。

根據(jù)相關(guān)系數(shù)矩陣的特征值，選擇特征值大于1的主成分。

另一種較為先進(jìn)的方法是平行分析（parallel analysis）。該方法首先生成若干組與原始數(shù)據(jù)結(jié)構(gòu)相同的隨機(jī)矩陣，求出其特征值并進(jìn)行平均，然后和真實數(shù)據(jù)的特征值進(jìn)行比對，根據(jù)交叉點的位置來選擇主成分個數(shù)。我們選擇USJudgeRatings數(shù)據(jù)集舉例，首先加載psych包，然后使用fa.parallel函數(shù)繪制下圖，從圖中可見第一主成分位于紅線上方，第二主成分位于紅線下方，因此主成分?jǐn)?shù)目選擇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 

從上面的結(jié)果觀察到，PC1即觀測變量與主成分之間的相關(guān)系數(shù)，h2是變量能被主成分解釋的比例，u2則是不能解釋的比例。主成分解釋了92%的總方差。注意此結(jié)果與princomp函數(shù)結(jié)果不同，princomp函數(shù)返回的是主成分的線性組合系數(shù)，而principal函數(shù)返回原始變量與主成分之間的相關(guān)系數(shù)，這樣就和因子分析的結(jié)果意義相一致。

3 旋轉(zhuǎn)主成分

旋轉(zhuǎn)是在保持累積方差貢獻(xiàn)率不變條件下，將主成分負(fù)荷進(jìn)行變換，以方便解釋。成分旋轉(zhuǎn)這后各成分的方差貢獻(xiàn)率將重新分配，此時就不可再稱之為“主成分”而僅僅是“成分”。旋轉(zhuǎn)又可分為正交旋轉(zhuǎn)和斜交旋轉(zhuǎn)。正交旋轉(zhuǎn)的流行方法是方差最大化，需要在principal中增加rotate='varimax'參數(shù)加以實現(xiàn)。也有觀點認(rèn)為主成分分析一般不需要進(jìn)行旋轉(zhuǎn)。

4 計算主成分得分

主成分得分是各變量的線性組合，在計算出主成分得分之后，還可以將其進(jìn)行回歸等做進(jìn)一步分析處理。但注意如果輸入數(shù)據(jù)不是原始數(shù)據(jù)時，則無法計算主成分得分。我們需要在principal中增加score=T的參數(shù)設(shè)置，結(jié)果將存放在結(jié)果的score元素中。

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

探索性因子分析（Exploratory Factor Analysis，EFA）是一項用來找出多元觀測變量的本質(zhì)結(jié)構(gòu)、并進(jìn)行處理降維的技術(shù)。 因而EFA能夠?qū)⒕哂绣e綜復(fù)雜關(guān)系的變量綜合為少數(shù)幾個核心因子。EFA和PCA的區(qū)別在于：PCA中的主成分是原始變量的線性組合，而EFA中的原始變量是公共因子的線性組合，因子是影響變量的潛在變量，變量中不能被因子所解釋的部分稱為誤差，因子和誤差均不能直接觀察到。進(jìn)行EFA需要大量的樣本，一般經(jīng)驗認(rèn)為如何估計因子的數(shù)目為N，則需要有5N到10N的樣本數(shù)目。

雖然EFA和PCA有本質(zhì)上的區(qū)別，但在分析流程上有相似之處。下面我們用ability.cov這個心理測量數(shù)據(jù)舉例，其變量是對人的六種能力，例如閱讀和拼寫能力進(jìn)行了測驗，其數(shù)據(jù)是一個協(xié)方差矩陣而非原始數(shù)據(jù)。R語言中stats包中的factanal函數(shù)可以完成這項工作，但這里我們使用更為靈活的psych包。

一、選擇因子個數(shù)

一般選擇因子個數(shù)可以根據(jù)相關(guān)系數(shù)矩陣的特征值，特征值大于0則可選擇做為因子。我們?nèi)允褂闷叫蟹治龇ǎ╬arallel analysis）。該方法首先生成若干組與原始數(shù)據(jù)結(jié)構(gòu)相同的隨機(jī)矩陣，求出其特征值并進(jìn)行平均，然后和真實數(shù)據(jù)的特征值進(jìn)行比對，根據(jù)交叉點的位置來選擇因子個數(shù)。根據(jù)下圖我們可以觀察到特征值與紅線的關(guān)系，有兩個因子都位于紅線上方，顯然應(yīng)該選擇兩個因子。

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函數(shù)來提取因子，將nfactors參數(shù)設(shè)定因子數(shù)為2，rotate參數(shù)設(shè)定了最大化方差的因子旋轉(zhuǎn)方法，最后的fm表示分析方法，由于極大似然方法有時不能收斂，所以此處設(shè)為迭代主軸方法。從下面的結(jié)果中可以觀察到兩個因子解釋了60%的總方差。Reading和vocabulary這兩個變量于第一項因子有關(guān)，而picture、blocks和maze變量與第二項因子有關(guān)，general變量于兩個因子都有關(guān)系。

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 

如果采用基本函數(shù)factanal進(jìn)行因子分析，那么函數(shù)形式應(yīng)該是factanal(covmat=correlations,factors=2,rottion='varimax')，這會得到相同的結(jié)果。此外，我們還可以用圖形來表示因子和變量之間的關(guān)系

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

三、因子得分

得到公共因子后，我們可以象主成分分析那樣反過來考察每個樣本的因子得分。如果輸入的是原始數(shù)據(jù)，則可以在fa函數(shù)中設(shè)置score=T參數(shù)來獲得因子得分。如果象上面例子那樣輸入的是相關(guān)矩陣，則需要根據(jù)因子得分系數(shù)來回歸估計。

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語言多元分析系列之三：多維標(biāo)度分析

多維標(biāo)度分析(MDS)是一種將多維空間的研究對象簡化到低維空間進(jìn)行定位、分析和歸類，同時又保留對象間原始關(guān)系的數(shù)據(jù)分析方法。

設(shè)想一下如果我們在歐氏空間中已知一些點的座標(biāo)，由此可以求出歐氏距離。那么反過來，已知距離應(yīng)該也能得到這些點之間的關(guān)系。這種距離可以是古典的歐氏距離，也可以是廣義上的“距離”。MDS就是在盡量保持這種高維度“距離”的同時，將數(shù)據(jù)在低維度上展現(xiàn)出來。從這種意義上來講，主成分分析也是多維標(biāo)度分析的一個特例。

一、距離的度量

多元分析中常用有以下幾種距離，即絕對值距離、歐氏距離(euclidean)、馬氏距離(manhattan)、 兩項距離(binary)、明氏距離(minkowski)。在R中通常使用disk函數(shù)得到樣本之間的距離。MDS就是對距離矩陣進(jìn)行分析，以展現(xiàn)并解釋數(shù)據(jù)的內(nèi)在結(jié)構(gòu)。

在經(jīng)典MDS中，距離是數(shù)值數(shù)據(jù)表示，將其看作是歐氏距離。在R中stats包的cmdscale函數(shù)實現(xiàn)了經(jīng)典MDS。它是根據(jù)各點的歐氏距離，在低維空間中尋找各點座標(biāo)，而盡量保持距離不變。

非度量MDS方法中，“距離"不再看作數(shù)值數(shù)據(jù)，而只是順序數(shù)據(jù)。例如在心理學(xué)實驗中，受試者只能回答非常同意、同意、不同意、非常不同意這幾種答案。在這種情況下，經(jīng)典MDS不再有效。Kruskal在1964年提出了一種算法來解決這個問題。在R中MASS包的isoMDS函數(shù)可以實現(xiàn)這種算法，另一種流行的算法是由sammon函數(shù)實現(xiàn)的。

二、經(jīng)典MDS

下面我們以HSAUR2包中的watervoles數(shù)據(jù)來舉例。該數(shù)據(jù)是一個相似矩陣，表示了不同地區(qū)水田鼠的相似程度。首先加載數(shù)據(jù)然后用cmdscales進(jìn)行分析。

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

下面計算前兩個特征值在所有特征值中的比例，這是為了檢測能否用兩個維度的距離來表示高維空間中距離，如果達(dá)到了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) 

然后從結(jié)果中提取前兩個維度的座標(biāo)，用ggplot包進(jìn)行繪圖。 

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 
從以上結(jié)果可觀察到判斷錯誤的樣本只有三個。在判別函數(shù)建立后，還可以類似主成分分析那樣對判別得分進(jìn)行繪圖 
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))