sacoche

公衆衛生学 / 疫学 / R. 知識をサッと持ち出せると良いよね

稀な疾病ってどの程度"稀"なんですか？

疫学・生物統計学

疾病の有無やイベント発生の有無など2値を取るアウトカムに関する分析の中で，リスク比 (Risk Ratio; RR)とオッズ比 (Odds Ratio; OR) はとても頻繁に登場する効果・関連の指標です．

リスク比・オッズ比を勉強する際に，必ず聞くフレーズがあります．「稀な疾病・イベントの場合，オッズ比はリスク比に近似できる」です．

しかしながら，稀ってどの程度？とか，稀でないときどうする？はそこまで有名ではないかなと．

そこで，このトピックについて，2回の記事で簡単に（？）まとめていきましょう．今回の記事では，「稀とは？」を説明します．

 

目次

 

稀とはイベント発生割合が10％以内

ORおよびRRが1以上の状況について，考えます．

「稀」の判断を行う場合，一般に10％が基準となります．10％を上回るイベント発生割合の場合はORがRRを上回ることが知られています．

イベント発生のあり/なしなど2値アウトカムの分析では，ロジスティック回帰が頻繁に使用されますね．イベントの発生割合が10％以下であれば，ロジスティック回帰によるオッズ比を効果の指標として用いることは適切であると考えられます．一方，10%を上回る場合は要注意です．OR > RRとなり，効果を誇張した (exaggerated) 結果になるかもしれません（稀でないイベントの回帰分析は次回）．

 

稀でない場合，ORとRRにどの程度の差が生まれるのか

稀とは10％が基準になる，という話は提示しましたが，「具体的にどの程度困るのか」にはまだ触れられていません．そこで，様々なイベント発生割合でデータを作り，ORとRRの関係を図示してみましょう．

データ作り

（結果だけ見たい人は データの可視化 まで飛ばしてください）

縦断的なコホート研究を通じて，データが得られたとします．この解析では交絡は存在しないと仮定し，シンプルな2×2表の分析をします．

観察集団は曝露群1000人・非曝露群1000人，合計2000人とします．また，集団全体のイベント発生割合 p_{overall} は次の5通りを用意します．

  1. 5%
  2. 10%
  3. 20%
  4. 30%
  5. 50%

曝露群・非曝露群のイベントありがそれぞれ a 人，b 人とすると，2×2表は次の通りとなります．

  イベントあり(人) イベントなし(人) 合計(人)
曝露あり a 1000-a 1000
曝露なし b 1000-b 1000

 

さて，イベント発生割合 p_{overall} は次の式を満たします．

p_{overall} = \frac{a+b}{2000}

したがって，b = 2000\times p_{overall} - a と変形できますね．p_{overall}は自分たちで5%から50%と決めているので， aを与えれば，2×2表のすべての値が定まります．

 

RRとORは次の通りです．

\mathrm{RR} = \frac{a/1000}{b/1000} =\frac{a}{b} = \frac{a}{2000p_{overall} - a}

\mathrm{OR} = \frac{a\times (1000-b)}{c\times (1000-a)}=\frac{a(1000+a-2000p_{overall})}{(1000-a)(2000p_{overall}-a)}

以上のロジックの下，様々なイベント発生割合を設定した上で 1\le a \le 999 で a を動かしてデータを作成し，ORとRRの関係を可視化します．

 

なお，Rコードは記事の最後にまとめておきます．

 

データの可視化

グラフ化すると次のとおりです．

f:id:saco_che:20220109163913p:plain

様々なイベント発生割合における
リスク比とオッズ比の関係

各イベント発生割合ごとに，RRとORの関係を描きました．グラフ中のグレーで引かれている y=x の直線に近いほど，RRとORの値の乖離が小さいことを示しています．

例えば， p_{overall}が5%の曲線を見ると， y=x とあまり離れていません．したがって， \mathrm{OR} \simeq \mathrm{RR} と言っても大きな問題はなさそうです．一方， p_{overall}が50%の曲線を見ると， y=x から大きく離れていることがわかります．これほどのイベント発生割合になると，OR > RR となることがはっきりわかりますね．

 

具体的な数値を確認しましょう．例えばRR=1.5とRR=3.0の場合，それぞれのイベント発生割合におけるORは次のとおりです．

イベント発生割合 RR=1.50の場合のOR RR=3.00の場合のOR
5% 1.53 3.16
10% 1.57 3.35
20% 1.66 3.86
30% 1.78 4.64
50% 2.25 9.00

30%以上になると，OR≒RRと言ってしまうとまずいレベルに乖離していますね．

 

