統計ソフトRでコレスポンデンス分析はどうやる?

にほんブログ村 科学ブログ 数学へ

コレスポンデンス分析とは、 人の属性の類似性(もしくは対応)を図に表した、 結果を見るといろいろと解釈したくなる、 面白い分析。

統計ソフトRでは、caというライブラリで分析可能。

最初の一回だけインストールをする。

install.packages("ca")

caライブラリを呼び出す。

library(ca)

例題のデータを読み込む。

こちらから、CSVファイルをダウンロードできる。

dat <- read.csv("C:/data/secondhand-smoke-place-ca-data.csv")

2番さんから30番さんに 受動喫煙の被害にあった場所を 回答してもらったデータ。

27人のデータ。14番と28番が無回答。

> table(dat$NO, dat$basho)
    
     カラオケ バー パチンコ レストラン 駅 家 学校 居酒屋 車内 職場 食堂
  2         1    0        1          1  0  1    1      1    0    1    1
  3         0    0        0          0  0  0    0      0    1    0    0
  4         0    0        0          1  1  0    0      1    0    1    0
  5         1    1        0          1  0  0    0      1    1    1    1
  6         1    0        0          1  0  0    0      1    1    0    0
  7         1    1        1          1  1  1    0      1    1    1    1
  8         1    1        1          1  0  0    0      1    1    1    1
  9         1    1        1          1  0  1    0      1    1    1    1
  10        0    0        0          1  1  1    0      1    0    0    1
  11        1    1        1          1  0  1    0      1    1    1    0
  12        0    0        0          1  0  1    0      0    0    0    0
  13        0    0        0          1  1  0    0      0    0    0    0
  15        1    0        0          1  0  0    0      1    1    1    0
  16        1    1        1          1  1  0    0      1    0    1    0
  17        0    0        1          0  0  0    0      1    0    0    0
  18        0    0        0          0  0  1    1      0    0    1    0
  19        1    0        0          1  1  1    0      1    1    1    0
  20        1    0        0          1  1  1    0      1    0    1    1
  21        1    1        0          1  0  1    0      1    1    0    0
  22        0    0        0          0  0  0    0      0    1    1    0
  23        0    0        0          1  0  1    0      0    0    1    0
  24        1    1        0          1  0  0    0      1    0    0    0
  25        0    0        0          0  0  0    1      0    0    1    0
  26        0    0        0          1  1  0    0      0    0    1    0
  27        0    0        1          1  0  1    0      0    0    0    0
  29        1    1        0          1  0  1    0      1    1    0    0
  30        1    1        1          1  1  0    0      1    1    1    1
> 

table()でオブジェクトを作り、 ca()をかければよい。

tab.NO.basho <- table(dat$NO, dat$basho)

ca.NO.basho <- ca(tab.NO.basho)

summary(ca.NO.basho)

plot(ca.NO.basho)

summary()で結果の概要が見られる。 Principal Inertia(固有値)1次元目と2次元目を使って、 図に表すと属性の広がりと、個人の散らばりが 視覚化できる。plot()を使う。

> summary(ca.NO.basho)

Principal inertias (eigenvalues):

 dim    value      %   cum%   scree plot               
 1      0.370047  30.3  30.3  ********                 
 2      0.225220  18.5  48.8  *****                    
 3      0.171709  14.1  62.9  ****                     
 4      0.142078  11.6  74.5  ***                      
 5      0.096621   7.9  82.4  **                       
 6      0.075889   6.2  88.7  **                       
 7      0.054686   4.5  93.2  *                        
 8      0.043836   3.6  96.7  *                        
 9      0.025404   2.1  98.8  *                        
 10     0.014295   1.2 100.0                           
        -------- -----                                 
 Total: 1.219786 100.0                                 


Rows:
     name   mass  qlt  inr    k=1 cor ctr     k=2 cor ctr  
1  |    2 |   58  620   42 |  729 612  84 |    83   8   2 |
2  |    3 |    7  344   57 | -588  36   7 |  1712 307  95 |
3  |    4 |   29  536   32 |  -25   0   0 |  -841 535  92 |
4  |    5 |   51  348   20 | -223 104   7 |   341 243  26 |
5  |    6 |   29  341   26 | -420 161  14 |   443 179  25 |
6  |    7 |   73  198    8 | -157 185   5 |   -42  13   1 |
7  |    8 |   58  447   18 | -220 131   8 |   342 317  30 |
8  |    9 |   66  397   11 | -141  94   4 |   253 303  19 |
9  |   10 |   36  484   38 | -125  12   2 |  -774 472  97 |
10 |   11 |   58  422   13 | -142  76   3 |   304 347  24 |
11 |   12 |   15  119   38 |  113   4   1 |  -605 115  24 |
12 |   13 |   15  702   52 | -285  19   3 | -1727 683 193 |
13 |   15 |   36  207   20 | -160  38   3 |   335 169  18 |
14 |   16 |   51  158   22 | -192  69   5 |  -218  89  11 |
15 |   17 |   15   35   56 | -304  20   4 |   269  15   5 |
16 |   18 |   22  954  110 | 2416 951 346 |   135   3   2 |
17 |   19 |   51  140   16 |  -88  21   1 |  -211 119  10 |
18 |   20 |   51  449   21 |  -23   1   0 |  -477 448  52 |
19 |   21 |   44  371   22 | -304 153  11 |   363 218  26 |
20 |   22 |   15  184   44 |  147   6   1 |   807 178  42 |
21 |   23 |   22  186   32 |  369  78   8 |  -436 108  18 |
22 |   24 |   29  185   31 | -431 144  15 |   230  41   7 |
23 |   25 |   15  934  149 | 3380 919 451 |   429  15  12 |
24 |   26 |   22  620   41 |  104   5   1 | -1184 615 136 |
25 |   27 |   22   32   46 |   11   0   0 |  -287  32   8 |
26 |   29 |   44  371   22 | -304 153  11 |   363 218  26 |
27 |   30 |   66  195   15 | -229 195   9 |     4   0   0 |

Columns:
         name   mass  qlt  inr    k=1 cor ctr     k=2 cor ctr  
1  | カラオケ |  109  378   36 | -253 159  19 |   297 219  43 |
2  |     バー |   73  324   58 | -385 152  29 |   409 172  54 |
3  | パチンコ |   66   23   99 | -118   8   2 |   166  15   8 |
4  | レストラ |  161  446   45 | -160  74  11 |  -358 372  92 |
5  |       駅 |   66  743  121 | -186  15   6 | -1281 728 478 |
6  |       家 |   95  108   98 |  298  71  23 |  -216  37  20 |
7  |     学校 |   22  956  244 | 3576 941 757 |   454  15  20 |
8  |   居酒屋 |  131  191   40 | -252 169  23 |    90  22   5 |
9  |     車内 |   95  510  120 | -358  83  33 |   813 427 278 |
10 |     職場 |  124  380   78 |  536 377  96 |   -47   3   1 |
11 |     食堂 |   58    9   61 |  -80   5   1 |   -71   4   1 |
Warning message:
In abbreviate(cnames.temp, 4) : abbreviate used with non-ASCII chars
> 

plot()で図が現れる。

2つの固有値をX軸とY軸にした散布図。

mass=T(TRUE)とすると、

plot(ca.NO.basho, mass=T)

同じ回答の複数人は、 大きな〇で描かれる。

行(回答者)だけプロットするには、 以下の通りにインプット。

plot(ca.NO.basho, mass=T, what=c("active","none"))

列(受動喫煙被害場所)をプロットするには、 以下の通り。

plot(ca.NO.basho, mass=T, what=c("none","active"))

参考図書