前回の記事の続きです。 bob3.hatenablog.com
- 結論
- 使える同時布置図の描き方
- コレスポンデンス分析の基本的な流れ
- 標準化残差の算出
- 標準化残差(Z)を特異値分解する
- ちょっと脱線 カイ二乗検定と残差分析
- 座標の重みづけ
- 座標の組み合わせ
- 指標化残差が角度と長さで表現できてるか確認
- 大きな集計表で再確認
- 結論（再掲）
- 参考リンク
結論
今回も最初に結論を。
- 縦横のスケールを合わせるのが大前提です。
- そうしないと見かけ上の角度が歪んでしまいます。
- コレスポンデンス分析の同時布置図を描くときは、対称バイプロットがおすすめです。
- 指標化残差を正確に角度として表現できて、なおかつ見やすいので。
- 従来の同時布置図はフレンチプロットが多いと思いますが、正確でないのであえて選ぶ理由はないと思います。
- ただ、「だいたい合ってる」ので、角度で読む分には大きく間違えることはなさそう。
なお、tarotanさんの記事も大変勉強になるので参照してください。 tarotan.hatenablog.com
使える同時布置図の描き方
前回の記事では指標化残差を角度で表現したのがコレポンの同時布置図だよ、という説明をしたのですが、肝心の描画に使う座標の作り方は関数任せでした。
今回はいくつかの描画手法について検討し、どれを使えばいいのかを考えてみたいと思います。
まず今回のサンプルデータですが、『対応分析入門』に載っているノルウェイの犯罪統計のデータを引用します。
ミニマムなデータで計算過程を追いたい意図です。
(dat <- c( 395, 2456, 1758, 147, 153, 916, 694, 327, 1347 ) |> matrix( nrow = 3, byrow = TRUE, dimnames = list( 地域 = c("オスロ", "中部", "北部"), 犯罪 = c("強盗", "詐欺", "破壊") ) )) ## 犯罪 ## 地域 強盗 詐欺 破壊 ## オスロ 395 2456 1758 ## 中部 147 153 916 ## 北部 694 327 1347
コレスポンデンス分析の基本的な流れ
以下にコレスポンデンス分析の基本的な計算過程を示します。
- クロス集計表の観測割合(P)を算出する。
- 各セルの値を総度数で割った値。
- 行と列の質量(masses)を算出する。
- 行、列の各カテゴリ毎の観測割合の合計。
- 各セルの期待割合(E)を算出する。
- 行と列の質量をセルごとにかけ合わせた値。
- 行と列に関連が無かった場合に想定される各セルの割合。
- 割合の残差(R)を算出する。
- 各セルの観測割合(P)から期待割合(E)を引いた値。
- 指標化残差(I)を算出する。
- 各セルの割合の残差を期待割合で割った値。
- 標準化残差（Z、ピアソン残差）を算出
- 指標化残差に期待値の平方根を乗じた値。
- （観測値-期待値）/sqrt(期待値)
- 『コレスポンデンス分析の利用法』でいうところの標準化、行と列の純粋な関連性。
- 特異値分解(SVD)で次元縮約。
- 標準化残差(Z)に対してSVDで次元縮約し、2次元のスコアを得る
- 固有値、寄与率もここで算出
- スコアに重みづけを行い同時布置図のための座標を得る。
- 座標の重みづけ
- 標準座標
- 主座標
- 対称バイプロット
- 座標の組み合わせ
- 対称プロット（主座標×主座標）
- 非対称プロット（主座標×標準座標）
- 対称バイプロット
- 座標の重みづけ
1から5までは前回の記事で扱ったのでRのコードを示すにとどめ、今回は6の指標化残差から標準化残差を出すところから始めます
n <- sum(dat) # 総度数(n) P <- dat/n # 観測割合(P) row.masses <- rowSums(P) # 行質量（row masses、平均列プロファイル） column.masses <- colSums(P) # 列質量（column masses） E <- row.masses %o% column.masses # 期待値 R <- P - E # 残差 I <- R / E # 指標化残差 list( `総度数(n)` = n, `観測割合(P)` = P, `行質量(row.masses)` = row.masses, `列質量(column.masses)` = column.masses, `期待値(E)` = E, `残差(R)` = R, `指標化残差(I)` = I ) |> print() ## $`総度数(n)` ## [1] 8193 ## ## $`観測割合(P)` ## 犯罪 ## 地域 強盗 詐欺 破壊 ## オスロ 0.04821189 0.29976809 0.2145734 ## 中部 0.01794215 0.01867448 0.1118028 ## 北部 0.08470646 0.03991212 0.1644086 ## ## $`行質量(row.masses)` ## オスロ 中部 北部 ## 0.5625534 0.1484194 0.2890272 ## ## $`列質量(column.masses)` ## 強盗 詐欺 破壊 ## 0.1508605 0.3583547 0.4907848 ## ## $`期待値(E)` ## 強盗 詐欺 破壊 ## オスロ 0.08486708 0.20159365 0.27609267 ## 中部 0.02239062 0.05318678 0.07284198 ## 北部 0.04360279 0.10357426 0.14185017 ## ## $`残差(R)` ## 犯罪 ## 地域 強盗 詐欺 破壊 ## オスロ -0.036655194 0.09817444 -0.06151925 ## 中部 -0.004448475 -0.03451230 0.03896078 ## 北部 0.041103669 -0.06366214 0.02255847 ## ## $`指標化残差(I)` ## 犯罪 ## 地域 強盗 詐欺 破壊 ## オスロ -0.4319130 0.4869917 -0.2228210 ## 中部 -0.1986758 -0.6488887 0.5348671 ## 北部 0.9426844 -0.6146521 0.1590303
標準化残差の算出
割合の標準化残差（Z、ピアソン残差）は指標化残差(I)に期待値(E)の平方根を乗じることで得られます。
(Z <- I * sqrt(E)) ## 犯罪 ## 地域 強盗 詐欺 破壊 ## オスロ -0.12582469 0.2186553 -0.11708024 ## 中部 -0.02972885 -0.1496484 0.14435664 ## 北部 0.19684458 -0.1978132 0.05989557
標準化残差(Z)を特異値分解する
標準化残差(Z)を特異値分解(singular value decomposition、SVD)に掛け、特異値と特異ベクトルを取り出します。
この特異ベクトルが同時布置図の座標のもとになります。
特異値分解とは何か？の説明は私の手にあまりますので各自で勉強してみてください（すいません……
res.SVD <- svd(Z) # 特異値分解 res.SVD$d # 特異値 ## [1] 4.212385e-01 1.596575e-01 7.558903e-18 U <- res.SVD$u[,1:2] # 左特異行列。行に対応。2次元だけ取り出す。 rownames(U) <- rownames(P) U ## [,1] [,2] ## オスロ -0.6600450 0.04227544 ## 中部 0.3840335 -0.83910599 ## 北部 0.6456461 0.54232272 V <- res.SVD$v[,1:2] # 右特異ベクトル。列に対応。2次元だけ取り出す。 rownames(V) <- colnames(P) V ## [,1] [,2] ## 強盗 0.4717636 0.7915672 ## 詐欺 -0.7822401 0.1724693 ## 破壊 0.4068654 -0.5862387
特異値を二乗したものが固有値になります。割合を出せば寄与率になります。
eigenvalues <- res.SVD$d^2 # 固有値 round(eigenvalues, 4) ## [1] 0.1774 0.0255 0.0000 eigen_prop <- prop.table(eigenvalues) # 寄与率 round(eigen_prop, 4) ## [1] 0.8744 0.1256 0.0000 cumsum(eigen_prop) # 累積寄与率 ## [1] 0.8743891 1.0000000 1.0000000
ちょっと脱線 カイ二乗検定と残差分析
コレスポンデンス分析と良く併用されるカイ二乗検定も、コレスポンデンス分析の計算途中の数値から計算できます。
カイ二乗値は残差の二乗を期待値で割って総度数を掛ければ出ます。
カイ二乗値が出ればp値も計算できます。
(chisq <- sum(R^2 / E * n)) # カイ二乗値 ## [1] 1662.625 pchisq(chisq, (ncol(dat) - 1) * (nrow(dat) - 1), lower.tail=FALSE) # p値 ## [1] 0
ちなみに残差分析の調整済み標準化残差はこう。
(ASR <- Z / sqrt((1 - row.masses) %o% (1 - column.masses)) * sqrt(n)) #調整済み標準化残差 ## 犯罪 ## 地域 強盗 詐欺 破壊 ## オスロ -18.686816 37.35695 -22.453905 ## 中部 -3.164442 -18.32454 19.842422 ## 北部 22.931313 -26.50958 9.010296
調整済み標準化残差からp値を出して、さらに多重比較の調整。
# 残差分析のｐ値を算出 p.value.matrix <- ASR |> abs() |> pnorm(lower.tail=FALSE) * 2 round(p.value.matrix, 4) ## 犯罪 ## 地域 強盗 詐欺 破壊 ## オスロ 0.0000 0 0 ## 中部 0.0016 0 0 ## 北部 0.0000 0 0 # ホルムの方法による有意水準の調整（多重比較） Dim <- dim(p.value.matrix) p.value.matrix.holm <- matrix(p.adjust(p.value.matrix), Dim[1], Dim[2]) dimnames(p.value.matrix.holm) <- dimnames(dat) round(p.value.matrix.holm, 4) ## 犯罪 ## 地域 強盗 詐欺 破壊 ## オスロ 0.0000 0 0 ## 中部 0.0016 0 0 ## 北部 0.0000 0 0
便利ですね。
座標の重みづけ
ここからが本番です。
特異値分解で得られた特異ベクトルに様々な考え方で重みづけを行い同時布置図を描きます。
同時散布図の描き方には例えば以下のようなものがあります。
| 名前 | 座標 | 行と列の関係 |
|---|---|---|
| 対称プロット（フレンチプロット） | 主座標 | 角度△ |
| 非対称プロット | 主座標と標準座標 | 角度○ |
| 対称バイプロット | 標準座標×特異値の平方根 | 角度○ |
tarotanさんの記事、特に「3. さまざまな座標，さまざまな解釈」の節も参照してください。
ここから座標に重み付けをします。
標準座標
行（または列）のベクトルを行（または列）の質量の平方根で割る。
(standard.coordinates.rows <- sweep(U, 1, sqrt(row.masses), "/")) ## [,1] [,2] ## オスロ -0.8800182 0.05636457 ## 中部 0.9968363 -2.17806838 ## 北部 1.2009506 1.00876133 (standard.coordinates.columns <- sweep(V, 1, sqrt(column.masses), "/")) ## [,1] [,2] ## 強盗 1.2146096 2.0379805 ## 詐欺 -1.3067231 0.2881079 ## 破壊 0.5807713 -0.8368139
主座標
標準座標に特異値を乗じたもの。
(principal.coordinates.rows <- sweep(standard.coordinates.rows, 2, res.SVD$d[1:2], "*")) ## [,1] [,2] ## オスロ -0.3706976 0.008999028 ## 中部 0.4199058 -0.347744988 ## 北部 0.5058866 0.161056329 (principal.coordinates.columns <- sweep(standard.coordinates.columns, 2, res.SVD$d[1:2], "*")) ## [,1] [,2] ## 強盗 0.5116403 0.3253789 ## 詐欺 -0.5504421 0.0459986 ## 破壊 0.2446432 -0.1336036
対称バイプロット
標準座標に特異値の平方根を乗じたもの。
(sympc.coordinates.rows <- sweep(standard.coordinates.rows, 2, sqrt(res.SVD$d[1:2]), "*")) ## [,1] [,2] ## オスロ -0.5711573 0.02252169 ## 中部 0.6469756 -0.87029441 ## 北部 0.7794516 0.40307245 (sympc.coordinates.columns <- sweep(standard.coordinates.columns, 2, sqrt(res.SVD$d[1:2]), "*")) ## [,1] [,2] ## 強盗 0.7883167 0.8143192 ## 詐欺 -0.8481010 0.1151198 ## 破壊 0.3769374 -0.3343671
座標の組み合わせ
これらをもとに、3種類の同時布置図を描いてみましょう。
- 対称プロット（フレンチプロット）
- 非対称プロット
- 対称バイプロット
library(ca) library(factoextra) library(patchwork) res.ca <- ca(dat) p1 <- fviz_ca(res.ca, map = "symmetric", arrows = c(TRUE, TRUE), title="フレンチプロット") + scale_y_continuous(limits = c(-0.6, 0.7)) + scale_x_continuous(limits = c(-0.6, 0.7)) p2 <- fviz_ca(res.ca, map = "colprincipal", arrows = c(TRUE, TRUE), title="非対称プロット（列主座標）") + scale_y_continuous(limits = c(-2.5, 1.5)) + scale_x_continuous(limits = c(-1.75, 2.25)) p3 <- fviz_ca(res.ca, map = "symbiplot",arrows = c(TRUE, TRUE), title="対称バイプロット") + scale_y_continuous(limits = c(-1.0, 1.0)) + scale_x_continuous(limits = c(-1.0, 1.0)) p1 + p2+ p3 + plot_layout(ncol = 3)
非対称プロットでは列が原点近くに集まりすぎてちょっと読み取りにくいですね。
フレンチプロットと対称バイプロットはかなり似ていますが、対称バイプロットの方が縦方向のバラツキが大きくなっています。
指標化残差が角度と長さで表現できてるか確認
問題は指標化残差が角度としてちゃんと反映されているかどうかです。
指標化残差を再確認しておきます。
I ## 犯罪 ## 地域 強盗 詐欺 破壊 ## オスロ -0.4319130 0.4869917 -0.2228210 ## 中部 -0.1986758 -0.6488887 0.5348671 ## 北部 0.9426844 -0.6146521 0.1590303
角度で正しく表現できていれば、行の点Aの座標を(x1, y1)、列の点Bの座標を(x2, y2)とすると、x1 * x2 + y1 * y2が指標化残差と一致するはずです。
一致の度合いは相関係数で確認してみましょう。
フレンチプロット（主座標）
相関係数0.94。一致はしていませんが、「だいたい合ってる」とは言っても良さそう？
COL <- principal.coordinates.columns ROW <- principal.coordinates.rows # ここはもう少し効率的なやり方があるはず…… mat1 <- cbind( ROW[c(1,1,1,2,2,2,3,3,3), ], COL[c(1,2,3,1,2,3,1,2,3), ] ) mat2 <- expand.grid(cn = rownames(COL), rn = rownames(ROW)) row.names(mat1) <- paste(mat2[,2],mat2[,1]) # 行と列の内積 Ip <- mat1[,1] * mat1[,3] + mat1[,2] * mat1[,4] # 指標化残差のベクトル化 Iv <- I |> t() |> as.vector() cor(Ip, Iv) ## [1] 0.9444723
非対称プロット（列主座標）
相関係数1.00。完全に一致します。
COL <- standard.coordinates.columns ROW <- principal.coordinates.rows mat1 <- cbind( ROW[c(1,1,1,2,2,2,3,3,3), ], COL[c(1,2,3,1,2,3,1,2,3), ] ) mat2 <- expand.grid(cn = rownames(COL), rn = rownames(ROW)) row.names(mat1) <- paste(mat2[,2],mat2[,1]) # 行と列の内積 Ip <- mat1[,1] * mat1[,3] + mat1[,2] * mat1[,4] # 指標化残差のベクトル化 Iv <- I |> t() |> as.vector() cor(Ip, Iv) ## [1] 1
対称バイプロット
こちらも相関係数1.00。完全に一致します。
COL <- sympc.coordinates.columns ROW <- sympc.coordinates.rows mat1 <- cbind( ROW[c(1,1,1,2,2,2,3,3,3), ], COL[c(1,2,3,1,2,3,1,2,3), ] ) mat2 <- expand.grid(cn = rownames(COL), rn = rownames(ROW)) row.names(mat1) <- paste(mat2[,2],mat2[,1]) # 行と列の内積 Ip <- mat1[,1] * mat1[,3] + mat1[,2] * mat1[,4] # 指標化残差のベクトル化 Iv <- I |> t() |> as.vector() cor(Ip, Iv) ## [1] 1
大きな集計表で再確認
tarotanさんの記事には
カテゴリー数が少ないのであれば，わざわざ単純対応分析をする必要はないと思う
と書かれています。
それでは大きな集計表で再確認してみましょう。
日経リサーチさんの解説のデータを引用させてもらいます。
18ブランド、20属性の集計表です。
ブランド名はアルファベットに置き換えました。
dat.nr <- c( 18, 43, 175, 249, 52, 133, 78, 34, 93, 40, 207, 165, 108, 35, 78, 23, 115, 97, 30, 52, 124, 248, 59, 29, 48, 19, 113, 75, 165, 161, 24, 25, 47, 177, 36, 126, 16, 38, 21, 11, 36, 48, 101, 97, 86, 63, 50, 57, 50, 54, 58, 35, 43, 63, 42, 45, 37, 53, 32, 49, 51, 284, 26, 15, 7, 6, 52, 11, 110, 144, 10, 17, 27, 118, 20, 87, 8, 9, 4, 10, 284, 7, 452, 188, 76, 25, 177, 223, 181, 19, 156, 39, 272, 86, 52, 85, 28, 85, 185, 318, 60, 2, 265, 100, 72, 55, 64, 66, 33, 10, 115, 45, 111, 20, 32, 12, 39, 48, 62, 123, 85, 4, 329, 132, 79, 51, 94, 105, 67, 27, 118, 32, 121, 47, 32, 42, 25, 72, 94, 182, 189, 40, 233, 36, 27, 16, 80, 176, 107, 19, 57, 24, 88, 73, 13, 67, 7, 49, 89, 104, 106, 112, 87, 20, 16, 11, 63, 60, 55, 65, 24, 10, 43, 69, 10, 43, 8, 31, 38, 26, 37, 297, 31, 23, 14, 14, 60, 20, 102, 94, 24, 14, 19, 87, 26, 71, 4, 21, 6, 10, 148, 28, 197, 46, 25, 21, 70, 151, 80, 10, 78, 22, 67, 55, 16, 53, 12, 43, 73, 73, 77, 100, 159, 81, 32, 26, 118, 125, 126, 22, 110, 25, 81, 77, 29, 78, 15, 79, 48, 45, 352, 40, 253, 18, 24, 13, 125, 203, 164, 111, 30, 22, 177, 184, 12, 98, 14, 79, 119, 152, 24, 240, 58, 120, 44, 117, 95, 18, 125, 87, 76, 88, 69, 55, 96, 67, 161, 61, 17, 24, 42, 16, 149, 42, 41, 38, 30, 37, 13, 8, 58, 30, 49, 18, 10, 12, 28, 29, 28, 43, 168, 121, 140, 75, 50, 35, 140, 97, 184, 98, 39, 25, 88, 157, 30, 153, 57, 68, 64, 77, 72, 155, 106, 43, 25, 27, 64, 42, 113, 59, 49, 19, 48, 77, 13, 71, 18, 48, 39, 36, 49, 155, 74, 24, 34, 27, 59, 35, 64, 63, 27, 17, 28, 66, 20, 47, 18, 30, 23, 18 ) |> matrix( nrow = 18, byrow = TRUE, dimnames = list( Brands = LETTERS[1:18], Attributes = c( "伝統や格式を重んじる", "庶民的", "価格が高い", "流行を作り出している", "価値上昇", "新鮮味", "顧客の心つかんでいる", "気品を感じさせる", "ファン層幅広い", "機能性重視", "デザイン重視", "既成概念囚れない", "存在感", "ずっと付合っていける", "社会動きに敏感", "誰からも好かれる", "躍動感", "さりげなく特徴をアピール", "豊かな気持ちにさせる", "自慢になる" ) ) )
同時布置図ですが、カテゴリ数が多いので矢印は行（ブランド）のみに引いています。
res.ca <- ca(dat.nr) p1 <- fviz_ca(res.ca, map = "symmetric", arrows = c(TRUE, FALSE), title="フレンチプロット") + scale_y_continuous(limits = c(-1.25, 0.75)) + scale_x_continuous(limits = c(-1.25, 0.75)) p2 <- fviz_ca(res.ca, map = "colprincipal", arrows = c(TRUE, FALSE), title="非対称プロット（列主座標）") + scale_y_continuous(limits = c(-2.5, 1.5)) + scale_x_continuous(limits = c(-2.5, 1.5)) p3 <- fviz_ca(res.ca, map = "symbiplot",arrows = c(TRUE, FALSE), title="対称バイプロット") + scale_y_continuous(limits = c(-1.75, 1.0)) + scale_x_continuous(limits = c(-1.75, 1.0)) p1 + p2 + p3 + plot_layout(ncol = 2) + plot_annotation(title = "同時布置図の座標比較")
結論（再掲）
- コレスポンデンス分析の同時布置図を描くときは、対称バイプロットがおすすめです。
- 指標化残差を正確に角度として表現できて、なおかつ見やすいので。
- 従来の同時布置図はフレンチプロットが多いと思いますが、正確でないのであえて選ぶ理由はありません。
- ただ、「だいたい合ってる」ので、角度で読む分には大きく間違えることはなさそうです。
- 縦横のスケールを合わせるのが大前提です。
- そうしないと見かけ上の角度が歪んでしまいます。
参考リンク
Understanding the Math of Correspondence Analysis
Normalization and Scaling in Correspondence Analysis
対応分析のグラフを適切に解釈する条件 : Standard Coordinate, Principal Coordinate を理解する
