臨床予測モデルの開発に必要なサンプルサイズの算出（Calculating the sample size required for developing a clinical prediction model）

概要

ヘルスケア領域のデータ分析に携わっていると、よく質問されるのが

「どれくらいサンプルサイズが必要なの？」

という質問です。統計的仮説検定であれば、

「有意水準、検出力…等を指定して、必要サンプル数を求めることが可能だと思います」

とそれらしく回答しているのですが、予測モデルを構築する場合は、

「なるべく多くのサンプル集める。特徴量の次元数の10倍のデータが必要みたいなことは、経験則として提示されていますね…」

というような感じで、もごもごと回答しています。

そこで、google検索してみたところ「Calculating the sample size required for developing a clinical prediction model」という論文が存在していたようなので、流し読みしてみました。

www.bmj.com

この論文では、binary outcome（二値の目的変数）に対しての、logistic regression（ロジスティック回帰モデル）を用いた場合に限り、確率的な性質を利用してのサンプルサイズの見積もり手法が提案されていました。しかも、Rのpmsampsizeパッケージも紹介されていました。

そこで、今回の記事では、背景や細かい理論は気にせず、個人的なメモとしてRのpmsampsizeパッケージの使い方を書いています。

pmsampsizeの使い方

まずは必要なlibraryを読み込みます。最後にちょっとしたデータ加工と可視化をするので、tidyverseも読みます。

library(pmsampsize)
library(tidyverse)

pmsampsizeのpmsampsize関数を用いることで、サンプルサイズを見積もることが出来ます。

例えば、目的変数が0/1の2値変数、ロジスティック回帰モデルで、説明変数が5個、1の出現割合が0.3、AUC=0.70とすると、以下のような結果が出力され、サンプルサイズ412と見積もることができました。

> pmsampsize(type = "b", cstatistic = 0.70, parameters = 5, prevalence = 0.3)

Given input C-statistic = 0.7  & prevalence = 0.3  
Cox-Snell R-sq = 0.1029  
 
NB: Assuming 0.05 acceptable difference in apparent & adjusted R-squared 
NB: Assuming 0.05 margin of error in estimation of intercept 
NB: Events per Predictor Parameter (EPP) assumes prevalence = 0.3  
 
           Samp_size Shrinkage Parameter    Rsq Max_Rsq   EPP
Criteria 1       412     0.900         5 0.1029    0.71 24.72
Criteria 2       131     0.743         5 0.1029    0.71  7.86
Criteria 3       323     0.900         5 0.1029    0.71 19.38
Final            412     0.900         5 0.1029    0.71 24.72
 
 Minimum sample size required for new model development based on user inputs = 412, 
 with 124 events (assuming an outcome prevalence = 0.3) and an EPP = 24.72

また、1の出現割合を、0.1, 0.2, ..., 0.4と動かしたり、説明変数の個数を1,2,...,20と動かしたりすることで、様々なパターンのサンプルサイズを見積もることが出来ます。

df_pm <- tibble()
p_s <- c(0.1, 0.2, 0.3, 0.4)

for (p in p_s) {
  for (para in seq(1, 20, 1)) {
    pm_tmp <- pmsampsize(type = "b", cstatistic = 0.70, parameters = para, prevalence = p)
    df_tmp <- tibble(１の出現割合=p, パラメータ数=para, サンプル数=pm_tmp$sample_size)
    df_pm <- bind_rows(df_pm, df_tmp)
  }
}

df_pm$１の出現割合 <- factor(as.character(df_pm$１の出現割合))

g <- ggplot(df_pm ,aes(パラメータ数, サンプル数, colour=１の出現割合, group=１の出現割合)) + 
  geom_line(size=1) + geom_point() + ggtitle("必要なサンプル")

plot(g)

f:id:mikutaifuku:20211015232738p:plain

さいごに

clinical prediction model（臨床予測モデル）、しかも0/1の2値のロジスティック回帰モデルに限定しての話で、おそらく一般的な機械学習だとこの限りではないと思います。
また、論文内では、train/testのようなデータの分割（ホールドアウト検証）は推奨されておらず、ブートストラップ法による検証が推奨されており、このあたりの感覚は少し違和感があります。
雰囲気でブログを書いただけなので、おかしな事をやっていたり、間違っていたらこっそりご連絡ください。