R語言應用統計1 主成分分析

這個系列就討論應用基礎，爭取一條公式都不用寫。當原始數據集比較龐大，并且不同變量之間存在一些相關性時，我們希望可以用更少的變量來表示原始數據集，用到的變量越少的同時，能夠表示的原始數據集中的信息越多自然就更好。主成分分析就可以實現這樣的目標，在主成分分析中用來表示原始數據集中的信息的變量被稱為主成分。下面我們用一個例子說明R語言中進行簡單的主成分分析的方法。

---

數據
使用HSAUR2包中的美國城市污染數據，代碼如下

install.packages("HSAUR2") data("USairpollution", package = "HSAUR2")

這個數據集有七個變量，分別是代表空氣污染的SO2（空氣中的二氧化硫濃度）、代表人類活動的popul、manu以及代表天氣情況的wind、temp、precip、 predays，這個數據集可以用多元線性回歸來分析代表人類活動與天氣情況的六個變量對二氧化硫濃度的影響，但現在我們嘗試用PCA對六個解釋變量進行分析。

數據的定性分析
首先分析一下這六個解釋變量之間的相關性，使用cor函數，代碼為

cor(USairpollution[,-1])

代碼輸出如下

temp manu popul wind precip predays temp 1.00000000 -0.19004216 -0.06267813 -0.34973963 0.38625342 -0.43024212 manu -0.19004216 1.00000000 0.95526935 0.23794683 -0.03241688 0.13182930 popul -0.06267813 0.95526935 1.00000000 0.21264375 -0.02611873 0.04208319 wind -0.34973963 0.23794683 0.21264375 1.00000000 -0.01299438 0.16410559 precip 0.38625342 -0.03241688 -0.02611873 -0.01299438 1.00000000 0.49609671 predays -0.43024212 0.13182930 0.04208319 0.16410559 0.49609671 1.00000000

可以發現代表人類活動的popul與manu的相關性系數高達0.9553，代表天氣情況的四個變量之間也存在一定的相關性，這說明如果直接用多元線性回歸分析這個數據集可能存在多重共線性的問題。

接下來我們畫出這六個變量兩兩之間的散點圖，代碼如下

panel.hist <- function(x, ...) {usr <- par("usr"); on.exit(par(usr))par(usr = c(usr[1:2], 0, 1.5) )h <- hist(x, plot = FALSE)breaks <- h$breaks; nB <- length(breaks)y <- h$counts; y <- y/max(y)rect(breaks[-nB], 0, breaks[-1], y, col="grey", ...) } USairpollution$negtemp <- USairpollution$temp * (-1) USairpollution$temp <- NULL pdf("Lot0.pdf") #將圖片以pdf形式導出到當前工作路徑 pairs(USairpollution[,-1], diag.panel = panel.hist,pch = ".", cex = 1.5) dev.off() #結束導出

從這些散點圖中可以看出，至少popul與manu之間，precip與negtemp之間存在比較明顯的線性相關性。negtemp就是取temp的相反數，可以這樣做但沒必要，只是取相反數畫出來圖看上去線性性要明顯一點而已。

主成分分析
使用princomp函數做主成分分析并用summary查看結果，代碼如下，

usair_pca <- princomp(USairpollution[,-1], cor = TRUE) summary(usair_pca, loadings = TRUE)

cor=TRUE代表在計算中需要用到相關性矩陣，因為在定性分析的時候我們發現六個變量之間確實存在一定的相關性，所以需要用到，如果大家在嘗試其他數據集，發現相關性系數矩陣非常接近單位矩陣，就可以用cor=FALSE，如果變量中存在常值變量，就只能用cor=FALSE；loadings=TRUE代表在展示結果時同時展示各個主成分的載荷。代碼輸出如下

Importance of components:Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Standard deviation 1.4819456 1.2247218 1.1809526 0.8719099 0.33848287 0.185599752 Proportion of Variance 0.3660271 0.2499906 0.2324415 0.1267045 0.01909511 0.005741211 Cumulative Proportion 0.3660271 0.6160177 0.8484592 0.9751637 0.99425879 1.000000000Loadings:Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 temp 0.330 0.128 0.672 0.306 0.558 0.136 manu -0.612 0.168 0.273 -0.137 -0.102 0.703 popul -0.578 0.222 0.350 -0.695 wind -0.354 -0.131 -0.297 0.869 0.113 precip -0.623 0.505 0.171 -0.568 predays -0.238 -0.708 -0.311 0.580

第二行Proportion of Variance代表了每一個主成分包含的原始數據集的信息，第三行Cumulative Proportion代表前幾個主成分累積包含的原始數據集的信息，比如假設選擇主成分的標準是被選擇的主成分包含的信息超過原始數據集的95%，則我們應該選前四個主成分，這樣就實現了用四個變量代表原始的六個變量的簡化目標，并且根據Loadings，我們可以用原始的六個變量的線性組合來表示這四個主成分，以第一主成分為例：
$$\begin{gathered} Comp1=0.330temp?0.612manu?0.578popul \\ ?0.354wind?0.238predays \end{gathered}$$

接下來我們分析這個主成分分析的score，score是每一個主成分在單個數據點上的得分，得分絕對值越高說明這個主成分在這個數據點上的表現越差，我們用MVA包中的bvbox函數畫出score的散點圖，

install.packages('MVA') library(MVA)

代碼如下，

pdf("Lot1.pdf") pairs(usair_pca$scores[,1:4], ylim = c(-6, 4), xlim = c(-6, 4), panel = function(x,y, ...) {text(x, y, abbreviate(row.names(USairpollution)),cex = 0.6)bvbox(cbind(x,y), add = TRUE)}) #[,1:4]代表分析前四個主成分 dev.off()

bvbox的作用是畫二維的箱線圖，在外圈虛線以外的點被認為是異常值，可以發現大部分數據點的信息都是可以被這四個主成分所代表的。

分析主成分與二氧化硫濃度的關系
現在我們用這些結果討論主成分與二氧化硫濃度之間的關系，此時被解釋變量是二氧化硫濃度，解釋變量是主成分的score，代碼如下

usair_reg <- lm(SO2 ~ usair_pca$scores, data = USairpollution) summary(usair_reg)

輸出如下

Call: lm(formula = SO2 ~ usair_pca$scores, data = USairpollution)Residuals:Min 1Q Median 3Q Max -23.004 -8.542 -0.991 5.758 48.758 Coefficients:Estimate Std. Error t value Pr(>|t|) (Intercept) 30.049 2.286 13.146 6.91e-15 *** usair_pca$scoresComp.1 -9.942 1.542 -6.446 2.28e-07 *** usair_pca$scoresComp.2 -2.240 1.866 -1.200 0.23845 usair_pca$scoresComp.3 0.375 1.935 0.194 0.84752 usair_pca$scoresComp.4 -8.549 2.622 -3.261 0.00253 ** usair_pca$scoresComp.5 -15.176 6.753 -2.247 0.03122 * usair_pca$scoresComp.6 39.271 12.316 3.189 0.00306 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 14.64 on 34 degrees of freedom Multiple R-squared: 0.6695, Adjusted R-squared: 0.6112 F-statistic: 11.48 on 6 and 34 DF, p-value: 5.419e-07

可以發現這個多元線性回歸模型是顯著的，第1個主成分是顯著的，在更弱一點的顯著性要求下，第4、5、6個主成分也是顯著的。R方說明用這六個主成分可以解釋二氧化硫濃度變化信息的61.12%。

總結

以上是生活随笔為你收集整理的R语言应用统计1 主成分分析的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。