R語言應(yīng)用統(tǒng)計(jì)1 主成分分析

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

---

數(shù)據(jù)
使用HSAUR2包中的美國城市污染數(shù)據(jù)，代碼如下

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

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

數(shù)據(jù)的定性分析
首先分析一下這六個(gè)解釋變量之間的相關(guān)性，使用cor函數(shù)，代碼為

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

可以發(fā)現(xiàn)代表人類活動(dòng)的popul與manu的相關(guān)性系數(shù)高達(dá)0.9553，代表天氣情況的四個(gè)變量之間也存在一定的相關(guān)性，這說明如果直接用多元線性回歸分析這個(gè)數(shù)據(jù)集可能存在多重共線性的問題。

接下來我們畫出這六個(gè)變量兩兩之間的散點(diǎn)圖，代碼如下

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形式導(dǎo)出到當(dāng)前工作路徑 pairs(USairpollution[,-1], diag.panel = panel.hist,pch = ".", cex = 1.5) dev.off() #結(jié)束導(dǎo)出

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

主成分分析
使用princomp函數(shù)做主成分分析并用summary查看結(jié)果，代碼如下，

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

cor=TRUE代表在計(jì)算中需要用到相關(guān)性矩陣，因?yàn)樵诙ㄐ苑治龅臅r(shí)候我們發(fā)現(xiàn)六個(gè)變量之間確實(shí)存在一定的相關(guān)性，所以需要用到，如果大家在嘗試其他數(shù)據(jù)集，發(fā)現(xiàn)相關(guān)性系數(shù)矩陣非常接近單位矩陣，就可以用cor=FALSE，如果變量中存在常值變量，就只能用cor=FALSE；loadings=TRUE代表在展示結(jié)果時(shí)同時(shí)展示各個(gè)主成分的載荷。代碼輸出如下

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代表了每一個(gè)主成分包含的原始數(shù)據(jù)集的信息，第三行Cumulative Proportion代表前幾個(gè)主成分累積包含的原始數(shù)據(jù)集的信息，比如假設(shè)選擇主成分的標(biāo)準(zhǔn)是被選擇的主成分包含的信息超過原始數(shù)據(jù)集的95%，則我們應(yīng)該選前四個(gè)主成分，這樣就實(shí)現(xiàn)了用四個(gè)變量代表原始的六個(gè)變量的簡化目標(biāo)，并且根據(jù)Loadings，我們可以用原始的六個(gè)變量的線性組合來表示這四個(gè)主成分，以第一主成分為例：
$$\begin{gathered} Comp1=0.330temp?0.612manu?0.578popul \\ ?0.354wind?0.238predays \end{gathered}$$

接下來我們分析這個(gè)主成分分析的score，score是每一個(gè)主成分在單個(gè)數(shù)據(jù)點(diǎn)上的得分，得分絕對(duì)值越高說明這個(gè)主成分在這個(gè)數(shù)據(jù)點(diǎn)上的表現(xiàn)越差，我們用MVA包中的bvbox函數(shù)畫出score的散點(diǎn)圖，

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]代表分析前四個(gè)主成分 dev.off()

bvbox的作用是畫二維的箱線圖，在外圈虛線以外的點(diǎn)被認(rèn)為是異常值，可以發(fā)現(xiàn)大部分?jǐn)?shù)據(jù)點(diǎn)的信息都是可以被這四個(gè)主成分所代表的。

分析主成分與二氧化硫濃度的關(guān)系
現(xiàn)在我們用這些結(jié)果討論主成分與二氧化硫濃度之間的關(guān)系，此時(shí)被解釋變量是二氧化硫濃度，解釋變量是主成分的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

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

總結(jié)

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

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。