R(R言語)で、主成分分析(寄与率、寄与率の求め方)

Rコード

#
# R(R言語)で、主成分分析(寄与率、寄与率の求め方)
#




#
# 使用するデータと解説は、『多変量解析がわかる』 第3章 p.81~82
#




# 使用するデータ 20人の中学生の5教科のテスト結果
出席番号_vec <- c(1:20)

数学x_vec <- c(71, 34, 58, 41, 69, 64, 16, 59, 57, 46,
            23, 39, 46, 52, 39, 23, 37, 52, 63, 39)

理科y_vec <- c(64, 48, 59, 51, 56, 65, 45, 59, 54, 54,
            49, 48, 55, 56, 53, 43, 45, 51, 56, 49)

社会u_vec <- c(83, 67, 78, 70, 74, 82, 63, 78, 84, 71,
            64, 71, 68, 82, 78, 63, 67, 74, 79, 73)

英語v_vec <- c(100, 57, 87, 60, 81, 100, 7, 59, 73, 43,
            33, 29, 42, 67, 52, 35, 39, 65, 91, 64)

国語w_vec <- c(71, 68, 66, 72, 66, 71, 59, 62, 72, 62,
            70, 66, 61, 60, 72, 59, 70, 69, 70, 60)

df1 <- data.frame(出席番号 = 出席番号_vec,
                数学x = 数学x_vec,
                理科y = 理科y_vec,
                社会u = 社会u_vec,
                英語v = 英語v_vec,
                国語w = 国語w_vec)

df1




summary(df1)




# 分散 (不偏分散ではない)
分散 <- function(data) {
    sum((data - mean(data))^2) / length(data)
}
# 全分散の合計 全分散 = sx^2 + sy^2 + su^2 + sv^2 + sw^2
sum(分散(df1$数学x), 分散(df1$理科y), 分散(df1$社会u), 分散(df1$英語v), 分散(df1$国語w))

sx2_sy2_su2_sv2_sw2 <- round(sum(分散(df1$数学x), 分散(df1$理科y), 分散(df1$社会u),
                                 分散(df1$英語v), 分散(df1$国語w)), 1)
sx2_sy2_su2_sv2_sw2

# 主成分Pの分散 sp^2
# 寄与率 C = sp^2 / (sx^2 + sy^2 + su^2 + sv^2 + sw^2)
# 0 <= C <= 1


df2 <- df1[, -1] # 出席番号を削除
df2


# prcomp()を使った場合
res1 <- prcomp(df2)
res1
summary(res1)
res1$rotation # 主成分
res1$rotation[, 1] # 第1主成分
sum(res1$rotation[, 1]^2) # 第1主成分(書籍では主成分負荷量)の平方和
pc <- round(abs(res1$rotation[, 1]), 2) # 第1主成分を、絶対値を取り変数に付値(代入)
pc


# 20名全員の主成分得点(合成変量)
df_mat <- as.matrix(df2)
df_mat
x <- df_mat[1:20, ] %*% pc
分散(round(apply(x, 1, sum), 1))
sp2 <- round(分散(round(apply(x, 1, sum), 1)), 1)
sp2 # 書籍では数値は844.5


# 寄与率
C <- sp2 / sx2_sy2_su2_sv2_sw2
C
round(C, 2) # 書籍では数値は0.91




# princomp()を使った場合
res2 <- princomp(df2)
res2
summary(res2) # Proportion of Varianceが寄与率、Cumulative Proportionが累積寄与率
res2$loadings # 主成分
res2$score # 主成分得点
par(family = "Osaka") # 文字化けするなら
biplot(res2)

R Console

