【R】モデルを用いたオッズ比の求め方

今回はRで、既存のデータセットを使用して、オッズ比計算をデモンストレーションしてみます。使用するのはEpiパッケージのdietデータセットです。

このデータは、冠状動脈性心疾患イベントと食事の関係性を評価するためのデータセットです。

まずはパッケージをインポートします。

library(tidyverse)
library(Epi)
library(broom)

dietデータセットをインポートします。

#オリジナルデータはClayton and Hills (1993).
data(diet)
head(diet)

#id: 被験者識別子、数値ベクトル。
#doe: フォローアップ研究への参加日、日付変数。
#dox：フォローアップ研究からの離脱日、日付変数。
#dob：生年月日、日付変数。
#y: 危険年数、数値ベクトル。(doxとdoeの差)
#fail: 終了時の状態、数値ベクトル（コード1,3,13はCHDイベントを表す）。
#job: 職業。レベル Driver Conductor Bank worker を持つ因子。
#month：食事調査の月、数値ベクトル。
#energy：総エネルギー摂取量（kCal/日/100）、数値ベクトル。
#height：（cm）、数値ベクトル。
#weight：(kg)、数値ベクトル。
#fat：脂肪摂取量（10g/日）、数値ベクトル。
#fibre：食物繊維の摂取量（10g/日）、数値ベクトル。
#energy.grp：1日のエネルギー摂取量、レベル<=2750 kCal, >2750 kCalの因子。
#chd：冠状動脈性心疾患イベント、数値ベクトル（1=CHDイベント、0=イベントなし）。

平均値を見た後に、全体を可視化してみます。

summary(diet)
#全体像を眺める
di <- select(diet, y, energy, height, weight, fat, fibre, chd)
plot(di)

energy~fibreは相関がみられることは容易に想像できますね。chdについては、傾向は見づらいです。

では、モデルを作成します。変数chdは0/1なので、二項分布を仮定したモデルを組みます。モデルには職業を表す変数job（driver, conductor, bank worker）も含めることにしました。（※本データはコホート研究なので、オッズ比ではなく相対リスクが求められる気もしますが、今回はORとします）

#ロジット-二項モデル
#説明変数に脂肪、食物繊維、摂取エネルギー、体重、調査月、職業
fit <- glm(chd ~ fat + fibre + energy + weight + month + job + y,
           data = diet,
           family = binomial(link="logit"))

#結果
summary(fit)

Call:
glm(formula = chd ~ fat + fibre + energy + weight + month + job + 
    y, family = binomial(link = "logit"), data = diet)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.0752  -0.4211  -0.2567  -0.1499   2.9383  

Coefficients:
               Estimate Std. Error z value Pr(>|z|)    
(Intercept)     2.41407    1.78929   1.349  0.17728    
fat            -0.25932    0.15530  -1.670  0.09496 .  
fibre          -0.99924    0.57815  -1.728  0.08393 .  
energy          0.09327    0.09391   0.993  0.32064    
weight          0.01472    0.02108   0.698  0.48509    
month          -0.05957    0.05790  -1.029  0.30357    
jobConductor    1.03778    0.54414   1.907  0.05650 .  
jobBank worker  1.55917    0.54841   2.843  0.00447 ** 
y              -0.32176    0.04614  -6.973  3.1e-12 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 258.88  on 328  degrees of freedom
Residual deviance: 174.04  on 320  degrees of freedom
  ( 8 個の観測値が欠損のため削除されました )
AIC: 192.04

Number of Fisher Scoring iterations: 6

p値で有意なのはjob（Bank worker）とy（危険年数≒観察期間）であることが分かりました。偏回帰係数estimateでは読み取りづらいため、指数変換してオッズ比を計算します。

#偏回帰係数estimateを指数変換して調整オッズ比を計算する
coef(summary(fit)) %>% as.data.frame() %>%
  mutate(OR = exp(Estimate),
         conf.low = exp(Estimate - 1.96*`Std. Error`),
         conf.high = exp(Estimate + 1.96*`Std. Error`)) %>%
  dplyr::select(OR, conf.low, conf.high)

このままでもよいのですが、あとで図にしたいのでdata.frame形式として出力するためbloomパッケージを使います。

#bloomパッケージのtidyを使う
#,c(1,2,6,7)によりORと95%信頼区間のみ表示する
OR <- tidy(fit, exponentiate = T, conf.int = T)[,c(1,2,6,7)]
OR

微妙に結果が異なるのは、信頼区間の計算に用いた数値(1.96)の桁数の違いと考えられます。では、この結果をフォレストプロットとして図示します。

#フォレストプロットで表示
forest <- OR %>%
  #intercept情報を削除
  filter(term != "(Intercept)") %>%
  #y軸の順番を指定する
  mutate(term = fct_relevel(
    term,
    "fat","fibre","energy","weight","month","jobConductor","jobBank worker","y")) %>%
  #グラフの構造説明
  ggplot(aes(x=estimate, y=term))+
  #estimateを点表示
  geom_point()+
  #信頼区間をエラーバーとして表示
  geom_errorbar(aes(xmin=conf.low, xmax=conf.high))+
  #RR=1を示す線を表示
  geom_vline(xintercept = 1, linetype="dashed")+
  #軸にラベル付け
  labs(x = "OR", y = "factors")
forest

結果が出ました。冠状動脈性心疾患と脂肪、食物繊維、摂取エネルギー、体重、調査月に関連はみられないことがわかりました。また、職業がDriverの人に対して、ConductorとBank workerは冠状動脈性心疾患のリスクが高いと言えることも見て取れます。

今回は、既存データを用いてオッズ比の計算を行ってみました。参考になれば幸いです。

Hazy Ideas

日々の勉強の気づきを書き出しています

【R】モデルを用いたオッズ比の求め方