クロス集計表のセルごとの有意確率を出してみよう
他人への質問を眺めていても学ぶことはあります。昨日は「クロス集計の独立性検定を示した表で、セルごとにアスタリスクがあるのは何故か?」と言う質問がされていましたが、私も理由が思いつきませんでした。以下のような頻度表、時折、見てきた気がするのですが、どうやって個々のセルの有意性を出しているのでしょうか?
JPN | USA | CHN | |
---|---|---|---|
MEAT | 10 | 13 | 6 |
FISH | 9 | 18** | 3** |
BEAN | 20 | 15** | 31** |
Note) ** and * are statistically significant at the 1% and 5% level, respectively.
独立性検定(全体)
まずはχ二乗検定をしてみましょう*1。
頻度が入った行列mを作ってchisq.test(m)で終わるのですが、今回は教科書的にチマチマと計算してみます。ただし、1%、5%棄却粋のχ二乗値を出して比較することはせずに、P値は計算してもらいます。
# 頻度の行列をつくる
m <- matrix(c(10,9,20,13,18,15,6,3,31), 3)# 見やすいように行名と列名をつけておく
colnames(m) <- c("JPN", "USA", "CHN")
rownames(m) <- c("MEAT", "FISH", "BEAN")# 列ごとの合計
apply(m, 2, sum)# 行ごとの合計
apply(m, 1, sum)# 全体の合計
sum(m)# 自由度を計算する
degf <- (ncol(m) - 1)*(nrow(m) - 1)# 行(MEAT, FISH, BEAN)ごとにの期待確率
exp_p <- apply(m, 1, sum)/sum(m)# 列(JPN, USA, CHN)ごとに期待確率exp_pをかけて、期待値の行列をつくる
exp_m <- sapply(apply(m, 2, sum), function(x){
x * exp_p
})# 観測値の行列と期待値の行列のセルごとの差分の二乗を合計して期待値で割り、χ二乗値を求める
chisq_v <- sum((m - exp_m)^2/exp_m)# P値を求める
pchisq(chisq_v, degf, lower.tail=FALSE)
当然ですが、0.0008233372と、一つしかP値は出ません。
セルごとに検定
残差の二乗値を期待値で割ったものの合計がχ二乗分布に従うことから検定をかけたわけですが、今度は個々のセルの調整済み*2残差が正規分布に従うことから検定をかけます。
これも2*pnorm(abs(chisq.test(m)$stdres), lower.tail=FALSE)で終わるのですが、チマチマとやっていきましょう。
# (1-列合計/合計)(1-行合計/合計)の修正用係数の行列
adj_m <- matrix(rep(1-apply(m, 2, sum)/sum(m), ncol(m)) * rep(1-apply(m, 1, sum)/sum(m), each=nrow(m)), nrow(m), byrow=TRUE)# 標準化残差
stdres <- (m - exp_m)/sqrt(c(adj_m) * c(exp_m))# セルごとのP値をまとめた行列をつくる
m_p_value <- 2*pnorm(abs(stdres), lower.tail=FALSE)
今回の結果は以下のようになりました。
JPN | USA | CHN | |
---|---|---|---|
MEAT | 0.6633 | 0.3064 | 0.1362 |
FISH | 0.8707 | 0.0025 | 0.0030 |
BEAN | 0.8189 | 0.0006 | 0.0001 |
多重比較のためのP値の調整
期待値1%の事象でも、100回繰り返すと64%の確率で1度は生じることになります。この理屈でクロス集計表が細かく切られている場合、偽陽性が出やすくなります。多重比較のために、P値を調整しましょう。ここはライブラリにお任せです。
m_adj_p_value <- matrix(p.adjust(m_p_value), nrow(m_p_value))
dimnames(m_adj_p_value) <- dimnames(m) # 見やすいように行名と列名をコピーしておく
p.adjustでとれる調整法は多彩ですが、デフォルトはHolm法になっています。
JPN | USA | CHN | |
---|---|---|---|
MEAT | 1 | 1.0000 | 0.6812 |
FISH | 1 | 0.0176 | 0.0183 |
BEAN | 1 | 0.0045 | 0.0013 |
今回は解釈を左右するような変化はないですが、P値が大きくなっていることに注意してください。FISH行の1%有意が5%有意になりました。なお、この手のP値の調整は偽陰性が大きくなる傾向がある*3そうです。
参考文献
本稿の作成において「クロス集計表の有意差検定 | 象と散歩」と「Rでクロス集計表の残差分析 - bob3’s blog」を参考にしました。
*1: サンプルサイズが小さく、セルの度数が一つでも5未満の場合は、多項分布を前提にしたFisherの正確検定を使うべきと言う話もありますが、ここでは無視します。
*2: 自由度を調整する作業になると思いますが、Agresti (2018) などで確認してください。Rのヘルプの雰囲気ではたぶん、理由が分かるはずです。たぶん(´・ω・`)ショボーン