R言語でオッズ比を計算する関数やその実行例を紹介します。
2つのグループにおける事象の起こりやすさの尺度であるオッズ比のR言語での計算例について見ていきます。
この記事で扱うRのコードは次のzipファイルに保存しています。
オッズ比について詳しく見たい方は次の記事を参照。
【統計学】標本オッズ比 オッズ比の計算 記述統計
2つのグループにおける事象の起こりやすさの尺度であるオッズ比について解説する。 オッズ比の定義や分割表が与えられたときのオッズの計算例を紹介する。 様々な分割表のオッズ比を計算し、オッズ比がどのように ...
続きを見る
オッズ比
オッズ比の概要とR言語でオッズ比を計算する関数をまとめました。
この記事で紹介するオッズ比の関数は下のオッズ比の概要で定義されるものとなります。
オッズ比の概要
オッズ比は2つの群におけるある事象の起こりやすさの尺度であり、2つの群のオッズの比として定義されます。以下にオッズ比の定義をまとめました。標本に関するオッズ比であることに注意してください。
標本オッズ比
\(x_1, x_2, \ldots, x_m\)をベルヌーイ分布\(Bernoulli(p)\)、\(y_1, y_2, \ldots , y_n\)をベルヌーイ分布\(Bernoulli(q)\)からの無作為標本とする。2つの母集団の標本比率をそれぞれ\(\hat{p} = m^{-1} \sum_{i = 1}^m x_i\)、\(\hat{q} = n^{-1} \sum_{i = 1}^n y_i\)とすると、標本オッズ比は次で定義される。
2つの群を\(A\)と\(B\)とし、それぞれの群の事象の発生に関するベルヌーイ分布からの標本を\(x_1, x_2, \ldots , x_m\)、\(y_1, y_2, \ldots, y_n\)とする。このとき、2×2分割表を用いることで事象の発生数をつぎのようにまとめることができます。
このとき、各群のオッズはそれぞれ\(\sum_{i = 1}^m x_i / (m - \sum_{i = 1}^m x_i )\)、\(\sum_{j = 1}^m y_j / (n - \sum_{j = 1}^m y_ij )\)であり、この比がオッズ比\(OR\)となります。2つの群の事象の発生数が同じであるときオッズ比は\(1\)となり、発生数が乖離するほど\(0\)または\(\infty\)になります。
関数odds.ratio
パッケージquestionrの関数odds.ratioを用いることでオッズ比を計算することができます。以下、関数とその引数です。
odds.ratio(x, ...)
# S3 method for glm odds.ratio(x, level = 0.95, ...)
# S3 method for multinom odds.ratio(x, level = 0.95, ...)
# S3 method for factor odds.ratio(x, fac, level = 0.95, ...)
# S3 method for table odds.ratio(x, level = 0.95, ...)
# S3 method for matrix odds.ratio(x, level = 0.95, ...)
# S3 method for numeric odds.ratio(x, y, level = 0.95, ...)
# S3 method for odds.ratio print(x, signif.stars = TRUE, ...)
x | オブジェクト。xに関するのオッズ比が計算される。テーブルや行列だけでなくglmなどのオブジェクトにも適用できる。 |
... | 他のメソッドに渡す引数。 |
level | オッズ比の信頼区間を計算する際の信頼水準。 |
fac | 第二群を表すfactor型のオブジェクト。 |
y | 第二群のnumeric型のオブジェクト。 |
signif.stars | logical型。TRUEのときp値はアスタリスクを用いて視覚的に表現される。 |
実行例
関数odds.ratioの実行例を紹介します。先ほどのオッズ比の概要で紹介した関数odds.ratioのためにパッケージquestiorn、グラフの作成のためにggplot2とgridExtraを用います。事前にインストールしておいてください。
1 2 3 4 5 6 | install.packages("questionr") install.packages("ggplot2") install.packages("gridExtra") library(questionr) library(ggplot2) library(gridExtra) |
次の2×2分割表を表す行列を例にオッズ比の計算を行っていきます。
1 2 3 | X <- matrix(c(22, 8, 11, 19), nrow = 2, dimnames = list(c("A", "B"), c("occurred", "not occurred"))) Y <- matrix(c(6, 24, 18, 12), nrow = 2, dimnames = list(c("A", "B"), c("occurred", "not occurred"))) Z <- matrix(rep(c(10, 5), each = 2), nrow = 2, dimnames = list(c("A", "B"), c("occurred", "not occurred"))) |
次のように行列X、Y、Zは、行にグループAとB、列にイベントが発生したかを表すクラスから成る2×2分割表となります。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 | > X occurred not occurred A 22 8 B 11 19 > > Y occurred not occurred A 6 24 B 18 12 > > Z occurred not occurred A 10 5 B 10 5 |
X、Y、Zのオッズ比を計算するには、次のように関数odds.ratioの引数にオッズ比を計算したいデータを渡します。odds.ratioを実行するとオッズ比に関するオブジェクトが返され、このオブジェクトに$ORを付けることでオッズ比を取得することができます。
1 2 3 4 5 6 7 | oddsRatio_X <- odds.ratio(X) oddsRatio_Y <- odds.ratio(Y) oddsRatio_Z <- odds.ratio(Z) OR_X <- oddsRatio_X$OR OR_Y <- oddsRatio_Y$OR OR_Z <- oddsRatio_Z$OR |
X、Y、Zのオッズ比はそれぞれ次の通りです。
1 2 3 4 5 6 7 8 | > OR_X [1] 4.61622 > > OR_Y [1] 0.1723106 > > OR_Z [1] 1 |
計算したオッズ比をcsvファイルや他のファイルに出力したい場合は、次のように関数write.csv(またはwrite.table)を用います。
1 2 3 4 5 | statsTable <- data.frame(odds = rbind(X[, "occurred"] / X[, "not occurred"], Y[, "occurred"] / Y[, "not occurred"], Z[, "occurred"] / Z[, "not occurred"]), OR = c(OR_X, OR_Y, OR_Z), row.names = c("X", "Y" ,"Z")) write.csv(statsTable, "オッズ比.csv") |
上記を実行することで作業ディレクトリに次の画像のようなcsvファイルを作成することができます。
余談ですが、パッケージggplot2などに実装されているヒートマップを用いることでオッズ比を視覚的に把握することができます。
次を実行すると、X、Y、Zの分割表の各群AとBのイベントの発生確率の推定値の大きさのヒートマップを描くことができます。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 | dataset_X_plot <- data.frame(XProb = c(X / rowSums(X)), row = as.factor(rev(rep(seq_len(nrow(X)), nrow(X)))), col = as.factor(rep(seq_len(nrow(X)), each = nrow(X)))) dataset_Y_plot <- data.frame(YProb = c(Y / rowSums(Y)), row = as.factor(rev(rep(seq_len(nrow(Y)), nrow(Y)))), col = as.factor(rep(seq_len(nrow(Y)), each = nrow(Y)))) dataset_Z_plot <- data.frame(ZProb = c(Z / rowSums(Z)), row = as.factor(rev(rep(seq_len(nrow(Z)), nrow(Z)))), col = as.factor(rep(seq_len(nrow(Z)), each = nrow(Z)))) tile_X <- ggplot(data = dataset_X_plot, aes(x = col, y = row, fill = XProb, label = round(XProb, 2))) + geom_tile() + scale_fill_gradient2(low = "white", high = "#832424", midpoint = 0, space = "Lab", guide = "colourbar", limits = c(0, 1)) + geom_text() + ggtitle(paste0("OR: ", round(OR_X, 2))) tile_Y <- ggplot(data = dataset_Y_plot, aes(x = col, y = row, fill = YProb, label = round(YProb, 2))) + geom_tile() + scale_fill_gradient2(low = "white", high = "#832424", midpoint = 0, space = "Lab", guide = "colourbar", limits = c(0, 1)) + geom_text() + ggtitle(paste0("OR: ", round(OR_Y, 2))) tile_Z <- ggplot(data = dataset_Z_plot, aes(x = col, y = row, fill = ZProb, label = round(ZProb, 2))) + geom_tile() + scale_fill_gradient2(low = "white", high = "#832424", midpoint = 0, space = "Lab", guide = "colourbar", limits = c(0, 1)) + geom_text() + ggtitle(paste0("OR: ", round(OR_Z, 2))) grid.arrange(tile_X, tile_Y, tile_Z) |
行方向に対して色の濃さが同じであるほどオッズ比が1に近づくのが理解できます。また、行方向について色の濃さ異なれば異なるほどオッズ比は無限大または0に向かうのが確認できます。
まとめ
R言語でオッズ比を計算する関数や例を解説しました。
パッケージquestionrの関数odds.ratioを用いることでオッズ比を計算することができます。
R言語ではオッズ比を計算する関数が標準で実装されていないため、オッズ比を計算するには自分で関数を実装するかこの記事で紹介したパッケージをインストールする必要があります。