CA_step2_CAandMCA.Rmd
/Users/kazuo/Dropbox/RStudio/計量分析セミナー2023/CA_MCA_Q16.Rmd
kazuo.fujimoto2007@gmail.com
Mon Sep 4 13:21:56 JST 2023
1 Package を読み込む
library(tidyverse)
library(gtsummary)
library(GDAtools)
library(vcd)
library(FactoMineR)
library(factoextra)
library(ggrepel)
library(showtext)
showtext_auto(TRUE)2 データ表を指示行列（indicator行列）に変換してCAを行う：MACの一つの方法
2.1 データを読み込む
- Le Roux＆Rouanetの『多重対応分析』のサンプルデータとして提供されているものを、日本語化したもの。
- もとデータはここからダウンロード可能。
- Le Roux先生のHome https://helios2.mi.parisdescartes.fr/~lerb/index_en.html
- ここの、「Data and Script」 https://helios2.mi.parisdescartes.fr/~lerb/Logiciels/software.html にTaste Exampleがある。https://helios2.mi.parisdescartes.fr/~lerb/Logiciels/Data/Taste_Example.xls
load("Data/taste_d_J.rda")
.d_J## # A tibble: 1,253 × 9
## ID Isup TV Film Art Eat Gender Age Income
## <int> <fct> <fct> <fct> <fct> <fct> <fct> <fct> <fct>
## 1 1 Active TV-メロドラマ 映画-アクション 芸術… 外食… 女性 55-64 £20-29
## 2 2 Active TV-メロドラマ 映画-ホラー 芸術… 外食… 女性 45-54 <£9
## 3 3 Active TV-自然 映画-アクション 芸術… 外食… 女性 55-64 <£9
## 4 4 Active TV-メロドラマ 映画-時代劇 芸術… 外食… 女性 65+ £10-19
## 5 5 Active TV-コメディー 映画-ホラー 芸術… 外食… 女性 35-44 £10-19
## 6 6 Active TV-コメディー 映画-ホラー 芸術… 外食… 女性 18-24 <£9
## 7 7 Active TV-ニュース 映画-アクション 芸術… 外食… 女性 25-34 £10-19
## 8 8 Active TV-ニュース 映画-ドキュメンタ… 芸術… 外食… 男性 65+ £10-19
## 9 9 Active TV-メロドラマ 映画-時代劇 芸術… 外食… 女性 65+ <£9
## 10 10 Active TV-ニュース 映画-アクション 芸術… 外食… 女性 65+ £10-19
## # ℹ 1,243 more rows2.2 指示行列化
dichotom(.d_J[1:1215,3:6]) -> .d_J.ind2.3 指示行列をCAする（FactoMineR::CAを使う）
res.CA.ind <- CA(.d_J.ind,graph=FALSE)2.4 resultを図示
fviz_ca_col(res.CA.ind) + coord_fixed(ratio = 1)
fviz_ca_row(res.CA.ind,alpha.row = "contrib") + coord_fixed(ratio = 1)
.d_J %>% summary(maxsum = 20)## ID Isup TV Film
## Min. : 1 Active:1215 TV-コメディー:164 映画-アクション :405
## 1st Qu.: 314 supp : 38 TV-ドラマ :135 映画-コメディ :247
## Median : 627 TV-映画 :122 映画-時代劇 :141
## Mean : 627 TV-自然 :162 映画-ドキュメンタリー:100
## 3rd Qu.: 940 TV-ニュース :228 映画-ホラー : 64
## Max. :1253 TV-警察もの : 85 映画-ミュージカル : 87
## TV-メロドラマ:216 映画-ロマンス :104
## TV-スポーツ :141 映画-SF :105
## Art Eat Gender
## 芸術-印象派 :126 外食先-フィッシュ＆チップス:113 男性:534
## 芸術-風景画 :645 外食先-フランス料理店 : 99 女性:719
## 芸術-現代美術 :117 外食先-インド料理店 :419
## 芸術-パフォーマンス・アート:110 外食先-イタリア料理店 :240
## 芸術-肖像画 :123 外食先-パブ :284
## 芸術-ルネッサンス美術 : 59 外食先-ステーキハウス : 98
## 芸術-静物画 : 73
##
## Age Income
## 18-24:102 ? :166
## 25-34:259 <£9 :232
## 35-44:269 >=£60 :130
## 45-54:196 £10-19:261
## 55-64:184 £20-29:206
## 65+ :243 £30-39:127
## £40-59:131
## 2.5 GDAtools::speMCAを使ってMCAを行った
.d_J[1:1215,3:6] %>% speMCA %>% flip.mca(dim=1) -> res.speMCA
ggcloud_variables(res.speMCA) + coord_fixed(ratio = 1)
ggcloud_indiv(res.speMCA,type = 'inames') + coord_fixed(ratio = 1)
3 MCAの基本的な使い方
3.1 留置A票の調査票
3.2 データ：SSM2005の留置A票の問16 男女性別役割意識のデータを使う。
- このデータは、SSJDAのオンライン集計で、SSM2005、留置A票の問16アイウと面接調査票の性別、年齢を、多元クロス集計し、それをRに取り込んで個票を復元したもの。
- 手順は、https://rpubs.com/kfj419/1077498 に公開。こうして作られたRオブジェクト、Q16ABC.rdaを読み込みます。
- なお、このようにして生成したデータの利用については問題なし、という確認を事務局にとってあります。
fname <- load("Data/Q16ABC.rda")
Q16ABC## # A tibble: 2,827 × 6
## Q16a Q16b Q16c Sex Age Age10
## <fct> <fct> <fct> <fct> <fct> <fct>
## 1 A A A 男性 20 20
## 2 A A A 男性 21 20
## 3 A A A 男性 22 20
## 4 A A A 男性 28 20
## 5 A A A 男性 29 20
## 6 A A A 男性 29 20
## 7 A A A 男性 30 30
## 8 A A A 男性 34 30
## 9 A A A 男性 34 30
## 10 A A A 男性 37 30
## # ℹ 2,817 more rows3.3 データの概要を確認
q16a <- "男性は外で働き、女性は家庭を守るべき"
q16b <- "男の子と女の子は違った育てかたをすべき"
q16c <- "家事や育児は、男性より女性が向いている"Q16ABC %>% summary## Q16a Q16b Q16c Sex Age Age10
## A : 246 A :277 A : 572 女性:1484 57 : 98 20:325
## B : 727 B :754 B :1141 男性:1343 65 : 96 30:516
## C : 677 C :647 C : 418 55 : 89 40:562
## D :1057 D :985 D : 545 61 : 87 50:684
## DKNA: 120 DKNA:164 DKNA: 151 64 : 87 60:726
## 58 : 81 70: 14
## (Other):22894 データの基本的分析
- 簡単に一覧表を出力するfunctionがあります。
- gtsummary::tbl_summary など
- explore::report など
使い方は、ネット上にある説明をみてください。
gtsummary::tbl_summary(Q16ABC[,-5])| Characteristic | N = 2,8271 |
|---|---|
| Q16a | |
| A | 246 (8.7%) |
| B | 727 (26%) |
| C | 677 (24%) |
| D | 1,057 (37%) |
| DKNA | 120 (4.2%) |
| Q16b | |
| A | 277 (9.8%) |
| B | 754 (27%) |
| C | 647 (23%) |
| D | 985 (35%) |
| DKNA | 164 (5.8%) |
| Q16c | |
| A | 572 (20%) |
| B | 1,141 (40%) |
| C | 418 (15%) |
| D | 545 (19%) |
| DKNA | 151 (5.3%) |
| Sex | |
| 女性 | 1,484 (52%) |
| 男性 | 1,343 (48%) |
| Age10 | |
| 20 | 325 (11%) |
| 30 | 516 (18%) |
| 40 | 562 (20%) |
| 50 | 684 (24%) |
| 60 | 726 (26%) |
| 70 | 14 (0.5%) |
| 1 n (%) | |
library(explore)
explore::report(Q16ABC,output_file = "Q16ABC.report",output_dir = ".")##
##
## processing file: template_report_variable.Rmd##
|
| | 0%
|
|...... | 11%
|
|........... | 22% (unnamed-chunk-1-2)
|
|................. | 33%
|
|...................... | 44% (unnamed-chunk-3-4)
|
|............................ | 56%
|
|................................. | 67% (unnamed-chunk-5-6)
|
|....................................... | 78%
|
|............................................ | 89% (unnamed-chunk-7-8)
|
|..................................................| 100% ## output file: /Users/kazuo/Dropbox/RStudio/計量分析セミナー2023/template_report_variable.knit.md## /Applications/RStudio.app/Contents/Resources/app/quarto/bin/tools/pandoc +RTS -K512m -RTS /Users/kazuo/Dropbox/RStudio/計量分析セミナー2023/template_report_variable.knit.md --to html4 --from markdown+autolink_bare_uris+tex_math_single_backslash --output /Users/kazuo/Dropbox/RStudio/計量分析セミナー2023/Q16ABC.report.html --lua-filter /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/library/rmarkdown/rmarkdown/lua/pagebreak.lua --lua-filter /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/library/rmarkdown/rmarkdown/lua/latex-div.lua --embed-resources --standalone --variable bs3=TRUE --section-divs --template /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/library/rmarkdown/rmd/h/default.html --no-highlight --variable highlightjs=1 --variable theme=bootstrap --mathjax --variable 'mathjax-url=https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML' --include-in-header /var/folders/_3/hq2q61m11yv4q4748clqs0rc0000gn/T//RtmpYPBISK/rmarkdown-str98b2199f1579.html##
## Output created: Q16ABC.report.html5 MCAを行う前の確認 mosaicで行います。
- Q16A,B,Cと年代and/or性別 の関係をみたい。
- Q16xと年代、Q16xと性別の関係を確認
- クロス集計、とmosaic 表示
5.1 クロス集計
- 保守的 A-B-C-D リベラル
Q16ABC %>% xtabs(~ Age10 + Q16a,.)## Q16a
## Age10 A B C D DKNA
## 20 20 63 77 151 14
## 30 28 131 126 204 27
## 40 27 147 158 212 18
## 50 50 166 165 279 24
## 60 117 214 151 207 37
## 70 4 6 0 4 06 mosaciで表示します。
6.1 各設問と年代のクロス集計。
- 青が残差が期待値よりプラス大きいセル、赤は残差が期待値よりマイナスに大きいセル
Q16ABC %>% xtabs(~ Age10 + Q16a,.) %>% mosaic(shade=TRUE,main=str_c("Q16a",q16a))
Q16ABC %>% xtabs(~ Age10 + Q16b,.) %>% mosaic(shade=TRUE,main=str_c("Q16b",q16b))
Q16ABC %>% xtabs(~ Age10 + Q16c,.) %>% mosaic(shade=TRUE,main=str_c("Q16c",q16c))
6.2 集計に性別を加えてみます
Q16ABC %>% xtabs(~ Sex + Age10+ Q16a,.) %>% mosaic(shade=TRUE,main=str_c("Q16a",q16a))
Q16ABC %>% xtabs(~ Sex + Age10 + Q16b,.) %>% mosaic(shade=TRUE,main=str_c("Q16b",q16b))
Q16ABC %>% xtabs(~ Sex + Age10 + Q16c,.) %>% mosaic(shade=TRUE,main=str_c("Q16c",q16c))
7 MCA を行う 使うのはGDAtools
- Q16A,B,Cで空間を生成し、あとから、そこに、年代、性別をplotします。
res.speMCA1 <- speMCA(Q16ABC[,1:3])7.1 変数空間と個体空間をplot
ggcloud_variables(res.speMCA1) # 変数空間（列空間） 
ggcloud_indiv(res.speMCA1) # 個体空間（行空間）
この変数空間でわかること。
- ３問とも、DKNAが右側にあつまっていて、それが、軸１をひっぱっている。
- そのため、Q16ABCの関係が左側に密集していてよみとれない。
7.2 specific MCAをつかい、junkカテゴリ（ここではDKNA）を除去する処理を行う。
- ただし、この除去は、最初のデータ表の周辺度数に影響がないように、残差行列を計算することで実現しないと、その先の処理にとってよろしくない。『多重対応分析』pp81、津田塾大学紀要55号。
- GDAtools::speMCAでは、このjunkカテゴリをカテゴリ番号で指定する。この仕様はFactoMineR::MCAでの処理と同じ。オプションexcl=で指定。そのカテゴリ番号は、GDAtools:: getindexcat で取得できる、
7.3 DKNAというラベルを有するカテゴリ番号を取得する。
getindexcat(Q16ABC[,1:3]) %>% str_detect("DKNA") %>% which(TRUE) -> excl_list7.4 これを指定したspeMCAのresultを、res.speMCA2に格納
res.speMCA2 <- speMCA(Q16ABC[,1:3],excl = excl_list)7.5 生成された空間を図示する
ggcloud_variables(res.speMCA2) + ggtitle("Q16abcの変数空間") -> pv12
ggcloud_indiv(res.speMCA2) + ggtitle("Q16abcの個体空間") -> pi12
pv12
pi12
8 MCAのresultの考察
8.1 固有値eigの確認
- mがついているのは修正分散率。『多重対応分析』p56「分散率および修正分散率」この章（３章）の訳註16には、Grenacreの修正分散率についてふれている（2007なので、CAiP2。『対応分析の理論と実践』2017=2020では、p251 Appendix A で説明されている。）。
- 修正前のeigは、12個生成されている。修正後は、3個。
res.speMCA2$eig## $eigen
## [1] 0.606317121 0.566282681 0.436087287 0.297435869 0.267769168 0.235737558
## [7] 0.211372330 0.195324230 0.183771451 0.036623948 0.007879875 0.006689602
##
## $rate
## [1] 19.8708382 18.5587890 14.2918938 9.7478693 8.7756021 7.7258298
## [7] 6.9273079 6.4013633 6.0227439 1.2002771 0.2582473 0.2192384
##
## $cum.rate
## [1] 19.87084 38.42963 52.72152 62.46939 71.24499 78.97082 85.89813
## [8] 92.29949 98.32224 99.52251 99.78076 100.00000
##
## $mrate
## [1] 53.48 38.94 7.58
##
## $cum.mrate
## [1] 53.48 92.42 100.008.2 グラフ表示
res.speMCA2$eig[1:3] %>% as_tibble() %>% ggplot(aes(x=1:12,y=rate)) + geom_col(fill="lightblue")
8.3 修正前の寄与率によるスクリープロット
res.speMCA2$eig[1:3] %>% as_tibble() %>% ggplot(aes(x=1:12,y=rate)) + geom_col(fill="lightblue") +
geom_point(aes(y=cum.rate)) + geom_line(aes(y=cum.rate)) + geom_text_repel(aes(y=cum.rate,label=round(cum.rate,2)))
8.4 修正寄与率の図示
res.speMCA2$eig[4:5] %>% as_tibble() %>% ggplot(aes(x=1:3,y=mrate)) + geom_col(fill="lightblue") +
geom_point(aes(y=cum.mrate)) + geom_line(aes(y=cum.mrate)) + geom_text_repel(aes(y=cum.mrate,label=cum.mrate))
9 変数空間の考察
9.1 変数カテゴリの軸への寄与を確認する
ggaxis_variables(res.speMCA2,axis = 1,prop = "ctr")
ggaxis_variables(res.speMCA2,axis = 2,prop = "ctr")
9.2 軸ごとに方向含めて確認する
tabcontrib(res.speMCA2,dim=1,best = FALSE) # best=TRUE だと、平均値より寄与しているものだけをリストアップ## Variable Category Weight Quality of representation Contribution (left)
## 4 Q16a D 1057 0.598
## 2 B 727 0.210 8.57
## 3 C 677 0.092 3.85
## 1 A 246 0.022 1.11
## 12 Q16c D 545 0.585
## 10 B 1141 0.143 4.69
## 9 A 572 0.043 1.89
## 11 C 418 0.011 0.52
## 8 Q16b D 985 0.578
## 6 B 754 0.159 6.41
## 7 C 647 0.104 4.4
## 5 A 277 0.027 1.34
## Contribution (right) Total contribution Cumulated contribution
## 4 20.58 34.1 34.1
## 2
## 3
## 1
## 12 25.95 33.05 67.15
## 10
## 9
## 11
## 8 20.71 32.85 100
## 6
## 7
## 5
## Contribution of deviation Proportion to variable
## 4 33.58 98.47
## 2
## 3
## 1
## 12 32.79 99.23
## 10
## 9
## 11
## 8 32.74 99.67
## 6
## 7
## 5第一軸への寄与は、大きものは、左側に、Q16a-B、右側へのQ16a-D、Q16c-D、Q16b-Dであることがわかる。つまり左側には、保守的、どちらかといえば、が位置し、右側には、リベラルが位置している。
これが第１軸への各カテゴリの寄与。
9.3 分析の順番
- まず、変数空間の分析から軸に名前をつける。
- 次に、個体空間の分析に入る-> 個体空間に、追加変数で、性別、年代を投影する。（詳細は、次の資料GDAで説明します）
ggadd_supvars(pi12,resmca = res.speMCA2,vars = Q16ABC[,c(4,6)],textsize = 4,shapes = TRUE,shapesize = 2)
9.3.1 年代-性別の合成変数をplotします。
ggadd_interaction(pi12,res.speMCA2,v1 = Q16ABC$Sex,v2 = Q16ABC$Age10) + coord_fixed(ratio = 1)## Warning: ggrepel: 2 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
9.4 集中楕円で、追加変数に対する個体の分布を表示します。
ggadd_ellipses(pi12,resmca = res.speMCA2,var = Q16ABC$Sex) + coord_fixed(ratio = 1)
ggadd_ellipses(pi12,resmca = res.speMCA2,var = Q16ABC$Age10) + coord_fixed(ratio = 1)