> #
> # R(R言語)で、主成分分析(寄与率、寄与率の求め方)
> #
> 
> 
> 
> 
> #
> # 使用するデータと解説は、『多変量解析がわかる』 第3章 p.81~82
> #
> 
> 
> 
> 
> # 使用するデータ 20人の中学生の5教科のテスト結果
> 出席番号_vec <- c(1:20)
> 
> 数学x_vec <- c(71, 34, 58, 41, 69, 64, 16, 59, 57, 46,
+             23, 39, 46, 52, 39, 23, 37, 52, 63, 39)
> 
> 理科y_vec <- c(64, 48, 59, 51, 56, 65, 45, 59, 54, 54,
+             49, 48, 55, 56, 53, 43, 45, 51, 56, 49)
> 
> 社会u_vec <- c(83, 67, 78, 70, 74, 82, 63, 78, 84, 71,
+             64, 71, 68, 82, 78, 63, 67, 74, 79, 73)
> 
> 英語v_vec <- c(100, 57, 87, 60, 81, 100, 7, 59, 73, 43,
+             33, 29, 42, 67, 52, 35, 39, 65, 91, 64)
> 
> 国語w_vec <- c(71, 68, 66, 72, 66, 71, 59, 62, 72, 62,
+             70, 66, 61, 60, 72, 59, 70, 69, 70, 60)
> 
> df1 <- data.frame(出席番号 = 出席番号_vec,
+                 数学x = 数学x_vec,
+                 理科y = 理科y_vec,
+                 社会u = 社会u_vec,
+                 英語v = 英語v_vec,
+                 国語w = 国語w_vec)
> 
> df1
   出席番号 数学x 理科y 社会u 英語v 国語w
1         1    71    64    83   100    71
2         2    34    48    67    57    68
3         3    58    59    78    87    66
4         4    41    51    70    60    72
5         5    69    56    74    81    66
6         6    64    65    82   100    71
7         7    16    45    63     7    59
8         8    59    59    78    59    62
9         9    57    54    84    73    72
10       10    46    54    71    43    62
11       11    23    49    64    33    70
12       12    39    48    71    29    66
13       13    46    55    68    42    61
14       14    52    56    82    67    60
15       15    39    53    78    52    72
16       16    23    43    63    35    59
17       17    37    45    67    39    70
18       18    52    51    74    65    69
19       19    63    56    79    91    70
20       20    39    49    73    64    60
> 
> 
> 
> 
> summary(df1)
    出席番号         数学x           理科y           社会u           英語v            国語w      
 Min.   : 1.00   Min.   :16.00   Min.   :43.00   Min.   :63.00   Min.   :  7.00   Min.   :59.00  
 1st Qu.: 5.75   1st Qu.:38.50   1st Qu.:48.75   1st Qu.:67.75   1st Qu.: 41.25   1st Qu.:61.75  
 Median :10.50   Median :46.00   Median :53.50   Median :73.50   Median : 59.50   Median :67.00  
 Mean   :10.50   Mean   :46.40   Mean   :53.00   Mean   :73.45   Mean   : 59.20   Mean   :66.30  
 3rd Qu.:15.25   3rd Qu.:58.25   3rd Qu.:56.00   3rd Qu.:78.25   3rd Qu.: 75.00   3rd Qu.:70.25  
 Max.   :20.00   Max.   :71.00   Max.   :65.00   Max.   :84.00   Max.   :100.00   Max.   :72.00  
> 
> 
> 
> 
> # 分散 (不偏分散ではない)
> 分散 <- function(data) {
+     sum((data - mean(data))^2) / length(data)
+ }
> # 全分散の合計 全分散 = sx^2 + sy^2 + su^2 + sv^2 + sw^2
> sum(分散(df1$数学x), 分散(df1$理科y), 分散(df1$社会u), 分散(df1$英語v), 分散(df1$国語w))
[1] 923.6575
> 
> sx2_sy2_su2_sv2_sw2 <- round(sum(分散(df1$数学x), 分散(df1$理科y), 分散(df1$社会u),
+                                  分散(df1$英語v), 分散(df1$国語w)), 1)
> sx2_sy2_su2_sv2_sw2
[1] 923.7
> 
> # 主成分Pの分散 sp^2
> # 寄与率 C = sp^2 / (sx^2 + sy^2 + su^2 + sv^2 + sw^2)
> # 0 <= C <= 1
> 
> 
> df2 <- df1[, -1] # 出席番号を削除
> df2
   数学x 理科y 社会u 英語v 国語w
