R Mantel-Haenszelの勉強

Program R

「Rで学ぶデータサイエンス  1カテゴリカルデータ解析」を参照して作成

UCBAdmissions

GAdata <- apply(UCBAdmissions, c(2, 1), sum)

mosaicplot(GAdata)

chisq.test(GAdata)

 

#データを次々と読み込み、異なる変数名を付与する

j = 1

for (i in 0:5) {

      Dept <- matrix(c(UCBAdmissions[4*i+1], UCBAdmissions[4*i+2], UCBAdmissions[4*i+3], UCBAdmissions[4*i+4]),nrow=2,

              dimnames = list(Admin = c("Admitted", "Rejected"),Gender = c("Male", "Female")))

      assign(paste("Dept",j,sep=""),Dept)

      j = j + 1

}

 

par(mfrow=c(2,3))

mosaicplot(Dept1)

mosaicplot(Dept2)

mosaicplot(Dept3)

mosaicplot(Dept4)

mosaicplot(Dept5)

mosaicplot(Dept6)

 

mantelhaen.test(UCBAdmissions)

mantelhaen.test(UCBAdmissions, exact = TRUE)   #☆

#Mantel and Haenszelは、層別された2×2表のデータに対して、cbind(c(ak, bk), c(ck,dk))を仮定、それぞれの層で2つのcategorical変数の間に関連が無い、すなわちオッズ比が1であることを帰無仮説としている。帰無仮説のもとでは、akの平均μkと分散vkは、μk＝r1k*t1k/nk, vk = r1k*r2k*t1k*t2k*/nk^2 * (nk - 1)　※rik, tikはそれぞれk番目のtableのi行目、i列目の合計。ただし、ここではそれぞれの層で周辺の度数を固定した超幾何分布モデルを用いている。このことを利用して、すべての層のオッズ比が1であると仮定した場合の平均的な結果と実際に観測されたデータとのずれを、統計量X^2 ＝ { | ∑ak - ∑μk | -0.5}^2/∑vk を使って評価する。ここで、0.5を引いている部分は、独立性のχ2乗検定のときと同じように、連続性の補正を行っているため(?)。すべての層のオズ比が1であるという仮説が正しいとすると、χ2の分布は近似的に自由度1のχ2乗分布となるので、それほど大きな値をとることは無い。上記の方法は、∑akの分布が近似的に正規分布であることを仮定して構成している。これとは別の方法として、各層のオッズ比が1である場合の∑akの分布を、各層のデータに対して超幾何分布モデルを仮定して計算することも出来る。このアイデアは二次元の場合のフィッシャーの直接確率法と同様。☆の式に対応。