Rでオッズ比と調整済み残差を出す

2x2のクロス集計表の場合はオッズ比と調整済み標準化残差は似たような統計量になる。どちらが好まれるかは分野によって異なるので、適宜、好まれる方を使用するのがよいと思う。社会学ではそれほど見かけることはないが使用されるはずである。オッズ比はなぜか知らないがあまり好まれない傾向にあるように感じる。ちなみに医学ではオッズ比が好まれ、調整済み標準化残差を見かけたことがない。

オッズ比は95%信頼区間を併記するのが普通である。調整済み標準化残差は5%基準、1%基準、0.1%基準をアスタリスク等で表記する。95%信頼区間はどちらかが1を越えていないかを一つ一つ確認しなければならないので、表を見るのがめんどくさい。一方で、調整済み標準化残差のアスタリスク表記は、ぱっと見で分かるので便利だと思う。

有名な統計パッケージでは、SPSSでは調整済み標準化残差は計算できるが、オッズ比の計算ができない。stataはマクロを入れると計算できたように記憶している。とりあえず、計算が一番簡単で効率よくできるのはRであるのは間違いない。

下準備

クロス集計表のマトリックスを作成しておく。

例1 データマトリックス

d1 <- matrix(c(2,15,30,12), nrow=2)

データセットからも作成しておく。

例2 CPS1985のデータを使う

library(AER) #AERパッケージの読み込み
data(CPS1985) #CPS1985の読み込み
d2 <- table(CPS1985$gender, CPS1985$married)

性別と結婚ステータスを選んだ。2x2のクロス表が作りたいのでどちらも二値データの変数である。実際の分析としてはあまり意味がないかもしれない。

調整済み残差

調整済み標準化残差はRのデフォルト機能で出力できる。

d1res<- chisq.test(d1) #カイ二乗検定
d1res$stdres # 調整済み残差

結果は以下。

   [,1]      [,2]
[1,] -4.166099  4.166099
[2,]  4.166099 -4.166099

絶対値4.17である。

d2res <- chisq.test(d2) #カイ二乗検定
d2res$stdres # 調整済み残差

結果は以下。

    no       yes
male    0.259397 -0.259397
female -0.259397  0.259397

こちらは絶対値0.26である。

調整済み標準化残差は標準正規分布区間推定値と同様の値をとる。つまり、95%信頼区間は-1.96~1.96であり、残差ののように期待値からは外れていることを示す数値はこの95%信頼区間の外でああればよい。したがって、1.96より大きいと5%水準で有意な差があることになる。どうように、1%水準は2.58、0.1%水準は3.29、0.01%水準は3.89以上になる。

したがって、例1は4.17であったので0.01%水準で有意だったが、例2は1.96に満たないため、有意ではないということになる。

標準正規分布(余談)

ちなみに1.96というような数字は標準化しているために使用できる基準である。標準正規分布表からこの数字は読み取れる。

正規分布の上側確率
http://aoki2.si.gunma-u.ac.jp/CGI-BIN/tgxp.html

リンク先のでは1.96が例に上がっている。1.96のセルは0.0250であるで、これを2倍(両側にするため)すると0.05=5%水準となる。

オッズ比

オッズ比を出すには、パッケージのインストールが必要である。いくつか存在しているが、一番いろいろできるのはepitoolsパッケージである。

library(epitools) #epitoolsパッケージの読み込み

mid-p法(median-unbiased estimation)

いわゆるオッズ比というとこの方法のことである。

oddsratio.midp(d1) #midpをつけた形で
oddsratio(d2) #oddsratioのみで

単純にオッズ比といったときにはmid-p法のことなので、oddsratio()で構わない。 d1の方を出力すると以下のようになる。

$data
          Outcome
Predictor  Disease1 Disease2 Total
  Exposed1        2       30    32
  Exposed2       15       12    27
  Total          17       42    59

$measure
          odds ratio with 95% C.I.
Predictor    estimate      lower     upper
  Exposed1 1.00000000         NA        NA
  Exposed2 0.05953432 0.00773934 0.2558031

$p.value
          two-sided
Predictor    midp.exact fisher.exact   chi.square
  Exposed1           NA           NA           NA
  Exposed2 3.413834e-05  3.46431e-05 3.098556e-05

$correction
[1] FALSE

attr(,"method")
[1] "median-unbiased estimate & mid-p exact CI"

必要なのは$measureのところなのでそこだけ取り出して出力する。

oddsratio(d1)$measure #後ろに$measureをつける

そうするとオッズ比と95%区間のみが表示される。

          odds ratio with 95% C.I.
Predictor    estimate      lower     upper
  Exposed1 1.00000000         NA        NA
  Exposed2 0.05953432 0.00773934 0.2558031

オッズ比が0.059、95%信頼区間が[0.008, 0.256]である。 クロス表の関連を総当たりで確かめることをよくするが、そういう時に便利だと思う。

その他の方法(参考)

詳しくはこちらを参照。

oddsratio.fisher(d1) # Fisherの方法
oddsratio.wald(d1) # Waldの方法
oddsratio.small(d1) # 正規分布近似+小標本補正