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變量于兩個因子都有關系。
fa = 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語言多元分析系列之三：多維標度分析

       多維標度分析(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)
R語言多元分析系列之四：判別分析
         判別分析（discriminant analysis）是一種分類技術。它通過一個已知類別的“訓練樣本”來建立判別準則，并通過預測變量來為未知類別的數據進行分類。
         判別分析的方法大體上有三類，即Fisher判別、Bayes判別和距離判別。Fisher判別思想是投影降維，使多維問題簡化為一維問題來處理。選擇一個適當的投影軸,使所有的樣品點都投影到這個軸上得到一個投影值。對這個投影軸的方向的要求是：使每一組內的投影值所形成的組內離差盡可能小，而不同組間的投影值所形成的類間離差盡可能大。Bayes判別思想是根據先驗概率求出后驗概率，并依據后驗概率分布作出統計推斷。距離判別思想是根據已知分類的數據計算各類別的重心，對未知分類的數據，計算它與各類重心的距離，與某個重心距離最近則歸于該類。
1.線性判別
          當不同類樣本的協方差矩陣相同時，我們可以在R中使用MASS包的lda函數實現線性判別。lda函數以Bayes判別思想為基礎。當分類只有兩種且總體服從多元正態分布條件下，Bayes判別與Fisher判別、距離判別是等價的。本例使用iris數據集來對花的品種進行分類。首先載入MASS包，建立判別模型，其中的prior參數表示先驗概率。然后利用table函數建立混淆矩陣，比對真實類別和預測類別。
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)
2.二次判別
           當不同類樣本的協方差矩陣不同時，則應該使用二次判別。
model2=qda(Species~.,data=iris,cv=T)
         這里將CV參數設置為T，是使用留一交叉檢驗（leave-one-out cross-validation），并自動生成預測值。這種條件下生成的混淆矩陣較為可靠。此外還可以使用predict(model)$posterior提取后驗概率。
          在使用lda和qda函數時注意：其假設是總體服從多元正態分布，若不滿足的話則謹慎使用。
R語言多元分析系列之五：聚類分析
            聚類分析（Cluster Analysis）是根據“物以類聚”的道理，對樣品或指標進行分類的一種多元統計分析方法，它是在沒有先驗知識的情況下，對樣本按各自的特性來進行合理的分類。
          聚類分析被應用于很多方面，在商業上，聚類分析被用來發現不同的客戶群，并且通過購買模式刻畫不同的客戶群的特征；在生物上，聚類分析被用來動植物分類和對基因進行分類，獲取對種群固有結構的認識；在因特網應用上，聚類分析被用來在網上進行文檔歸類來修復信息。
          聚類分析有兩種主要計算方法，分別是凝聚層次聚類（Agglomerative hierarchical method）和K均值聚類（K-Means）。
一、層次聚類
         層次聚類又稱為系統聚類，首先要定義樣本之間的距離關系，距離較近的歸為一類，較遠的則屬于不同的類?？捎糜诙x“距離”的統計量包括了歐氏距離(euclidean)、馬氏距離(manhattan)、 兩項距離(binary)、明氏距離(minkowski)。還包括相關系數和夾角余弦。
         層次聚類首先將每個樣本單獨作為一類，然后將不同類之間距離最近的進行合并，合并后重新計算類間距離。這個過程一直持續到將所有樣本歸為一類為止。在計算類間距離時則有六種不同的方法，分別是最短距離法、最長距離法、類平均法、重心法、中間距離法、離差平方和法。
        下面我們用iris數據集來進行聚類分析，在R語言中所用到的函數為hclust。首先提取iris數據中的4個數值變量，然后計算其歐氏距離矩陣。然后將矩陣繪制熱圖，從圖中可以看到顏色越深表示樣本間距離越近，大致上可以區分出三到四個區塊，其樣本之間比較接近。
data=iris[,-5]
dist.e=dist(data,method='euclidean')
heatmap(as.matrix(dist.e),labRow = F, labCol = F)
        然后使用hclust函數建立聚類模型，結果存在model1變量中，其中ward參數是將類間距離計算方法設置為離差平方和法。使用plot(model1)可以繪制出聚類樹圖。如果我們希望將類別設為3類，可以使用cutree函數提取每個樣本所屬的類別。 
model1=hclust(dist.e,method='ward')
result=cutree(model1,k=3)
        為了顯示聚類的效果，我們可以結合多維標度和聚類的結果。先將數據用MDS進行降維，然后以不同的的形狀表示原本的分類，用不同的顏色來表示聚類的結果?？梢钥吹絪etose品種聚類很成功，但有一些virginica品種的花被錯誤和virginica品種聚類到一起。 
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))
二、K均值聚類
            K均值聚類又稱為動態聚類，它的計算方法較為簡單，也不需要輸入距離矩陣。首先要指定聚類的分類個數N，隨機取N個樣本作為初始類的中心，計算各樣本與類中心的距離并進行歸類，所有樣本劃分完成后重新計算類中心，重復這個過程直到類中心不再變化。
            在R中使用kmeans函數進行K均值聚類，centers參數用來設置分類個數，nstart參數用來設置取隨機初始中心的次數，其默認值為1，但取較多的次數可以改善聚類效果。model2$cluster可以用來提取每個樣本所屬的類別。 
model2=kmeans(data,centers=3,nstart=10)
         使用K均值聚類時需要注意，只有在類的平均值被定義的情況下才能使用，還要求事先給出分類個數。一種方法是先用層次聚類以決定個數，再用K均值聚類加以改進?；蛘咭暂喞禂祦砼袛喾诸悅€數。改善聚類的方法還包括對原始數據進行變換，如對數據進行降維后再實施聚類。

cluster擴展包中也有許多函數可用于聚類分析，如agnes函數可用于凝聚層次聚類，diana可用于劃分層次聚類，pam可用于K均值聚類，fanny用于模糊聚類。