「ORはORとして解釈すればいいじゃないか！乖離して何が問題なんだ！」という意見もあるかと思います．しかし，もともとORはRRの代替指標として扱われている以上，近似できない場合はRRを捉えるべきだと思います．

 

今回はこれでおしまい

まとめ

  • 稀な疾病の場合OR≒RRと考えることができる
  • 「稀な」とは10%が一般的な基準であり，イベント発生割合が10%を超過する場合はOR > RRとなる (ORおよびRRが1以上の場合)
  • イベント発生割合が10%を超過する場合，ロジスティック回帰などORを算出する分析を適用する際には注意が必要であり，RRを導く解析の方が妥当だと考えられる．

そのうち，RRを導く解析として，修正ポアソン回帰 (modified poisson regression) と対数二項回帰 (log binomial regression) を簡単に紹介します．

 

論文読みたい人は...

この辺読めばいいと思います．

Schmidtらの論文の図は今回の図と同じはず．
Fuyamaさんの論文の話は聞いたことあるんだけど，やっぱり面白いなあと．多分読者の関心には沿っていないと思いますが，興味ある方はぜひ．

  • Schmidt, C. O., & Kohlmann, T. (2008). When to use the odds ratio or the relative risk?. International journal of public health, 53(3), 165.
  • Fuyama, K., Hagiwara, Y. & Matsuyama, Y. A simulation study of regression approaches for estimating risk ratios in the presence of multiple confounders. Emerg Themes Epidemiol 18, 18 (2021). 

 

付録: Rコード

データ作成
a <- c(1:999)
df <- data.frame(a)

df <- df %>% 
  # set overall risk
  mutate(p1 = 0.05) %>% 
  mutate(p2 = 0.1) %>% 
  mutate(p3 = 0.2) %>% 
  mutate(p4 = 0.3) %>% 
  mutate(p5 = 0.5)

df2 <- df %>% 
  pivot_longer(cols = p1:p5,
               names_to = "p_overall",
               values_to = "p") %>% 
  mutate(RR = a/(2000*p - a)) %>% 
  mutate(OR = a*(1000+a-2000*p)/((1000-a)*(2000*p-a)) ) %>% 
  
  mutate(p_overall = case_when(
    p_overall == "p1" ~ "5%",
    p_overall == "p2" ~ "10%",
    p_overall == "p3" ~ "20%",
    p_overall == "p4" ~ "30%",
    p_overall == "p5" ~ "50%"
  ),
  p_overall = factor(p_overall, levels = c("5%", "10%", "20%","30%", "50%")))

現在の状態では，RRやORが負の値となっていたり，無限大（0割りの結果）となる場合も含まれうるので，除外して整理します．

今回は，次の範囲で描写しますね．

  • 1\le \mathrm{RR} \le 10
  • 1\le \mathrm{OR} \le 10
df2 <- df2 %>% 
  filter(RR >= 1 & OR >=1 & RR <= 10 & OR <= 10)

以上でデータ作りは終了です．

 

グラフの描画
plot <- ggplot(data = df2, aes(x = RR, y = OR, group = p_overall, colour = p_overall)) + 
  geom_line(size = 1) + 
  geom_abline(slope = 1, intercept = 0, alpha = 0.8, colour = "darkgray") + 
  geom_point(x = 1, y = 1, colour = "black") +
  scale_x_continuous(limits = c(1, 10), breaks = c(1, 2, 4, 6, 8, 10)) +
  scale_y_continuous(limits = c(1, 10), breaks = c(1, 2, 4, 6, 8, 10)) +
  xlab("Risk Ratio") + 
  ylab("Odds Ratio") + 
  ggtitle("Risk Ratio vs Odds Ratio") +
  
  annotate("text", x = 3, y = 10, label = "p[overall] == 0.5", size = 3, parse = TRUE) + 
  annotate("text", x = 5.5, y = 10, label = "p[overall] == 0.3", size = 3, parse = TRUE) + 
  annotate("text", x = 6.7, y = 10, label = "p[overall] == 0.2", size = 3, parse = TRUE) +
  annotate("text", x = 8, y = 10, label = "p[overall] == 0.1", size = 3, parse = TRUE) + 
  annotate("text", x = 9.3, y = 10, label = "p[overall] == 0.05", size = 3, parse = TRUE) + 
 
  annotate("text", x = 1.8, y = 1, label = "RR = OR = 1", size = 3) + 

  theme_light()

plot

 

RR=1.5の場合のORの算出（RR=3.0でも同様）

df2 %>% 
  filter(round(RR, digits = 2) == 1.5)

 

おつかれさまでした〜