1     71    64    83   100    71
2     34    48    67    57    68
3     58    59    78    87    66
4     41    51    70    60    72
5     69    56    74    81    66
6     64    65    82   100    71
7     16    45    63     7    59
8     59    59    78    59    62
9     57    54    84    73    72
10    46    54    71    43    62
11    23    49    64    33    70
12    39    48    71    29    66
13    46    55    68    42    61
14    52    56    82    67    60
15    39    53    78    52    72
16    23    43    63    35    59
17    37    45    67    39    70
18    52    51    74    65    69
19    63    56    79    91    70
20    39    49    73    64    60
> 
> 
> # prcomp()を使った場合
> res1 <- prcomp(df2)
> res1
Standard deviations:
[1] 29.814345  6.796644  4.231235  3.517360  2.627993

Rotation:
              PC1        PC2         PC3         PC4        PC5
数学x -0.49151300  0.7559881  0.10647703 -0.39445342  0.1413013
理科y -0.17305310  0.2378355 -0.01310165  0.34925878 -0.8895693
社会u -0.19588597  0.1986404  0.22761783  0.83516616  0.4157623
英語v -0.82781535 -0.5257361 -0.19384988 -0.02333993  0.0141702
国語w -0.06941201 -0.2367695  0.94821060 -0.15613309 -0.1250652
> summary(res1)
Importance of components:
                           PC1     PC2     PC3     PC4    PC5
Standard deviation     29.8143 6.79664 4.23124 3.51736 2.6280
Proportion of Variance  0.9143 0.04751 0.01841 0.01272 0.0071
Cumulative Proportion   0.9143 0.96176 0.98017 0.99290 1.0000
> res1$rotation # 主成分
              PC1        PC2         PC3         PC4        PC5
数学x -0.49151300  0.7559881  0.10647703 -0.39445342  0.1413013
理科y -0.17305310  0.2378355 -0.01310165  0.34925878 -0.8895693
社会u -0.19588597  0.1986404  0.22761783  0.83516616  0.4157623
英語v -0.82781535 -0.5257361 -0.19384988 -0.02333993  0.0141702
国語w -0.06941201 -0.2367695  0.94821060 -0.15613309 -0.1250652
> res1$rotation[, 1] # 第1主成分
      数学x       理科y       社会u       英語v       国語w 
-0.49151300 -0.17305310 -0.19588597 -0.82781535 -0.06941201 
> sum(res1$rotation[, 1]^2) # 第1主成分(書籍では主成分負荷量)の平方和
[1] 1
> pc <- round(abs(res1$rotation[, 1]), 2) # 第1主成分を、絶対値を取り変数に付値(代入)
> pc
数学x 理科y 社会u 英語v 国語w 
 0.49  0.17  0.20  0.83  0.07 
> 
> 
> # 20名全員の主成分得点(合成変量)
> df_mat <- as.matrix(df2)
> df_mat
      数学x 理科y 社会u 英語v 国語w
 [1,]    71    64    83   100    71
 [2,]    34    48    67    57    68
 [3,]    58    59    78    87    66
 [4,]    41    51    70    60    72
 [5,]    69    56    74    81    66
 [6,]    64    65    82   100    71
 [7,]    16    45    63     7    59
 [8,]    59    59    78    59    62
 [9,]    57    54    84    73    72
[10,]    46    54    71    43    62
[11,]    23    49    64    33    70
[12,]    39    48    71    29    66
[13,]    46    55    68    42    61
[14,]    52    56    82    67    60
[15,]    39    53    78    52    72
[16,]    23    43    63    35    59
[17,]    37    45    67    39    70
[18,]    52    51    74    65    69
[19,]    63    56    79    91    70
[20,]    39    49    73    64    60
> x <- df_mat[1:20, ] %*% pc
> 分散(round(apply(x, 1, sum), 1))
[1] 846.3133
> sp2 <- round(分散(round(apply(x, 1, sum), 1)), 1)
> sp2 # 書籍では数値は844.5
[1] 846.3
> 
> 
> # 寄与率
> C <- sp2 / sx2_sy2_su2_sv2_sw2
> C
[1] 0.9162066
> round(C, 2) # 書籍では数値は0.91
[1] 0.92
> 
> 
> 
> 
> # princomp()を使った場合
> res2 <- princomp(df2)
> res2
Call:
princomp(x = df2)

