主成分分析是一門降維的技術，即將多個指標用少數幾個綜合指標表示出來。主成分分析可以用於變數的降維和數據的解釋。前言主成分分析（PCA)是一種經典的線性維歸約方法，基於高維空間的數據尋找主成分將其轉換至低維空間，並保證低維空間下每個維度具有最大方差且不相關性。目錄1. 什麼是主成分分析2. 為什麼降維3. 主成分滿足條件4. 分析過程及數據模型5. 案例一、什麼是主成分分析主成分分析，是一種降維的分析方法，其考察多個變數間相關性的一種多元統計方法，研究如何通過少數幾個主成分來揭示多個變數間的內部結構，即從原始變數中導出少數幾個主成分，使它們儘可能多地保留原始變數的信息，且彼此間互不相關。舉個栗子好理解若在數據收集過程中有許多變數,這些變數可能有幾十上百個,那麼怎麼去理解這些變數間的關係了?如果兩兩去看,那得有幾百個相關關係了,另外我們還會遇到這樣的問題:1、 比如拿到一個汽車的樣本，裡面既有以「千米/每小時」度量的最大速度特徵，也有「英里/小時」的最大速度特徵，顯然這兩個特徵有一個多餘。2、 拿到一個數學系的本科生期末考試成績單，裡面有三列，一列是對數學的興趣程度，一列是複習時間，還有一列是考試成績。我們知道要學好數學，需要有濃厚的興趣，所以第二項與第一項強相關，第三項和第二項也是強相關。那是不是可以合併第一項和第二項呢？3、 拿到一個樣本，特徵非常多，而樣例特別少，這樣用回歸去直接擬合非常困難，容易過度擬合。比如北京的房價：假設房子的特徵是（大小、位置、朝向、是否學區房、建造年代、是否二手、層數、所在層數），搞了這麼多特徵，結果只有不到十個房子的樣例。要擬合房子特徵->房價的這麼多特徵，就會造成過度擬合。4、 這個與第二個有點類似，假設在IR中我們建立的文檔-詞項矩陣中，有兩個詞項為「learn」和「study」，在傳統的向量空間模型中，認為兩者獨立。然而從語義的角度來講，兩者是相似的，而且兩者出現頻率也類似，是不是可以合成為一個特徵呢？主成分分析便是一種降維的技巧,就是將大量相關的變數變成一組很少的不相關的變數,這些無關變數稱之為主成分.二、為什麼要降維脫水版：1）多重共線性--預測變數之間存在一定程度的相關性。多重共線性會導致解空間的不穩定，從而可能導致結果的不連貫。2）高維空間本身具有稀疏性。3）過多的變數會妨礙查找規律的建立。4）僅在變數層面上分析可能會忽略變數之間的潛在聯繫。例如幾個預測變數的綁定才可以反映數據某一方面特徵。拖沓版：上面例子中有幾個關鍵詞,大量相關的變數,很少不相關的變數.學過線性代數的應該了解這叫求最大線性無關組.其實把每個變數當做一個人,相關就是指兩個人認識比較熟,不相關就是比較陌生.我們認為熟悉的人之間可以互相代表,所以若一組人之間都認識那麼只需要一個人就可以代表這個組,那麼最大線性無關組就是變成組裡面只剩下相互陌生的人了,這個小組就能代表之前的大組.而PCA的思想與之有些區別,PCA模型中的那個代表是另外構造的,並不是來自原先組中原本的特徵,如果我們將每個特徵看做一個維度的話,那麼構造出的代表其實就是將原先的多維變成少量新的維度.也就是說PCA的思想是將n維特徵映射到k維上（k<n），這k維是全新的正交特徵。這k維特徵稱為主元，是重新構造出來的k維特徵，而不是簡單地從n維特徵中去除其餘n-k維特徵。
三、主成分滿足的條件：1）每個主成分P都是原變數的線性組合，有多少個原變數就有多少個主成分，任意主成分可以表示成：
2）公式中的未知係數aij滿足平方和為1；3）P1是線性組合中方差最大，依次是P2，P3，...Pm，並且各主成分之間互不相關。四、主成分計算過程及數據模型接下來我們來看看主成分分析基本步驟：1.將原始數據標準化，用scale()函數2.求標準化數據的協方差陣，用cov()函數：或者求數據的相關陣用cor()函數3.求協方差陣或者相關矩陣的特徵值和單位特徵向量，用eigen()函數，其中$values是按從達到小對應的特徵值，$vectors是對應的單位特徵向量4.主成分分析，用princomp(x,cor...)函數,x為矩陣,cor為確定x是否為相關係數矩陣5.確定主成分個數,可以用screeplot()函數，用可視化的方法來確定主成分個數，選取一個拐彎點對應的序號6.解釋主成分，用PCA$loadings顯示主成分載荷矩陣，PCA為主成分分析賦值的變數。7.確定各樣本的主成分得分，用PCA$scores 來確定，並根據樣本各主成分的分值來對樣本進行解釋。對於一個樣本資料，觀測p 個變數（即特徵）x1, x2,……,xp， n 個樣品（觀測）的數據資料陣為：
主成分分析就是將p 個觀測變數綜合成為p 個新的變數（綜合變數），即：
要求模型滿足以下條件：
於是，稱F1為第一主成分，F2為第二主成分，依此類推，有第p 個主成分。主成分又叫主分量。這裡aij我們稱為主成分係數。
五、主成分案例：要用R語言結合一個小案例來實戰解析主成分分析的運用及結果分析。1. 函數總結#R中作為主成分分析最主要的函數是princomp()函數#princomp()主成分分析 可以從相關陣或者從協方差陣做主成分分析#summary()提取主成分信息#loadings()顯示主成分分析或因子分析中載荷的內容#predict()預測主成分的值#screeplot()畫出主成分的碎石圖#biplot()畫出數據關於主成分的散點圖和原坐標在主成分下的方向2. 案例1#現有30名中學生身高、體重、胸圍、坐高數據，對身體的四項指標數據做主成分分析。#1）載入原始數據test<-data.frame(X1=c(148, 139, 160, 149, 159, 142, 153, 150, 151, 139,140, 161, 158, 140, 137, 152, 149, 145, 160, 156,151, 147, 157, 147, 157, 151, 144, 141, 139, 148),X2=c(41, 34, 49, 36, 45, 31, 43, 43, 42, 31,29, 47, 49, 33, 31, 35, 47, 35, 47, 44,42, 38, 39, 30, 48, 36, 36, 30, 32, 38),X3=c(72, 71, 77, 67, 80, 66, 76, 77, 77, 68,64, 78, 78, 67, 66, 73, 82, 70, 74, 78,73, 73, 68, 65, 80, 74, 68, 67, 68, 70),X4=c(78, 76, 86, 79, 86, 76, 83, 79, 80, 74,74, 84, 83, 77, 73, 79, 79, 77, 87, 85,82, 78, 80, 75, 88, 80, 76, 76, 73, 78))> test X1 X2 X3 X41 148 41 72 782 139 34 71 763 160 49 77 864 149 36 67 795 159 45 80 866 142 31 66 767 153 43 76 838 150 43 77 799 151 42 77 8010 139 31 68 7411 140 29 64 7412 161 47 78 8413 158 49 78 8314 140 33 67 7715 137 31 66 7316 152 35 73 7917 149 47 82 7918 145 35 70 7719 160 47 74 8720 156 44 78 8521 151 42 73 8222 147 38 73 7823 157 39 68 8024 147 30 65 7525 157 48 80 8826 151 36 74 8027 144 36 68 7628 141 30 67 7629 139 32 68 7330 148 38 70 78#2）作主成分分析並顯示分析結果test.pr<-princomp(test,cor=TRUE) #cor是邏輯變數 當cor=TRUE表示用樣本的相關矩陣R做主成分分析；當cor=FALSE表示用樣本的協方差陣S做主成分分析summary(test.pr,loadings=TRUE) #loading是邏輯變數 當loading=TRUE時表示顯示loading 的內容#loadings的輸出結果為載荷 是主成分對應於原始變數的係數即Q矩陣Importance of components: Comp.1 Comp.2 Comp.3 Comp.4Standard deviation 1.8817805 0.55980636 0.28179594 0.25711844Proportion of Variance 0.8852745 0.07834579 0.01985224 0.01652747Cumulative Proportion 0.8852745 0.96362029 0.98347253 1.00000000Loadings: Comp.1 Comp.2 Comp.3 Comp.4X1 0.497 0.543 -0.450 0.506X2 0.515 -0.210 -0.462 -0.691X3 0.481 -0.725 0.175 0.461X4 0.507 0.368 0.744 -0.232#3）分析結果含義#----Standard deviation 標準差 其平方為方差=特徵值#----Proportion of Variance 方差貢獻率#----Cumulative Proportion 方差累計貢獻率#由結果顯示 前兩個主成分的累計貢獻率已經達到96% 可以捨去另外兩個主成分達到降維的目的因此可以得到函數表達式：Z1=-0.497X"1-0.515X"2-0.481X"3-0.507X"4Z1= 0.543X"1-0.210X"2-0.725X"3-0.368X"4#4）畫主成分的碎石圖並預測screeplot(test.pr,type="lines")
由碎石圖可以看出 第二個主成分之後 圖線變化趨於平穩 因此可以選擇前兩個主成分做分析！> p<-predict(test.pr) # predict函數根據主成分生成最後的新變數(即將原始數據矩陣X與主成分係數矩陣A結合,生成最後的主成分新變數,用這個生成的新變數代替原始變數！)> p Comp.1 Comp.2 Comp.3 Comp.4 [1,] -0.06990950 -0.23813701 -0.35509248 -0.266120139 [2,] -1.59526340 -0.71847399 0.32813232 -0.118056646 [3,] 2.84793151 0.38956679 -0.09731731 -0.279482487 [4,] -0.75996988 0.80604335 -0.04945722 -0.162949298 [5,] 2.73966777 0.01718087 0.36012615 0.358653044 [6,] -2.10583168 0.32284393 0.18600422 -0.036456084 [7,] 1.42105591 -0.06053165 0.21093321 -0.044223092 [8,] 0.82583977 -0.78102576 -0.27557798 0.057288572 [9,] 0.93464402 -0.58469242 -0.08814136 0.181037746[10,] -2.36463820 -0.36532199 0.08840476 0.045520127[11,] -2.83741916 0.34875841 0.03310423 -0.031146930[12,] 2.60851224 0.21278728 -0.33398037 0.210157574[13,] 2.44253342 -0.16769496 -0.46918095 -0.162987830[14,] -1.86630669 0.05021384 0.37720280 -0.358821916[15,] -2.81347421 -0.31790107 -0.03291329 -0.222035112[16,] -0.06392983 0.20718448 0.04334340 0.703533624[17,] 1.55561022 -1.70439674 -0.33126406 0.007551879[18,] -1.07392251 -0.06763418 0.02283648 0.048606680[19,] 2.52174212 0.97274301 0.12164633 -0.390667991[20,] 2.14072377 0.02217881 0.37410972 0.129548960[21,] 0.79624422 0.16307887 0.12781270 -0.294140762[22,] -0.28708321 -0.35744666 -0.03962116 0.080991989[23,] 0.25151075 1.25555188 -0.55617325 0.109068939[24,] -2.05706032 0.78894494 -0.26552109 0.388088643[25,] 3.08596855 -0.05775318 0.62110421 -0.218939612[26,] 0.16367555 0.04317932 0.24481850 0.560248997[27,] -1.37265053 0.02220972 -0.23378320 -0.257399715[28,] -2.16097778 0.13733233 0.35589739 0.093123683[29,] -2.40434827 -0.48613137 -0.16154441 -0.007914021[30,] -0.50287468 0.14734317 -0.20590831 -0.1220788193. 案例2：變數分類問題#案例：已知16項指標的相關矩陣 從相關矩陣R出發 進行主成分分析 對16個指標進行分類#1）載入相關係數 下三角表示（輸入R之後，其實是以向量存儲的）x<-c(1.00,0.79, 1.00,0.36, 0.31, 1.00,0.96, 0.74, 0.38, 1.00,0.89, 0.58, 0.31, 0.90, 1.00,0.79, 0.58, 0.30, 0.78, 0.79, 1.00,0.76, 0.55, 0.35, 0.75, 0.74, 0.73, 1.00,0.26, 0.19, 0.58, 0.25, 0.25, 0.18, 0.24, 1.00,0.21, 0.07, 0.28, 0.20, 0.18, 0.18, 0.29,-0.04, 1.00,0.26, 0.16, 0.33, 0.22, 0.23, 0.23, 0.25, 0.49,-0.34, 1.00,0.07, 0.21, 0.38, 0.08,-0.02, 0.00, 0.10, 0.44,-0.16, 0.23, 1.00,0.52, 0.41, 0.35, 0.53, 0.48, 0.38, 0.44, 0.30,-0.05, 0.50,0.24, 1.00,0.77, 0.47, 0.41, 0.79, 0.79, 0.69, 0.67, 0.32, 0.23, 0.31,0.10, 0.62, 1.00,0.25, 0.17, 0.64, 0.27, 0.27, 0.14, 0.16, 0.51, 0.21, 0.15,0.31, 0.17, 0.26, 1.00,0.51, 0.35, 0.58, 0.57, 0.51, 0.26, 0.38, 0.51, 0.15, 0.29,0.28, 0.41, 0.50, 0.63, 1.00,0.21, 0.16, 0.51, 0.26, 0.23, 0.00, 0.12, 0.38, 0.18, 0.14,0.31, 0.18, 0.24, 0.50, 0.65, 1.00))> x[1] 1.00 0.79 1.00 0.36 0.31 1.00 0.96 0.74 0.38 1.00 0.89 0.58 0.31 0.90 1.00 0.79[17] 0.58 0.30 0.78 0.79 1.00 0.76 0.55 0.35 0.75 0.74 0.73 1.00 0.26 0.19 0.58 0.25[33] 0.25 0.18 0.24 1.00 0.21 0.07 0.28 0.20 0.18 0.18 0.29 -0.04 1.00 0.26 0.16 0.33[49] 0.22 0.23 0.23 0.25 0.49 -0.34 1.00 0.07 0.21 0.38 0.08 -0.02 0.00 0.10 0.44 -0.16[65] 0.23 1.00 0.52 0.41 0.35 0.53 0.48 0.38 0.44 0.30 -0.05 0.50 0.24 1.00 0.77 0.47[81] 0.41 0.79 0.79 0.69 0.67 0.32 0.23 0.31 0.10 0.62 1.00 0.25 0.17 0.64 0.27 0.27[97] 0.14 0.16 0.51 0.21 0.15 0.31 0.17 0.26 1.00 0.51 0.35 0.58 0.57 0.51 0.26 0.38[113] 0.51 0.15 0.29 0.28 0.41 0.50 0.63 1.00 0.21 0.16 0.51 0.26 0.23 0.00 0.12 0.38[129] 0.18 0.14 0.31 0.18 0.24 0.50 0.65 1.00#2）輸入16個指標的名稱names<-c("X1", "X2", "X3", "X4", "X5", "X6", "X7", "X8","X9","X10", "X11", "X12", "X13", "X14", "X15", "X16")#3）將矩陣生成相關矩陣 即原始下三角矩陣的擴展R<-matrix(0,nrow=16,ncol=16,dimnames=list(names,names))<-for(i in 1:16){for(j in 1:i){R[i,j]<-x[(i-1)*i / 2+j]; R[j,i]<-R[i,j] #這個公式略D，有木有！} }> RX1 X2 X3 X4 X5 X6 X7 X8 X9 X10 X11 X12 X13 X14 X15 X16X1 1.00 0.79 0.36 0.96 0.89 0.79 0.76 0.26 0.21 0.26 0.07 0.52 0.77 0.25 0.51 0.21X2 0.79 1.00 0.31 0.74 0.58 0.58 0.55 0.19 0.07 0.16 0.21 0.41 0.47 0.17 0.35 0.16X3 0.36 0.31 1.00 0.38 0.31 0.30 0.35 0.58 0.28 0.33 0.38 0.35 0.41 0.64 0.58 0.51X4 0.96 0.74 0.38 1.00 0.90 0.78 0.75 0.25 0.20 0.22 0.08 0.53 0.79 0.27 0.57 0.26X5 0.89 0.58 0.31 0.90 1.00 0.79 0.74 0.25 0.18 0.23 -0.02 0.48 0.79 0.27 0.51 0.23X6 0.79 0.58 0.30 0.78 0.79 1.00 0.73 0.18 0.18 0.23 0.00 0.38 0.69 0.14 0.26 0.00X7 0.76 0.55 0.35 0.75 0.74 0.73 1.00 0.24 0.29 0.25 0.10 0.44 0.67 0.16 0.38 0.12X8 0.26 0.19 0.58 0.25 0.25 0.18 0.24 1.00 -0.04 0.49 0.44 0.30 0.32 0.51 0.51 0.38X9 0.21 0.07 0.28 0.20 0.18 0.18 0.29 -0.04 1.00 -0.34 -0.16 -0.05 0.23 0.21 0.15 0.18X10 0.26 0.16 0.33 0.22 0.23 0.23 0.25 0.49 -0.34 1.00 0.23 0.50 0.31 0.15 0.29 0.14X11 0.07 0.21 0.38 0.08 -0.02 0.00 0.10 0.44 -0.16 0.23 1.00 0.24 0.10 0.31 0.28 0.31X12 0.52 0.41 0.35 0.53 0.48 0.38 0.44 0.30 -0.05 0.50 0.24 1.00 0.62 0.17 0.41 0.18X13 0.77 0.47 0.41 0.79 0.79 0.69 0.67 0.32 0.23 0.31 0.10 0.62 1.00 0.26 0.50 0.24X14 0.25 0.17 0.64 0.27 0.27 0.14 0.16 0.51 0.21 0.15 0.31 0.17 0.26 1.00 0.63 0.50X15 0.51 0.35 0.58 0.57 0.51 0.26 0.38 0.51 0.15 0.29 0.28 0.41 0.50 0.63 1.00 0.65X16 0.21 0.16 0.51 0.26 0.23 0.00 0.12 0.38 0.18 0.14 0.31 0.18 0.24 0.50 0.65 1.00> x<-c(1,-1,1,-2,-3,1)> names<-c("x1","x2","x3")> R<-matrix(0,nrow=3,ncol=3,dimnames=list(names,names))> R x1 x2 x3x1 0 0 0x2 0 0 0x3 0 0 0> for(i in 1:3){ + for(j in 1:i){ + R[i,j]<-x[(i-1)*i/2+j];R[j,i]<-R[i,j]#可見該公式可將相關係數下三角轉化成相關矩陣+ } + }> R x1 x2 x3x1 1 -1 -2x2 -1 1 -3x3 -2 -3 1#4）主成分分析> x.pr<-princomp(covmat=R) #covmat協方差陣（指定使用相關矩陣來進行PCA分析！）> x.prCall:princomp(covmat = R)Standard deviations:Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7 Comp.8 Comp.9 Comp.10 Comp.11 Comp.12 Comp.13 Comp.14 Comp.15 Comp.16 2.6526359 1.6167971 1.2775386 0.9217124 0.8763498 0.8037803 0.6873924 0.6753023 0.5975465 0.5559751 0.4901757 0.4666029 0.4112178 0.3675549 0.2607323 0.169998216 variables and NA observations.> load<-loadings(x.pr) #載荷
#5）散點圖plot(load[,1:2]) #取前兩個主成分,畫兩個主成分的散點圖！text(load[,1],load[,2],adj=c(-0.4,0.3))
通過散點圖可以得到分析結果：指標1 2 4 5 6 7 13可以歸為一類 ；指標9 10 12可以歸為一類；指標3 8 11 14 15 16可以歸為一類。附錄特徵值和特徵向量隱藏的秘密：主成分變數對應的特徵向量的每個元素，與對應的特徵值的平方根的乘積，等於該主成分變數，與該元素列標籤對應的原始變數之間的相關係數。這是特徵值與特徵向量隱藏的秘密，可以用矩陣代數嚴格推導出來。不過這句話讀起來比較費勁，我們用圖8來表示這一關係。圖中的eigVec1至eigVec4是4個特徵向量，對應的特徵值分別為eigVal1至eigVal4。我們在每個列中進行操作，用特徵向量每個元素分別乘以對應特徵值的平方根，得到該主成分變數與所有原始變數的相關係數！
4. 案例3：主成分回歸問題#考慮進口總額Y與三個自變數：國內總產值，存儲量，總消費量之間的關係。現收集了1949-1959共11年的數據，試做線性回歸和主成分回歸分析。economy<-data.frame(x1=c(149.3, 161.2, 171.5, 175.5, 180.8, 190.7, 202.1, 212.4, 226.1, 231.9, 239.0),x2=c(4.2, 4.1, 3.1, 3.1, 1.1, 2.2, 2.1, 5.6, 5.0, 5.1, 0.7),x3=c(108.1, 114.8, 123.2, 126.9, 132.1, 137.7, 146.0, 154.1, 162.3, 164.3, 167.6),y=c(15.9, 16.4, 19.0, 19.1, 18.8, 20.4, 22.7, 26.5, 28.1, 27.6, 26.3))> economy x1 x2 x3 y1 149.3 4.2 108.1 15.92 161.2 4.1 114.8 16.43 171.5 3.1 123.2 19.04 175.5 3.1 126.9 19.15 180.8 1.1 132.1 18.86 190.7 2.2 137.7 20.47 202.1 2.1 146.0 22.78 212.4 5.6 154.1 26.59 226.1 5.0 162.3 28.110 231.9 5.1 164.3 27.611 239.0 0.7 167.6 26.3#1）線性回歸lm.sol<-lm(y~x1+x2+x3, data=economy)summary(lm.sol)@結果> lm.sol<-lm(y~x1+x2+x3, data=economy)> summary(lm.sol)Call:lm(formula = y ~ x1 + x2 + x3, data = economy)Residuals: Min 1Q Median 3Q Max -0.52367 -0.38953 0.05424 0.22644 0.78313 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -10.12799 1.21216 -8.355 6.9e-05 ***x1 -0.05140 0.07028 -0.731 0.488344 x2 0.58695 0.09462 6.203 0.000444 ***x3 0.28685 0.10221 2.807 0.026277 * ---Signif. codes: 0 『***』 0.001 『**』 0.01 『*』 0.05 『.』 0.1 『 』 1Residual standard error: 0.4889 on 7 degrees of freedomMultiple R-squared: 0.9919, Adjusted R-squared: 0.9884 F-statistic: 285.6 on 3 and 7 DF, p-value: 1.112e-07###2）主成分回歸# 主成分分析economy.pr<-princomp(~x1+x2+x3, data=economy, cor=T)summary(economy.pr, loadings=TRUE)Importance of components: Comp.1 Comp.2 Comp.3Standard deviation 1.413915 0.9990767 0.0518737839Proportion of Variance 0.666385 0.3327181 0.0008969632Cumulative Proportion 0.666385 0.9991030 1.0000000000Loadings: Comp.1 Comp.2 Comp.3x1 -0.706 0.707x2 -0.999 x3 -0.707 -0.707pre<-predict(economy.pr)economy$z1<-pre[,1]; economy$z2<-pre[,2]> economy x1 x2 x3 y z1 z21 149.3 4.2 108.1 15.9 2.2296493 -0.669830322 161.2 4.1 114.8 16.4 1.6979452 -0.582654453 171.5 3.1 123.2 19.0 1.1695976 0.076541754 175.5 3.1 126.9 19.1 0.9379462 0.086390365 180.8 1.1 132.1 18.8 0.6756511 1.370463036 190.7 2.2 137.7 20.4 0.1996423 0.691319687 202.1 2.1 146.0 22.7 -0.3771746 0.779972368 212.4 5.6 154.1 26.5 -1.0192344 -1.420148829 226.1 5.0 162.3 28.1 -1.6354243 -1.0110995310 231.9 5.1 164.3 27.6 -1.8532401 -1.0647686411 239.0 0.7 167.6 26.3 -2.0253583 1.74381457lm.sol<-lm(y~z1+z2, data=economy)summary(lm.sol)Call:lm(formula = y ~ z1 + z2, data = economy)Residuals: Min 1Q Median 3Q Max -0.89838 -0.26050 0.08435 0.35677 0.66863 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 21.8909 0.1658 132.006 1.21e-14 ***z1 -2.9892 0.1173 -25.486 6.02e-09 ***z2 -0.8288 0.1660 -4.993 0.00106 ** ---Signif. codes: 0 『***』 0.001 『**』 0.01 『*』 0.05 『.』 0.1 『 』 1Residual standard error: 0.55 on 8 degrees of freedomMultiple R-squared: 0.9883, Adjusted R-squared: 0.9853 F-statistic: 337.2 on 2 and 8 DF, p-value: 1.888e-08
推薦閱讀：
