R語言多元分析系列

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ù)交叉點(diǎn)的位置來選擇主成分個數(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)。也有觀點(diǎ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ù)交叉點(diǎn)的位置來選擇因子個數(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)系。
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
         如果采用基本函數(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語言多元分析系列之三：多維標(biāo)度分析

       多維標(biāo)度分析(MDS)是一種將多維空間的研究對象簡化到低維空間進(jìn)行定位、分析和歸類，同時又保留對象間原始關(guān)系的數(shù)據(jù)分析方法。
       設(shè)想一下如果我們在歐氏空間中已知一些點(diǎn)的座標(biāo)，由此可以求出歐氏距離。那么反過來，已知距離應(yīng)該也能得到這些點(diǎn)之間的關(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ù)各點(diǎn)的歐氏距離，在低維空間中尋找各點(diǎn)座標(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)

三、非度量MDS

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

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