Standard deviations:
   Comp.1    Comp.2    Comp.3    Comp.4    Comp.5 
29.059429  6.624549  4.124098  3.428298  2.561450 

 5  variables and  20 observations.
> summary(res2) # Proportion of Varianceが寄与率、Cumulative Proportionが累積寄与率
Importance of components:
                           Comp.1     Comp.2     Comp.3     Comp.4      Comp.5
Standard deviation     29.0594289 6.62454895 4.12409778 3.42829844 2.561450450
Proportion of Variance  0.9142463 0.04751182 0.01841395 0.01272466 0.007103313
Cumulative Proportion   0.9142463 0.96175807 0.98017202 0.99289669 1.000000000
> res2$loadings # 主成分

Loadings:
      Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
数学x -0.492  0.756  0.106  0.394 -0.141
理科y -0.173  0.238        -0.349  0.890
社会u -0.196  0.199  0.228 -0.835 -0.416
英語v -0.828 -0.526 -0.194              
国語w        -0.237  0.948  0.156  0.125

               Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
SS loadings       1.0    1.0    1.0    1.0    1.0
Proportion Var    0.2    0.2    0.2    0.2    0.2
Cumulative Var    0.2    0.4    0.6    0.8    1.0
> res2$score # 主成分得点
           Comp.1      Comp.2     Comp.3      Comp.4     Comp.5
 [1,] -49.9666177   0.5476662  1.1964818 -0.42803473  2.3483821
 [2,]   9.9266846 -11.0905498 -0.6845141  2.45597161  0.2297421
 [3,] -30.6235937  -3.4441428 -3.4813050 -0.71788899  1.3751507
 [4,]   2.6181823  -7.0134912  3.9156664  2.35842287  1.1198042
 [5,] -29.7606415   6.5180747 -2.0181247  7.86950005 -1.0998014
 [6,] -46.5031938  -4.7050556  0.2104231 -2.70330130  4.6428230
 [7,]  62.0920974   2.2113247 -2.3136686 -2.82794310  1.3504303
 [8,]  -7.6586289  12.9795333 -1.7398738 -1.60148586  1.1303542
 [9,] -19.2691882   1.7422224  6.2465950 -3.76700589 -4.4771938
[10,]  14.2125530   8.9838043 -1.5502937  0.48963785  1.6564842
[11,]  35.4766768  -7.6203780  3.9971016  0.02533161  4.5111283
[12,]  30.0032295   8.6781017  4.2897177  0.12178998 -1.9931785
[13,]  15.5243852   9.3882241 -3.0006096  2.46640453  3.6824453
[14,] -10.9661212   4.0363220 -4.9826570 -6.78109486 -2.5757857
[15,]   8.3105371  -2.2547843  7.0482506 -5.99705019 -0.0311912
[16,]  35.8187827  -7.6930397 -6.9699226  1.28526634 -1.8145834
[17,]  23.7331571  -0.5463815  5.0599409  4.57925572 -2.3576756
[18,]  -7.5028454   0.1785679  2.1835038  3.00504624 -2.5436063
[19,] -36.3467948  -3.2290905  0.3354458  2.18488035 -1.9722459
[20,]   0.8813393  -7.6669279 -7.7421576 -2.01770223 -3.1814824
> par(family = "Osaka") # 文字化けするなら
> biplot(res2)

biplotのスクリーンショット

f:id:my_notes:20170719201626p:plain

参考文献

多変量解析がわかる (ファーストブック)

多変量解析がわかる (ファーストブック)