1 概要

このスクリプトは、(Le Roux et al. 2021)にあるspeMCA、（および(Greenacre and 藤本一男 2020)にあるサブセットMCA）の数理的な処理を、MCA自体を、指標行列化した対象データに対するCAとして実行することまでさかのぼり、juncカテゴリを１）残差行列Sをつくるまでは行和に変更を加えず、２）S'としてjuncカテゴリ列を除去した残差行列S'へのSVDによって空間生成を行うことで、speMCAの処理を確かめている。SVDの結果から、標準座標を計算する際には、列和ベクトルに、juncカテゴリを除去したc’をもちいる必要がある。

2 環境設定

library(ca)
library(FactoMineR)
library(factoextra)
##  要求されたパッケージ ggplot2 をロード中です
## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
library(tidyverse)
## ── Attaching packages
## ───────────────────────────────────────
## tidyverse 1.3.2 ──
## ✔ tibble  3.1.8      ✔ dplyr   1.0.10
## ✔ tidyr   1.2.1      ✔ stringr 1.4.1 
## ✔ readr   2.1.3      ✔ forcats 0.5.2 
## ✔ purrr   0.3.5      
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
library(GDAtools)
library(showtext)
##  要求されたパッケージ sysfonts をロード中です 
##  要求されたパッケージ showtextdb をロード中です
showtext_auto(TRUE)

3 データの読み込み確認

3.1 GDAtoolsに付属しているMusicを使う

data("Music")

3.2 内容確認

Music %>% head()
##      FrenchPop Rap Rock Jazz Classical Gender   Age OnlyMus Daily
## 2124        NA  No  Yes   No        No    Men 25-49    Rare    No
## 4485        No Yes   No   NA        No  Women 15-24   Often   Yes
## 4731        No  No  Yes   No        No  Women 15-24   Daily   Yes
## 1411        No  No  Yes   No        No  Women 15-24    Rare   Yes
## 3986       Yes  No   No   No       Yes  Women 25-49   Never    No
## 3110       Yes  No  Yes   No        No    Men   50+   Never    No

4 speMCAを実行

Music[,1:5] %>% speMCA -> res.speMCA

5 rpm行列をindicator行列に変換しCAする

#ca::caconv(Music,from = "rpm",to="ind") -> Music.ind  # ca::caconv はバグがある
GDAtools::dichotom(Music) -> Music.ind
Music.ind %>% summarize_all(sum)
##   FrenchPop.No FrenchPop.Yes FrenchPop.NA Rap.No Rap.Yes Rap.NA Rock.No
## 1          194           301            5    417      77      6     360
##   Rock.Yes Rock.NA Jazz.No Jazz.Yes Jazz.NA Classical.No Classical.Yes
## 1      135       5     395       95      10          351           142
##   Classical.NA Gender.Men Gender.Women Age.15-24 Age.25-49 Age.50+
## 1            7        253          247        78       204     218
##   OnlyMus.Daily OnlyMus.Often OnlyMus.Rare OnlyMus.Never Daily.No Daily.Yes
## 1            58           121           64           257      303       197

5.1 FactoMineR::CAでMusic.indをCAする

res.CA <- CA(Music.ind)

5.2 変数グラフ（CAの列変数）

fviz_ca_col(res.CA) + coord_fixed(ratio = 1)

6 行列計算のstep by step で speMCAをやる

6.1 junc指定するカテゴリ列を確認する

getindexcat(Music[,1:5])
##  [1] "FrenchPop.No"  "FrenchPop.Yes" "FrenchPop.NA"  "Rap.No"       
##  [5] "Rap.Yes"       "Rap.NA"        "Rock.No"       "Rock.Yes"     
##  [9] "Rock.NA"       "Jazz.No"       "Jazz.Yes"      "Jazz.NA"      
## [13] "Classical.No"  "Classical.Yes" "Classical.NA"

6.2 そのカテゴリ列番号を取得する

getindexcat(Music[,1:5]) %>% str_which(".NA") -> excl_list
excl_list
## [1]  3  6  9 12 15

6.2.1 まず、indicator行列化したMusicデータを行列処理でCAを行う

\(S=D_r^{-1/2}(P-rc^t)D_c^{-1/2}\) (1)

\(S=UD_{\alpha}V^t 　ここで　U^tU＝V^tV=I\) (2)

これをもとに、標準座標（\(\Phi\), \(\Gamma\)）、主座標(\(F\),\(G\))を求める。(3)

\(\Phi= D_r^{-1/2}U\) (3.1)

\(\Gamma= D_c^{-1/2}V\) (3.2)

\(F=D_r^{-1/2}U D_\alpha = \Phi D_\alpha\) (3.3)

\(G=D_c^{-1/2}V D_\alpha = \Gamma D_\alpha\) (3.4)

6.2.2 junc カテゴリ番号をexclにベクトルとして取得しているとする。

残差行列\(S\)の計算までは同じ。同じ行和\(r\)、同じ列和\(c\)のままである。

この\(S\)から、exclで示されるカテゴリ列をはずして\(S’\)を作成する。

この\(S'\)に対して、SVDを行う。得られた三つの行列を、\(S'=U'D'_{\alpha}V'^t\) とする。

この三つ行列をもちいて、(3)と同様に、標準座標、主座標を求めればよい。ただし、

\(\Phi'= D_r^{-1/2}U\) の行和\(r\)はそのままでいいが、列和\(c\)は、juckカテゴリを除去した\(P\)に対する列和\(c'\)でなければ、大きさが合わない。

そのため、 \(\Gamma '= D_{c'}^{-1/2}V\)　 (3.2’)

