Within_Between
設定
Org. 2020.3.12
2021.1.3, 2.13, 2.16
授業での文面のままである。 但し、朝野氏の論説を2021.2.16に確認して、その検討とデータを最後に追加した。
library(FactoMineR)
library(anacor)
library(ca)
library(CAinterprTools)
library(corrplot)## corrplot 0.90 loadedlibrary(DescTools) # added
source("pcabi.R")
source("superd.R")
source("permutesEig.R")
source("bertin2.R")
source("assocplot3.R")
load("mexm.rdata")
load("asa12.rdata") # added
data(tocher)添付の mexm.rdata は、様々な実例からなる
- 西里 barley, rorshach
- 豊田 cosme, pastry
- 林 kanzume, kanzume0, kojidai
- 藤本（Clausen） norway (.tbl2.1)
- 矢野 (Greenacre 1984) kitsuen
- 大津 goodman.mental
- Haberman Hab1
- Guttman upper6, offdiag6
- Rdata children0, haireye super6
- その他 hashi, reconst
課題
対応分析には、様々な課題がある。参考のために添付したものは、西里氏によるもの以外に、作新学院大学 人間文化学部の 藤本一男教授による論説を添付した。日本語ではあるが、院生がすぐに理解するとは限らず、一方では出来の良い学生（定量・定性を受けていれば）理解するのもいるかと思われる。
有効な次元数
主成分分析の場合、平行分析が、色々意見はあるにせよ有効だと知られている。因子分析の場合もさらにMAP基準が知られている。対応分析の場合は、\(\chi^2\)値の分解に応じる Malinvaud’s test はある。また sig.dim.perm.scree は、一種の permutation test である。この２つは CAinterprTools で提供されている。
TIA
Total Information Analysis
既に説明したように、集計表の\(\chi^2\)距離は、principal coord.での布置のユークリッド距離である。それは、行変数座標と列変数座標各々について成立しており、within-variable あるいは intra-variable distance と呼ばれる。これに対して、行変数座標と列変数座標の間の距離は between-variable あるいは extra-variable distance と呼ばれ、その合理的な定め方は西里によって決定されている。その両方の距離を一括して扱う場合、全体を super-distance と呼び、それを用いる研究方法を TIA と称した。 Total Information Analysis である。
このTIAは、理論的に対応分析の全体を整合させているとはいえる。ところが、それを紹介している論説はあるにせよ、実際の状況をきちんと検討しているものが存在しない。西里氏自身は、極力低次元でのグラフ表示を使わないので、説明を解読するのは難しい。 さらには、実際に実例を色々チェックすると、固有値のバランスによって、かなり不思議な状況が発生することがわかる。
実際に色々確認してみよう。
このために、super-distance の計算と表示のスクリプト superd.R を作成した。その 70行目あたりからの spd が距離であり、その他の表示プログラムもすべて入れこんである。singv が対応分析の処理での特異値である。学生にも分かりやすいように、apply系は用いていない。
なお、superd によるbarley と rorshach の距離行列は、Clavel-Nishisato (2010) p.95 Table 1, p.97 Table 4 の下三角部分と各々一致する。また、それのMDSによる固有値も、論文の記述と一致する。以下の barleyst, rorshach のグラフにても読み取れる。
ところが、西里の一部も含めて紹介している 藤本一男(2017) にあるデータ norway（犯罪件数）については、多少雰囲気の異なる結果を得る。つまり、行と列別に直線に乗る。場合によってはそれがかなり極端になる場合もある。解釈に適切かどうかはまた別途考えねばならないが、同様の例は多く見いだせる。
norway data の Bertin’s plot
ここでは解釈を述べないが、藤本氏が mosaic.plot を提供しておられるので、比較のために Bertin’s plot(大津氏によるスクリプト)と、assocplot を用いた 調整済残差の表示と、同じデータの corrplot を用いた表示を行う。
par(mfrow=c(1,2), mar=c(0,0,1,0))
bertin.q(norway)
bertin.q(norway, bardirection="v")
bns <- assocplot3(norway)
corrplot(t(bns$stdres), is.corr=FALSE, method="ellipse")
実例について
以下では、屡々使われる実例について、次のグラフを提示する。
- Malinvaud’s test
- permutation test
- CA eigenvalues
- MDS(cmdscale) eig. super-dist.
- TIA mds-plot
- CA biplot (principal coord.)
そして、TIA mds-plot においては、行変数（青）と列変数 （赤）の各々について、その２次元座標での回帰直線を追加してある。
西里氏は、このような低次元のグラフ表示は好まれず、kmeans や クラスター分析が好ましいとされる。しかし、ここではあえて mds のプロットを提示している。
CA eigenvalues の棒グラフでは、Malinvaud’s test の結果で 有意でない となったものを、濃い色にしている。
類似した結果の場合
CA の Dim1, Dim2(, Dim3) が比較的似た寄与率を持っている場合、CAの biplot と super-distance でのMDSの配置は比較的似ている。この場合、残りの次元は寄与が小さい。
par(mfrow=c(1,1), mar=c(2.5,2.5,2,0.5), mgp=c(1.5,0.5,0))
superd(barleyst)##
## 次のパッケージを付け加えます: 'plotrix'## 以下のオブジェクトは 'package:CAinterprTools' からマスクされています:
##
## rescale
superd(rorshach)## Warning in chisq.test(x): Chi-squared approximation may be incorrect## Warning in chisq.test(t(x)): Chi-squared approximation may be incorrect
superd(hashi0)## Warning in chisq.test(x): Chi-squared approximation may be incorrect
## Warning in chisq.test(x): Chi-squared approximation may be incorrect
superd(kanzume0)## Warning in chisq.test(x): Chi-squared approximation may be incorrect
## Warning in chisq.test(x): Chi-squared approximation may be incorrect
superd(offdiag6)## Warning in chisq.test(x): Chi-squared approximation may be incorrect
## Warning in chisq.test(x): Chi-squared approximation may be incorrect
かなり異なる見かけとなる場合
CAで Dim1 が比較的優勢だが、高次元の寄与がだらだらとあって総計ではかなりになｒる場合、列変数と行変数が別々に直線上に乗る状況が生じる。
これが総体としては正しいと言われても、中々納得しない人が多いであろう。どう考えるか、一通り見てください。
superd options
x, # data
addrc=FALSE, # 変数名に r,c 付加
cex=0.8,
k=min(dim(x))-1, # 有効最大次元
axes=1:2,
para=TRUE, # ２分割で併記
tiaright=FALSE, # TIA を右に
addline=TRUE, # TIA に回帰線
docorr=FALSE, # corrplot 実行 doperm=TRUE, # permutation 実行
dname=NULL # data 名称
superd(norway)
superd(kitsuen)## Warning in chisq.test(x): Chi-squared approximation may be incorrect## Warning in chisq.test(t(x)): Chi-squared approximation may be incorrect
superd(Hab1)
superd(cosme)
superd(pastry)
superd(haireye)
superd(goodman.mental)
朝野氏の提示した PSVD
入手の遅れていた、朝野煕彦氏の論説は、雑誌そのものが個人購入できたので確認した。
朝野煕彦(2008) コレスポンデンス分析の空間表現ーPSVD の提案。Marketing Resercher, No.107, 43-53.
この提案 Pseudo-SVD とは、SVD の疑似ではなく、元の contingency table: tbl を DescTools::Untable し、それの indicator matrix: indm(Untable(as.matrix(tbl))) を pseudo-indicator matrix と呼ぶことに由来している。詳細は、原論説を参照されたい。 なお、関数 indm() は単独の script としてもよいが、pseudo-contingency table そのものを求める関数は FactoMineR::MCA を用いて
pct <- function(tbl) MCA(Untable(as.matrix(tbl)), graph=FALSE)$call$Xtot と定義すれば求めることができる。
pct(asa1)## a b c d y1 y2 y3
## 1 1 0 0 0 1 0 0
## 2 0 0 1 0 1 0 0
## 3 0 0 1 0 1 0 0
## 4 0 0 0 1 1 0 0
## 5 1 0 0 0 0 1 0
## 6 0 1 0 0 0 1 0
## 7 0 1 0 0 0 1 0
## 8 0 1 0 0 0 1 0
## 9 0 0 1 0 0 1 0
## 10 0 0 1 0 0 1 0
## 11 1 0 0 0 0 0 1
## 12 1 0 0 0 0 0 1
## 13 0 1 0 0 0 0 1
## 14 0 0 1 0 0 0 1
## 15 0 0 1 0 0 0 1
## 16 0 0 1 0 0 0 1
## 17 0 0 1 0 0 0 1
## 18 0 0 0 1 0 0 1
## 19 0 0 0 1 0 0 1
## 20 0 0 0 1 0 0 1但し、朝野氏は転置した t(pct(tbl)) を pseudo-contingency table と呼んでいる。 いずれにせよ、この pct の、Burt matrix を, weight 付きで求める。単なる burt 行列は t(pct(tbl)) %*% as.matrix(pct(tbl)) である。たが、これ以降の処理を見ていると、結局は Untable した２列の df を 多重対応分析していると思われる。つまり、朝野氏の手続きをまとめれば、次の関数となる。なお、２次元表示の座標を反転させるために、 crev というパラメータを入れた。
asano <- function(tbl, axes=1:2, crev=c(1,1)){
library(DescTools)
library(ca)
adat <- Untable(as.matrix(tbl))
colnames(adat) <- c("r","c")
mjres <- mjca(adat)
mcrev <- rep(1,ncol(mjres$colcoord))
mcrev[axes] <- crev
mjres$colcoord[,] <- mjres$colcoord %*% diag(mcrev)
plot(mjres, dim=axes)
return(invisible(list(unt=adat, mjres=mjres)))
} すると、スケーリングはともかく、結局は通常のCAと同じ結果となる。
asano(asa2, crev=c(-1,1))
CA(asa2)
## **Results of the Correspondence Analysis (CA)**
## The row variable has 6 categories; the column variable has 7 categories
## The chi square of independence between the two variables is equal to 4545.54 (p-value = 0 ).
## *The results are available in the following objects:
##
## name description
## 1 "$eig" "eigenvalues"
## 2 "$col" "results for the columns"
## 3 "$col$coord" "coord. for the columns"
## 4 "$col$cos2" "cos2 for the columns"
## 5 "$col$contrib" "contributions of the columns"
## 6 "$row" "results for the rows"
## 7 "$row$coord" "coord. for the rows"
## 8 "$row$cos2" "cos2 for the rows"
## 9 "$row$contrib" "contributions of the rows"
## 10 "$call" "summary called parameters"
## 11 "$call$marge.col" "weights of the columns"
## 12 "$call$marge.row" "weights of the rows"提唱された距離については、別途検討する。
朝野氏の実例
朝野氏の実例２件について、TIA の状況を見る。
superd(asa1)## Warning in chisq.test(x): Chi-squared approximation may be incorrect## Warning in chisq.test(t(x)): Chi-squared approximation may be incorrect
superd(asa2)
