R(R言語)で、主成分分析(主成分の求め方、主成分得点)

Rコード

#
# R(R言語)で、主成分分析(主成分の求め方、主成分得点)
#




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




# 使用するデータ 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)

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

df




summary(df)




# 分散
分散 <- function(data) {
    sum((data - mean(data))^2) / length(data)
}
分散(df$数学x)
分散(df$理科y)
分散(df$社会u)
分散(df$英語v)
分散(df$国語w)
# 分散が大きい変量(変数)が主成分に大きく寄与する。
# 英語の分散に注目。




# 合成変量p(主成分(principal component))の数式
# p = ax + by + cu + dv + ew
# ただし、a^2 + b^2 + c^2 + d^2 + e^2 = 1




# prcomp()を使った場合
res1 <- prcomp(df)
res1
summary(res1)
res1$rotation # 主成分
res1$rotation[, 1] # 第1主成分
sum(res1$rotation[, 1]^2) # 第1主成分(書籍では主成分負荷量)の平方和
pc <- round(abs(res1$rotation[, 1]), 2) # 第1主成分を、絶対値をとり変数に付値(代入)
# 2番目と10番目の生徒の主成分得点と単純な合計点と比べてみる
sum(df[2, ] * pc)
sum(df[2, ])
sum(df[10, ] * pc)
sum(df[10, ])
# 主成分得点が高いほうが、学力は上ではないかと考えられる。


# 20名全員の主成分得点(合成変量)
df_mat <- as.matrix(df)
df_mat
x <- df_mat[1:20, ] %*% pc
round(apply(x, 1, sum), 1)
# これを各人の単純な合計点と比べてみる。




# princomp()を使った場合
res2 <- princomp(df)
res2
summary(res2)
res2$loadings # 主成分
res2$score # 主成分得点

R Console

> #
> # R(R言語)で、主成分分析(主成分の求め方、主成分得点)
> #
> 
> 
> 
> 
> #
> # 使用するデータと解説は、『多変量解析がわかる』 第3章 p.77~80
> #
> 
> 
> 
> 
> # 使用するデータ 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)
> 
> df <- data.frame(数学x = 数学x_vec,
+                 理科y = 理科y_vec,
+                 社会u = 社会u_vec,
+                 英語v = 英語v_vec,
+                 国語w = 国語w_vec)
> 
> df
   数学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
> 
> 
> 
> 
> summary(df)
     数学x           理科y           社会u           英語v            国語w      
 Min.   :16.00   Min.   :43.00   Min.   :63.00   Min.   :  7.00   Min.   :59.00  
 1st Qu.:38.50   1st Qu.:48.75   1st Qu.:67.75   1st Qu.: 41.25   1st Qu.:61.75  
 Median :46.00   Median :53.50   Median :73.50   Median : 59.50   Median :67.00  
 Mean   :46.40   Mean   :53.00   Mean   :73.45   Mean   : 59.20   Mean   :66.30  
 3rd Qu.:58.25   3rd Qu.:56.00   3rd Qu.:78.25   3rd Qu.: 75.00   3rd Qu.:70.25  
 Max.   :71.00   Max.   :65.00   Max.   :84.00   Max.   :100.00   Max.   :72.00  
> 
> 
> 
> 
> # 分散
> 分散 <- function(data) {
+     sum((data - mean(data))^2) / length(data)
+ }
> 分散(df$数学x)
[1] 231.24
> 分散(df$理科y)
[1] 34.4
> 分散(df$社会u)
[1] 44.3475
> 分散(df$英語v)
[1] 591.46
> 分散(df$国語w)
[1] 22.21
> # 分散が大きい変量(変数)が主成分に大きく寄与する。
> # 英語の分散に注目。
> 
> 
> 
> 
> # 合成変量p(主成分(principal component))の数式
> # p = ax + by + cu + dv + ew
> # ただし、a^2 + b^2 + c^2 + d^2 + e^2 = 1
> 
> 
> 
> 
> # prcomp()を使った場合
> res1 <- prcomp(df)
> 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主成分を、絶対値をとり変数に付値(代入)
> # 2番目と10番目の生徒の主成分得点と単純な合計点と比べてみる
> sum(df[2, ] * pc)
[1] 90.29
> sum(df[2, ])
[1] 274
> sum(df[10, ] * pc)
[1] 85.95
> sum(df[10, ])
[1] 276
> # 主成分得点が高いほうが、学力は上ではないかと考えられる。
> 
> 
> # 20名全員の主成分得点(合成変量)
> df_mat <- as.matrix(df)
> 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] 150.2  90.3 130.9  97.6 130.0 146.8  38.0 107.8 119.5  86.0  64.7  70.2  84.6 111.2  91.9  64.4
[17]  76.5 107.7 136.6  99.4
> # これを各人の単純な合計点と比べてみる。
> 
> 
> 
> 
> # princomp()を使った場合
> res2 <- princomp(df)
> res2
Call:
princomp(x = df)

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)
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

参考文献

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

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