\(G’=D_{c'}^{-1/2}V D_\alpha = \Gamma' D_\alpha\) (3.4’)

となる。

7 指標行列化してCAを行う

GDAtools::dichotom(Music[,1:5]) -> Music.ind
tab <- Music.ind
tab.P <- as.matrix(tab)/sum(tab)
tab.r <- apply(tab.P,1,sum) # 行和
tab.c <- apply(tab.P,2,sum) # 列和
tab.S <- diag(1/sqrt(tab.r))%*%(tab.P -tab.r %*% t(tab.c)) %*% diag(1/sqrt(tab.c))
tab.svd <- svd(tab.S)

7.1 普通のMCA

7.2 行と列の標準座標

tab.rsc <- diag(1/sqrt(tab.r)) %*% tab.svd$u
tab.csc <- diag(1/sqrt(tab.c)) %*% tab.svd$v

7.3 行と列の主座標

tab.rpc <- tab.rsc %*% diag(tab.svd$d)
tab.cpc <- tab.csc %*% diag(tab.svd$d)

7.4 図示する

showtext_auto(FALSE)
tab.cpc[,1:2]
##              [,1]        [,2]
##  [1,]  0.15517193 -0.07210583
##  [2,] -0.08101899 -0.07115303
##  [3,] -1.14332804  7.08111858
##  [4,] -0.11157641 -0.20256396
##  [5,]  0.65496235  0.62665822
##  [6,] -0.65078957  6.03608175
##  [7,]  0.25849045 -0.13642956
##  [8,] -0.65708449  0.42172048
##  [9,] -0.87003074 -1.56352432
## [10,]  0.38323247  0.01120467
## [11,] -1.59452604  0.05231874
## [12,]  0.01031479 -0.93961242
## [13,]  0.42791761  0.14791017
## [14,] -1.13491140 -0.34461398
## [15,]  1.56547692 -0.42589792
plot(tab.cpc[,1],tab.cpc[,2])

showtext_auto(TRUE)

8 speMCAを行う

GDAtools::dichotom(Music[,1:5]) -> Music.ind
tab <- Music.ind
tab.P <- as.matrix(tab)/sum(tab)
tab.r <- apply(tab.P,1,sum) # 行和
tab.c <- apply(tab.P,2,sum) # 列和
tab.c2 <- tab.c[-excl_list]
tab.S <- diag(1/sqrt(tab.r))%*%(tab.P -tab.r %*% t(tab.c)) %*% diag(1/sqrt(tab.c))
tab.S2 <- tab.S[,-excl_list] # 【speMCAなので、juncカテゴリ例をSから外す】
tab.svd <- svd(tab.S2)

8.1 行と列の標準座標（standard coordinate 平均ゼロ、分散1）

tab.rsc <- diag(1/sqrt(tab.r)) %*% tab.svd$u # 行和は、変更なし
tab.csc <- diag(1/sqrt(tab.c2)) %*% tab.svd$v　# 列和の総数は、juncカテゴリ列を除去して再計算したもの：C2

8.2 行と列の主座標：標準座標が計算されていれば、主座標との関係は特異値でスケールするだけ

tab.rpc <- tab.rsc %*% diag(tab.svd$d)
tab.cpc <- tab.csc %*% diag(tab.svd$d)

8.3 変数ポイントの値を表示

tab.cpc[,1:2]
##              [,1]        [,2]
##  [1,]  0.13619669 -0.03863579
##  [2,] -0.08386882  0.03246438
##  [3,] -0.11780763  0.31425633
##  [4,]  0.64211212 -1.69801595
##  [5,]  0.24788917  0.41151917
##  [6,] -0.65151904 -1.11088162
##  [7,]  0.38738927  0.08521378
##  [8,] -1.60952273 -0.36109891
##  [9,]  0.46286382 -0.17772100
## [10,] -1.16596062  0.43836406

8.4 plotしてポイントの位置を確かめる

showtext_auto(FALSE)
plot(tab.cpc[,1],tab.cpc[,2],asp=1,main="行列演算でもとめたMusicのspeMCA、変数空間")

showtext_auto(TRUE)

9 speMCAの検算：GDAtoolsのsepMCAで確認

res.speMCA <- speMCA(Music[,1:5],excl = excl_list)
ggcloud_variables(res.speMCA) + ggtitle("GDAtools::speMCAで求めたMusicの変数空間")

9.1 座標値を比べる

res.speMCA$var$coord[,1:2]
##                     dim.1       dim.2
## FrenchPop.No   0.13619669 -0.03863579
## FrenchPop.Yes -0.08386882  0.03246438
## Rap.No        -0.11780763  0.31425633
## Rap.Yes        0.64211212 -1.69801595
## Rock.No        0.24788917  0.41151917
## Rock.Yes      -0.65151904 -1.11088162
## Jazz.No        0.38738927  0.08521378
## Jazz.Yes      -1.60952273 -0.36109891
## Classical.No   0.46286382 -0.17772100
## Classical.Yes -1.16596062  0.43836406

9.2 行列「手計算」によるspeMCAのリザルト

tab.cpc[,1:2]
##              [,1]        [,2]
##  [1,]  0.13619669 -0.03863579
##  [2,] -0.08386882  0.03246438
##  [3,] -0.11780763  0.31425633
##  [4,]  0.64211212 -1.69801595
##  [5,]  0.24788917  0.41151917
##  [6,] -0.65151904 -1.11088162
##  [7,]  0.38738927  0.08521378
##  [8,] -1.60952273 -0.36109891
##  [9,]  0.46286382 -0.17772100
## [10,] -1.16596062  0.43836406

参考文献

Greenacre, Michael, and 藤本一男. 2020. 対応分析の理論と実践: 基礎・応用・展開. 東京: オーム社.
Le Roux, Brigitte, Henry Rouanet, 大隅昇, 小野裕亮, and 鳰真紀子. 2021. 多重対応分析. 東京: オーム社.