はじめに
この記事はStan Advent calendar12月12日の記事です。
普段よりもこのページを見る人が多いと思います、私のことはこちらで紹介しています。匿名主義なので何の情報もありませんが…
ここでは多変量解析の一つである対応分析を題材にします。
対応分析は多変量のカテゴリカルデータに用いられる分析で、クロス集計表の結果を視覚的に表現するために使われる統計的手法です。
実は対応分析は私の修士研究で使用した解析手法なのですが、対応分析の付置はデータ数に依存しないことから、座標を得るだけでその信頼性が確認できないところに当時問題を感じていました。
大学を卒業してからベイズ統計に関心を持ち、ある程度勉強したところで当時のことを思い出し、rstanを使ってこの問題に取り組もうと思ったのが本記事の執筆動機です。
本記事の構成は以下のとおりです。
参考にした論文は前回と同様のこちらです。
この論文はもうだいぶ古いですが、既に対応分析の付置の信頼性の問題に着目しています。ただし解析的な正規分布への近似をもとに信頼区間を求めているようなので、本記事ではこの論文とは別のベイズ的なアプローチをとることにします。
ベイズ的なアプローチの部分は先行文献が見当たらなかったのでオリジナルです。そんな変なことはしていないと思うのですが…いや、どうなんだろう。
対応分析について
対応分析とはクロス集計表など、カテゴリカルな行と列のデータを視覚的に表現することで、項目間の関係を把握しようとする手法です。コレスポンデンス分析とも呼ばれ、数理的には林先生の数量化理論第Ⅲ類と同等のものとされています。
対応分析の数理的な内容については、前回の記事で説明しているのでそちらを覗いてください。
対応分析の同時付置図は信頼できるのか
ここでは対応分析の同時付置図の信頼性について考えてみます。
付置は標本数に依らない？
前回の記事を振り返り、対応分析について復習します。
$i$行$j$列が$f_{ij}$（$i=1 … M$、$j = 1 … N$）である$M$行$N$列の分割表について考えたとき、対応分析の同時付置図作成までの過程は以下のようになっていました。
ここで各行・列および全項目の観測数の総和を前回と同様以下のように表記します。
-
確率行列
$$ \boldsymbol{P} _ I = \mathrm{diag} \left[ p_{1.},\ldots, p_{M.} \right] ~~~ p _{i.} = f _{i.}/Sum $$ $$ \boldsymbol{P} _ J = \mathrm{diag} \left[ p_{.1}, …, p_{.N} \right]~~~p_{.j} = f_{.j}/Sum $$ $$ \boldsymbol{P} _ {IJ}=\left(f_{ij}/N\right) $$より$M$行$N$列の行列$\boldsymbol{X}$$$ \boldsymbol{X} = \boldsymbol{P}_{I}^{-1} \boldsymbol{P}_{IJ} \boldsymbol{P}_{J}^{-1/2} $$を作成する。 -
$\boldsymbol{X}$の主成分分析を計算し、最大固有値から順に$\lambda_{1},\lambda_{2}\ldots,\lambda_{N-1}$、およびこれらに対応する固有ベクトル$\boldsymbol{l} _1,\ldots,\boldsymbol{l} _{K}$ ($K = \rm{min}(\textit{M},\textit{N})-1$)を得る1。
-
第$k$主成分ベクトルによる項目$I$の数量化スコアを下記式より計算する。
$$ \boldsymbol{z} _k = \boldsymbol{X} \boldsymbol{l} _k = \boldsymbol{P} _I^{-1} \boldsymbol{P} _{IJ} \boldsymbol{P} _{J}^{-1/2}\boldsymbol{l} _k \tag{1} $$ -
項目$J$の数量化スコアを、下記式に基づき項目$I$の数量化スコアから求める。
$$ \boldsymbol{z} _k^* = \cfrac{1}{\sqrt{λ_k}} \boldsymbol{P}_J^{-1} \boldsymbol{P}_{IJ}\boldsymbol{z} _k \tag{2} $$ -
$\boldsymbol{z}_k$、$\boldsymbol{z}^{*}_k$をもとに行項目・列項目の同時付置図を作成する。
ここで、$(1)$$(2)$式を見ると、すべての数量化スコアが確率行列のみに依存しており、標本数$Sum$に影響されないことが確認できます。つまり、標本数が違う2つの分割表があったとして、それらの確率行列$P_{IJ}$が一致する場合、同時付置図は全く同じものになると…
よって、同時付置図はそれのみでは付置がどの程度信頼できるかが確認できないようになっているのです。
分割表の一様性の検定
上記問題の解決というわけではないですが、分割表の一様性の検定というものがありまして、分割表の各行の確生起確率$(p_{i1},\ldots,p_{iN})$が同じであるといえるかどうかを判定する検定があります。分割表の検定の結果下に述べる帰無仮説が棄却された場合、存在すると考えられる生起確率の差の傾向を対応分析で明瞭に把握できるだろう、というわけです2。
この検定の帰無仮説は
です。$T$が漸近的に自由度$(M-1)(N-1)$の$χ^2$分布に従うことを用いて検定を行います。
試しに前回の記事と同じ相撲データについて一様性の検定をしてみましょう。
# smo_dataの作成 ----------------------------------------------------------
smo_data <- matrix(data=c(15,15,0,0,1,0,1,19,9,7,12,1,3,3,4,22,0,11,2,0,4,
14,3,11,3,1,1,6,24,7,5,5,0,0,2,30,2,6,2,2,0,2,
11,13,5,6,2,1,4,25,8,10,7,5,3,1,13,31,3,3,0,2,3,
8,2,19,3,2,0,4,28,8,8,2,1,0,7,9,3,14,11,2,0,4,
5,3,12,11,4,0,7,7,1,18,0,1,0,11,12,6,1,12,7,2,2,
25,2,11,3,1,0,1,16,4,18,3,0,0,5,23,2,6,5,1,0,2,
7,17,2,16,3,0,1,2,8,1,11,1,0,2),byrow=TRUE, ncol=7)
colnames(smo_data) <- c("寄り","押し","投げ","引き","送り","突き","その他")
rownames(smo_data) <- c("小錦","琴錦","貴闘力","霧島","栃ノ和歌","水戸錦","若ノ花","貴ノ花","武蔵丸",
"大翔山","安芸ノ島","三杉里","旭道山","舞の海","寺尾","琴の若","貴ノ浪",
"琴富士","隆三杉","春日富士")
# smo_dataの一様性検定 ----------------------------------------------------------
my_uniformity_test <- function(data){
data_teststat <- matrix(NA, nrow=nrow(data), ncol=ncol(data))
p_j <- rep(NA, length=n)
for(j in 1:n){
p_j[j] <- sum(data[,j])/sum(data)
}
for(i in 1:m){
for(j in 1:n){
data_teststat[i,j] <- (data[i,j] - rowSums(data)[i]*p_j[j])^2 / (rowSums(data)[i]*p_j[j])
}
}
T_stat <- sum(data_teststat)
cat(sprintf("p-value:%.5f",pchisq(T_stat, df=((m-1)*(n-1)), lower.tail = FALSE)))
}
my_uniformity_test(smo_data)
## p-value:0.00000
帰無仮説は有意水準0.1%以下で棄却されました。では、このデータの各度数を4分の1にしたデータではどうでしょうか4。
my_uniformity_test(smo_data/4)
## p-value:0.39973
p値は約0.4ですので、とても帰無仮説を棄却できません。
そんな訳なので、同じ付置図が得られても標本数によって信頼性に差があるため、付置図だけからは正しい解釈ができない可能性があるんです。
同時付置図と一様性の帰無仮説の関係
対応分析の同時付置図がどのようになっているときに、一様性の帰無仮説が棄却されやすいかを考えます。
生起確率が全体の平均に等しく$(p_{.1},\ldots,p_{.n})$である個体$i$の第$k$数量化スコアは、
なので5、同時付置図の原点が帰無仮説に対応する点になっています。
また、$(3)$式を展開すると、
となることから、帰無仮説の点（同時付置図の原点）から離れている行$i$で、行和$f_i$が大きいものがある場合に$T$が大きくなり、一様性の帰無仮説が棄却されやすい6ことが分かります。
以上を踏まえ、同時付置図にそれぞれのポイントの信頼度を表示することができれば、分割表の情報を余すことなく付置図に表現でき、視覚的な結果の解釈や項目の分類の手助けになるのでは、というのが本記事の狙いです。
モデル作成
分割表の生成過程にモデルを設定して、MCMCによってそのパラメータを推定し、同時付置図における各ポイントのベイズ信頼区間（またはベイズ予測区間）を求めることを考えます7。
分割表のモデリング
分割表の生成過程は、分割表の性質により以下の2パターンのいずれかを仮定することができます。
対応分析にかける分割表の場合、後者の場合が多いかと思います。一様性の検定の際にも、後者の仮定を設定して、検定をしていました。
$(6)$式とデータによる尤度に基づいたMCMCにより、各行ごとの観測度数の生成確率$\boldsymbol{\theta} _ {i}$の事後分布及びMCMCサンプルが得られます。
観測度数の行の総和$F_{i.}$が小さい行の$\boldsymbol{\theta}_{i}$ほど、事後分布の推定幅が広くなります。このMCMCサンプルに対する主成分ベクトルの射影を得ることで、同時付置図における項目$I$（行項目）の各ポイントの信頼区間を表示できると考えられます。
項目$J$（列項目）についても、各列の観測度数に独立の多項分布を仮定することで、各ポイントの座標について信頼区間が得られると考えられますが、今回の実装では項目$J$を変量とみなし、列毎の観測度数の総和は試行ごとに一定でないと考えるので、列毎に独立の多項分布を仮定しない（列項目のポイントの座標の信頼区間は求めない）ことにします8。
事前分布
$\boldsymbol{\theta}_i$について以下の事前分布を設定しています。これは、相撲データの性質を踏まえて設定したもので、決まり手には比較的メジャーで頻繁に使われる技もあれば、めったに決まる、もしくは使われることのない技もあるだろう、ということを考慮しました。
$(9)$式のように事前分布を設定すると、ディリクレ分布の性質により$\theta_{ij}$の期待値と分散は以下のようになります。
なので、$(10)$式より各行毎のセルの値の生成確率$\theta_{ij}$の期待値は全観測度数に占める列毎の観測度数の総和の割合$\cfrac{f_{ij}}{f_{i.}}$に等しくなります。
分散については深く考えて決めたものではありませんが、ディリクレ分布の全パラメータの比はそのままに和を大きくしていくと、出力の期待値は一定のまま分散が減少するようになっているので、この性質を利用して$(9)$式の$(p_{.1},\ldots,p_{.J})$を等倍していって事前分布の強さを変化させていっても面白そうです。
MCMCサンプルより生成量を得る値
対応分析の分析過程のなかのどの数値に対し、MCMCサンプルより生成量を得るかを考えます。$(1)$式より、座標は$\boldsymbol{X}$及び固有ベクトル$\boldsymbol{l}_k$ により計算されますから、それぞれについて考えます。
主成分分析の対象$\boldsymbol{X}$
行列$\boldsymbol{X}$は、統計的に独立な$M$つの多項分布からの生成量より構成されるものです。
これの生成量を素直に得ようとすると、$M$個の独立な生成過程がかかわってしまいます。考え方は色々あると思いますが、ここでは、個体$i$の同時付置図の信頼区間には、個体$i$の従う多項分布のパラメータ$\boldsymbol{\theta}_i$の影響しか反映させないものとし、$(12)$の行列から$\boldsymbol{X}$を計算することにします9。
固有ベクトル$\boldsymbol{l}_k$
分割表の生成過程のモデリングにより、各セルの観測度数の期待値がMCMCサンプルから計算できます。それらをもとに新しい分割表を生成し、新たに主成分ベクトルを得ることができるにはできるのですが、そうすると軸がMCMCサンプル毎に異なってしまいます。最終的に作成する同時付置図は軸を固定して示したいのでこれは不都合です。
今回のモデルでは、観測値を真の値とみなすのと同様に、観測値から得られた主成分ベクトルを真のベクトルであるとみなし、同時付置図も観測値から計算された主成分ベクトル上で付置するものとします。固有値や寄与率についても同様に観測値から計算された値を真の値とします。
モデルの実装
以上のモデルを実装したStanコードを載せます。
data{}ブロックでは分割表に関する数値を指定するだけですが、通常の対応分析を実行するパートはStan外で行うことにしたので、事前に得られた固有ベクトルをvector[N] eigenvector[K]で指定しています。
transformed data{}ブロックでは、後のgenerated quantities{}ブロックでの計算に必要な値$p_{i.}$、$P_{i.}$や事前分布の計算に必要な$p_{.j}$を指定します。
ちなみにStanで実装されている関数eigenvalues_sym()やeigenvectors_sym()を使えば、このブロック内で対応分析を済ませることも可能です。私はStan内で主成分ベクトルを計算すると、外で計算したものと異なってくるのでクレバーではないと判断しました。
{parameters{}ブロックでは項目$I$について仮定した各行毎に独立の多項分布のパラメータ$\boldsymbol{\theta_{i}}$を指定しています。
model{}ブロックでは$(6)$式のモデル及び$(9)$式の事前分布を指定しています。
generated quantities{}ブロックでは、data{}ブロックで指定した固有ベクトルやtransformed data{}で指定した行列・ベクトル、及びMCMCサンプルを用いて座標の生成量を求めています。項目$I$のポイントについては2通り計算しており、観測度数の期待値から分割表を作成し求めた各ポイントの座標の生成量（Z_EST）と、観測度数の期待値から乱数を発生させて作成した分割表から求めた座標の生成量（Z_RNG）があります。
// model.stan
data {
int<lower=0> Sum; //全観測数の和
int<lower=0> N; //分割表の列数
int<lower=N> M; // 分割表の行数
int<lower=0> K; // min(N,M)-1 対応分析で得られる次元数
int d[M,N]; //分割表
vector[N] eigenvector[K];
}
transformed data{
matrix[M,N] d_numeric; //分割表 確率行列を作成するためにint型の分割表は使えない
vector[M] p_i; //行ごとの観測値全体に占める割合
vector[N] p_j; //列ごとの観測値全体に占める割合
matrix[M,M] P_I;
for(j in 1:N){
for(i in 1:M){
d_numeric[i,j] = d[i,j];
}
}
for(i in 1:M){
p_i[i] = sum(d_numeric[i,]) / Sum;
}
for(j in 1:N){
p_j[j] = sum(d_numeric[,j])/Sum;
}
P_I = diag_matrix(p_i);
}
parameters{
simplex[N] theta[M];
}
model {
//モデル部分
for(i in 1:M){
d[i,] ~ multinomial(theta[i]);
}
//事前分布
for(i in 1:M){
theta[i] ~ dirichlet(p_j);
}
}
generated quantities{
vector[M] Z_EST[K]; //求めたい座標（期待値）
vector[M] Z_RNG[K]; //求めたい座標（事後分布からの乱数生成）
//期待値の座標を計算
{
matrix[M,N] P_IJ_EST;
matrix[M,N] P_IJ_RNG;
int d_RNG[M,N];
matrix[M,N] d_numeric_RNG;
matrix[N,N] sq_P_J_EST;
matrix[N,N] sq_P_J_RNG;
matrix[M,N] X_EST;
matrix[M,N] X_RNG;
//以降は行毎に必要な操作
for(i in 1:M){
//着目する行によってP_IJは異なる
for(m in 1:M){
//P_IJの生成
if(m==i){
P_IJ_EST[m,] = (theta[i] * (sum(d_numeric[i,])) / Sum)'; //対象の行のみ期待値を用いる
}else{
P_IJ_EST[m,] = d_numeric[m,] / Sum;
}
}
//着目する行によってP_Jも異なる
{
vector[N] sq_p_j_EST;
for(n in 1:N){
sq_p_j_EST[n] = sqrt(sum(P_IJ_EST[,n]));
}
sq_P_J_EST = diag_matrix(sq_p_j_EST);
}
X_EST = inverse(P_I) * P_IJ_EST * inverse(sq_P_J_EST);
for(k in 1:K){
Z_EST[k,i] = X_EST[i] * eigenvector[k];
}
}
//乱数生成したときの座標を計算
for(i in 1:M){
for(m in 1:M){
//P_IJの生成
if(m==i){
d_RNG[m,] = multinomial_rng(theta[m],sum(d[m,]));
}else{
d_RNG[m,] = d[m,];
}
}
for(n in 1:N){
for(m in 1:M){
d_numeric_RNG[m,n] = d_RNG[m,n];
}
}
P_IJ_RNG = d_numeric_RNG / Sum;
{
vector[N] sq_p_j_RNG;
for(n in 1:N){
sq_p_j_RNG[n] = sqrt(sum(P_IJ_RNG[,n]));
}
sq_P_J_RNG = diag_matrix(sq_p_j_RNG);
}
X_RNG = inverse(P_I) * P_IJ_RNG * inverse(sq_P_J_RNG);
for(k in 1:K){
Z_RNG[k,i] = X_RNG[i] * eigenvector[k];
}
}
}
}
結果の確認
以下でMCMCを走らせ、$\hat{R}$で収束を確認します。モデル部分が簡素なのですぐ終わりますし、収束もばっちりですな。
library(rstan)
# 行列Xの作成 ------------------------------------------------------------
N <- sum(smo_data)
m <- nrow(smo_data)
n <- ncol(smo_data)
K <- min(n,m) - 1
p_i <- rep(NA, length=m)
for(i in 1:m){
p_i[i] <- sum(smo_data[i,])/N
}
P_i <- diag(p_i)
p_j <- rep(NA, length=n)
for(j in 1:n){
p_j[j] <- sum(smo_data[,j])/N
}
P_j <- diag(p_j)
P_ij <- smo_data/N
matpow <- function(x, pow=2) {
y <- eigen(x)
y$vectors %*% diag( (y$values)^pow ) %*% t(y$vectors)
}
X <- solve(P_i) %*% P_ij %*% matpow(P_j, pow=(-1/2))
# 論文の方法に則った対応分析 -----------------------------------------------------------
X_bar <- matrix(rep(colSums(X)/nrow(X),nrow(X)), byrow=T, nrow=nrow(X))
X_cov <- t((X - X_bar)) %*% P_i %*% (X - X_bar)
X_pca <- eigen(X_cov)
rikishi_coord <- (X %*% X_pca$vectors[,1:K])
waza_coord <- matrix(NA, nrow=n, ncol=K)
for(k in 1:K){
waza_coord[,k] <- 1/sqrt(X_pca$values[k]) * solve(P_j) %*% t(P_ij) %*% rikishi_coord[,k]
}
waza_coord <- data.frame(dim=waza_coord, target=c("寄り","押し","投げ","引き","送り","突き","その他"))
varsum <- sum(X_pca$values[1:6])
cat(paste(sprintf("\nContribution Rate of dim%.f: %.3f", 1:6,X_pca$values[1:6]/varsum)))
## Contribution Rate of dim1: 0.463
## Contribution Rate of dim2: 0.248
## Contribution Rate of dim3: 0.165
## Contribution Rate of dim4: 0.060
## Contribution Rate of dim5: 0.038
## Contribution Rate of dim6: 0.026
var_rate <- X_pca$values[1:6]/varsum
# stanを走らせる ------------------------------------------------------------
data <- list(Sum=N, N=n, M=m, d=smo_data, K=K, eigenvector=t(X_pca$vectors[,1:K]))
library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
fit <- stan(file="model.stan", data=data, chains = 4, iter=1300, warmup = 300)
print(fit, pars="theta")
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## theta[1,1] 0.47 0 0.08 0.30 0.41 0.47 0.53 0.63 6965 1
## theta[1,2] 0.46 0 0.08 0.30 0.40 0.46 0.52 0.62 6147 1
## theta[1,3] 0.01 0 0.01 0.00 0.00 0.00 0.00 0.04 4527 1
## theta[1,4] 0.00 0 0.01 0.00 0.00 0.00 0.00 0.03 4572 1
## theta[1,5] 0.03 0 0.03 0.00 0.01 0.02 0.04 0.11 3747 1
## theta[1,6] 0.00 0 0.00 0.00 0.00 0.00 0.00 0.00 4123 1
## theta[1,7] 0.03 0 0.03 0.00 0.01 0.02 0.05 0.11 4463 1
## ～（以下略、全部Rhat1でした）～
次にインスタ映えプロットを作ります。今回のインスタ映えプロットは観測値から得られる（真の）第1軸と第2軸をもとに付置した対応分析の同時付置図です。項目$I$の各ポイントについてそのベイズ信頼区間を示しました10。
library(tidyr)
library(dplyr)
library(stringr)
library(tibble)
library(ggplot2)
library(ggrepel)
# MAP 推定値を得る関数
MAP <- function(x){
dens <- density(x)
mode_i <- which.max(dens$y)
mode_x <- dens$x[mode_i]
return(mode_x)
}
# 集計 ----------------------------------------------------------------------
rikishi_label <- c("小錦","琴錦","貴闘力","霧島","栃ノ和歌","水戸錦","若ノ花","貴ノ花","武蔵丸",
"大翔山","安芸ノ島","三杉里","旭道山","舞の海","寺尾","琴の若","貴ノ浪",
"琴富士","隆三杉","春日富士")
fit %>% rstan::extract() %>% data.frame()%>%
select(starts_with("Z_EST")) %>% pivot_longer(everything(), names_to="key", values_to = "val") %>%
mutate(dim=as.numeric(str_sub(key, start=7,end=7)),
target=as.numeric(str_sub(key,start=9, end=10)),
key=NULL) %>% group_by(target, dim) %>% nest() %>%
pivot_wider(names_from = dim, values_from = data)%>% unnest() %>%
group_by(target) %>% rename(dim1=val, dim2=val1, dim3=val2, dim4=val3, dim5=val4, dim6=val5) %>%
summarise_all(funs(MAP=MAP(.),lower=quantile(., 0.25),upper=quantile(.,0.75))) %>% arrange(target) %>%
mutate(target=factor(seq(1,m),
labels=rikishi_label)) -> plot_row_EST.df
fit %>% rstan::extract() %>% data.frame() %>%
select(starts_with("Z_RNG")) %>% pivot_longer(everything(), names_to="key", values_to = "val") %>%
mutate(dim=as.numeric(str_sub(key, start=7,end=7)),
target=as.numeric(str_sub(key,start=9, end=10)),
key=NULL) %>% group_by(target, dim) %>% nest() %>%
pivot_wider(names_from = dim, values_from = data) %>% unnest() %>% group_by(target) %>%
rename(dim1=val, dim2=val1, dim3=val2, dim4=val3, dim5=val4, dim6=val5) %>% summarise_all(funs(MAP=MAP(.),lower=quantile(., 0.25),upper=quantile(.,0.75))) %>% arrange(target) %>%
mutate(target=factor(seq(1,m),
labels=rikishi_label)) -> plot_row_RNG.df
fit %>% rstan::extract() %>% data.frame() %>%
select(starts_with("Z_EST")) %>% pivot_longer(everything(), names_to="key", values_to = "val") %>%
mutate(dim=as.numeric(str_sub(key, start=7,end=7)),
target=str_sub(key,start=9, end=10),
key=NULL) %>% group_by(target, dim) %>% nest() %>%
pivot_wider(names_from = dim, values_from = data) %>% unnest() %>%
rename(dim1=val, dim2=val1, dim3=val2, dim4=val3, dim5=val4, dim6=val5) %>% ungroup() %>% mutate(target=factor(target,levels=seq(1,m), labels=rikishi_label))-> plot_row_EST_range.df
fit %>% rstan::extract() %>% data.frame() %>%
select(starts_with("Z_RNG")) %>% pivot_longer(everything(), names_to="key", values_to = "val") %>%
mutate(dim=as.numeric(str_sub(key, start=7,end=7)),
target=str_sub(key,start=9, end=10),
key=NULL) %>% group_by(target, dim) %>% nest() %>%
pivot_wider(names_from = dim, values_from = data) %>% unnest() %>%
rename(dim1=val, dim2=val1, dim3=val2, dim4=val3, dim5=val4, dim6=val5) %>% ungroup() %>% mutate(target=factor(target,levels=seq(1,m), labels=rikishi_label))-> plot_row_RNG_range.df
# Plot --------------------------------------------------------------------
library(RColorBrewer)
mycol <- c(brewer.pal(11,"Paired"),brewer.pal(9,"Set1"))
# Z_ESTを用いたプロット
p <- ggplot() +
theme_light(base_size=11) + guides(alpha=FALSE) +
stat_density_2d(data=plot_row_EST_range.df, aes(x=dim1, y=dim2, alpha=..nlevel.., fill=target), geom="polygon", bins=7)+
geom_errorbar(data=plot_row_EST.df, aes(x=dim1_MAP, y=dim2_MAP,ymin=dim2_lower, ymax=dim2_upper),col="grey30", width=0.1,alpha=1) +
geom_errorbarh(data=plot_row_EST.df, aes(y=dim2_MAP,xmin=dim1_lower, xmax=dim1_upper),col="grey30", height=0.1,alpha=1) +
geom_point(data=waza_axis, aes(x=dim.1, y=dim.2),size=2, shape=17) +
geom_point(data=plot_row_EST.df, aes(x=dim1_MAP, y=dim2_MAP,fill=target),shape=21,size=2) +
geom_text_repel(data=plot_row_EST.df,aes(x=dim1_MAP,y=dim2_MAP,label=target)) +
geom_text_repel(data=waza_axis, aes(x=dim.1, y=dim.2, label=target)) +
labs(title="Correspondence Analysis with Bayesian 50% confidence intervas & density",
x=sprintf("dim1(%.3f%%)", var_rate[1]*100), y=sprintf("dim2(%.3f%%)", var_rate[2]*100)) +
theme(legend.title = element_blank()) +
scale_fill_manual(values=mycol) +
scale_alpha_continuous(range=c(0.1,0.5))
p
.png)
今回は1次元・2次元だけでなく3次元以降も可視化できるようにしてます。
_dim3.png)
It’s so brilliant!
どちらの図も、付置が縄跳び状ですね。対応分析の付置図はこのようになわとびのような形になる（馬蹄形効果）ことが多いようです。それ自体が悪いことなのかよく知らないのですが、これに対する対応策もあるようです（今後の自分への宿題）。
対応分析の結果の解釈は前回の記事に譲り、項目$I$の各ポイントの信頼区間について確認してみます。
白星の最も多い力士は貴乃花で59、最も少ない力士は春日富士で25でした。そこで両者の付置の信頼区間をみると、春日富士の方がやはり貴乃花よりも推定幅が広いようです。がしかしその差は第1・第2主成分ベクトルの付置図ではそれほど明瞭でありません。主成分ベクトルによる写像の影響も絡んでいるのでしょう。
一応貴ノ花と春日富士の従う多項分布のパラメータの事後分布を確認します。
.png)
.png)
春日富士の方がやはり事後分布の幅が広い？
ほかの力士についても言及すると、第1・第2主成分ベクトルの付置図では小錦や武蔵丸などは他の力士と分布がほぼかぶっていないので、他の力士とは性格の異なる力士であるということが自信をもって言えそうです。逆に分布のかぶり具合で力士同士の類似性もわかる？
パラメータの事後分布から乱数生成して求めた各ポイントの座標の予測区間も示しておきます11。誤差を考慮したので当然着色部の幅が広くなっています。
.png)
おまけ～私が相撲に参加すると？～
標本数の違いによる信頼区間の幅の差があまり実感出来なかったので、極端な例で試してみました。
私(R.morta)が相撲に参加した体で、新たな行を作成します。ただ私は相撲などしたことがない一般人なので、当然力士に勝てる見込みもありません。それでは分割表が作れないので、八百長をした体にして、「引き」での勝ち数が多いことにしてみましょう12。
smo_data <- matrix(data=c(15,15,0,0,1,0,1,19,9,7,12,1,3,3,4,22,0,11,2,0,4,
14,3,11,3,1,1,6,24,7,5,5,0,0,2,30,2,6,2,2,0,2,
11,13,5,6,2,1,4,25,8,10,7,5,3,1,13,31,3,3,0,2,3,
8,2,19,3,2,0,4,28,8,8,2,1,0,7,9,3,14,11,2,0,4,
5,3,12,11,4,0,7,7,1,18,0,1,0,11,12,6,1,12,7,2,2,
25,2,11,3,1,0,1,16,4,18,3,0,0,5,23,2,6,5,1,0,2,
7,17,2,16,3,0,1,2,8,1,11,1,0,2,
3,0,1,5,0,0,2),byrow=TRUE, ncol=7)
colnames(smo_data) <- c("寄り","押し","投げ","引き","送り","突き","その他")
rownames(smo_data) <- c("小錦","琴錦","貴闘力","霧島","栃ノ和歌","水戸錦","若ノ花","貴ノ花","武蔵丸",
"大翔山","安芸ノ島","三杉里","旭道山","舞の海","寺尾","琴の若","貴ノ浪",
"琴富士","隆三杉","春日富士","R.morita")
～（以下略）～
.png)
私の従う多項分布のパラメータ事後分布幅はやはり広いですね。 インスタ映えプロットも見てみましょう。
_R.morita.png)
期待通りの結果です。私の信頼区間だけ他と比べて広いです。特に「引き」の寄与が大きい「組んで取る相撲」-「離れて取る相撲」の第2軸についてその傾向が大きく、座標の正負もあやしい状態です。
以上で対応分析の同時付置図に各ポイントの信頼区間を示すことが出来ていることを確認できました。
まとめ
対応分析の同時付置図にベイズ信頼区間を表示する試みを紹介しました。やったことは単純で、対応分析の同時付置図には抜け落ちてしまっている標本数や付置の信頼度の情報を、多項分布モデルの事後分布幅に読み替え、MCMCサンプルの分布を同時付置図上に加えただけです。
同時付置図を基に項目をパターン分けしたりする際、今回のような方法で同時付置図上にポイントの信頼区間を表示することで、誤った判断を下すことを防げるのではないかと思っています。
あとは図がとても美しいのがウリですね。
-
$N$個の列変量に対する次元削減を狙いとした主成分分析なので、$N$個の固有値・固有ベクトルが得られます。そのうち最小固有値は自明な解$\lambda_{N}≃0$、$l_{N}=(\sqrt{p_{.1}},\ldots,\sqrt{p_{.N}})$であり、それを除いた$N-1$次元が得られる次元になるようです。（論文では自明な解$\lambda_{N}≃1$と記載されていたが計算すると自明な固有ベクトルに対応する固有値は0に近かった。どちらが正しいのか検算・理解が追い付かなかった…）さらに$N>M$のときは$M-1$次元までを採用することになります ↩︎
-
行列$X$のユークリッド距離が$\cfrac{f_{ij}}{f_{i.}}$に関する行項目間の$χ^2$距離になっているので、分割表の一様性の検定における検定対象は対応分析において視覚的に強調しようとしている対象そのものになっています ↩︎
-
参考論文では帰無仮説：$\cfrac{p_{ij}}{p_{i.}} = p_{.j}$となっていますが同じことです。分かりにくいナ ↩︎
-
4で割ると整数にならない観測度数がでてきますが気にしないようにしましょう。整数じゃなくても検定はできます。 ↩︎
-
自明な固有ベクトルとの内積になっているので0です ↩︎
-
最右辺の$\cfrac{\left(\cfrac{p_{ij}}{p_{i.}}-p_{.j}\right)^2}{p_{.j}}$は、対応分析で各ポイントと原点からの距離として視覚的に表示される対象そのものです。分割表の一様性検定における検定統計量が$(5)$式のように分解されることからも、
対応分析の同時付置図には標本数の情報がすっぽり抜き取られていることが確認できます。 ↩︎ -
ベイズ信頼区間は、パラメータのMCMCサンプルより分割表の各セルの期待値を計算し、各ポイントの分布を求めたものとし、ベイズ予測区間はパラメータのMCMCサンプルから分割表の各セルの値を乱数生成した上で各ポイントの分布を求めたものとします ↩︎
-
列項目については$項目間のχ^2$距離を最大限説明できる軸を解析的に得るのみにして、生成過程を仮定してMCMCサンプルを得る、ということはしない、という考えです ↩︎
-
個体$i$の付置の予測区間を求める際は、$(12)$式の$i$行目に$\boldsymbol{\theta_{i}}$のMCMCサンプルを用いた乱数生成量をいれます ↩︎
-
エラーバーで各座標の50%信頼区間を示し、着色で二次元カーネル密度推定結果を、等高線の最大値を1としたときの0.2以上の範囲を描画しています。 ↩︎
-
Z_RNGをZ_ESTのかわりに用いることで得られます。集計までのコードも一緒に載せています ↩︎
-
八百長での決まり手には「引き」がよく採用されるようです ↩︎
