Rで解析：応答適応型ランダム化を臨床試験の設計前にシミュレーションできる「grouprar」パッケージ

「grouprar」パッケージは、反応を見ながら割り当て比率を変える臨床試験（応答適応型ランダム化）を、実施前にシミュレーションで比べられるパッケージです。渡すのは、各群の真の成功率や平均・分散と症例数で、割り当ての偏り、失敗率、検出力が返ってきます。進行中の試験なら、これまでの割り当てと反応のベクトルから次の患者の割り当て確率を出せます。反応が欠損したり遅れて届いたりする場合や、患者がまとまって登録される場合も扱えます。試験計画を立てる人が、候補のデザインを事前に比べる道具として使えるのではないかと考えます。

パッケージバージョンは0.2.0。Windows 11 x64 (build 26200)のR version 4.6.1で確認しています。

＜おすすめのRに関する書籍です＞

Rによるやさしいテキストマイニング: 機械学習編

テキストマイニングに機械学習を組み合わせる方法を扱う一冊。分類や特徴抽出を実践的に学べる。

www.amazon.co.jp

パッケージのインストール

下記コマンドを実行してください。

# パッケージのインストール
install.packages("grouprar")

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

コマンド例

詳細はコメント、パッケージのヘルプを確認してください。例は、京都市伏見区の診療所で宇治茶カテキン飲料（A）、丹波黒大豆茶（B）、麦茶（C）を比べる架空の試験です。

DBCDは「Doubly-adaptive Biased Coin Design（二重適応型バイアスコイン設計）」の略で、それまでの反応から目標の割り当て比率を推定し直し、実際の割り当てをその目標へ寄せていく方法です。ERADEは「Efficient Randomized Adaptive DEsign」の略で、目標への寄せ方をDBCDより強くした方法です。目標の比率には、失敗数を減らす”RSIHR”、検出力を重視する”Neyman”、壺モデルに合わせた”RPW”などを選べます。

シミュレーションの結果は、各群の割り当て割合（allocation）とその標準偏差、失敗率または反応の平均、検出力（power）の順に表示されます。

進行中の試験で次の割り当て確率を計算：nextAllocコマンド

# これまでに登録した24人の割り当て（1:宇治茶 2:黒大豆茶 3:麦茶）
alloc <- c(1, 2, 3, 1, 2, 3, 3, 1, 2, 1, 3, 2, 1, 2, 3, 3, 1, 2, 3, 1, 2, 1, 3, 2)

# 疲労感の改善（1:改善 0:変化なし）。最後の4人は判定待ちなのでNA
improve <- c(1, 0, 1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 0, 1, 1, 1, 1, 0, 1, 1, NA, NA, NA, NA)

# 次の1人の割り当て確率をRSIHRの目標で計算
set.seed(10)
nextAlloc(alloc, improve, k = 3, target.alloc = "RSIHR")
$prob
treatment A treatment B treatment C
  0.4364095   0.1271810   0.4364095

$target
treatment A treatment B treatment C
  0.3755172   0.2489655   0.3755172

$estimate
       pA        pB        pC
0.8125000 0.3571429 0.8125000

$assignment
[1] 3

# 次に来院する4人の割り当てをERADEで抽選
nextAlloc(alloc, improve,
  k = 3, target.alloc = "RSIHR",
  allocation = "ERADE", size = 4
)$assignment
[1] 1 1 3 1

# 睡眠時間の延び（分）を連続の反応として扱う
sleep_min <- c(
  22, 8, 15, 30, 5, 18, 12, 26, 10, 19, 14, 3,
  24, 9, 20, 16, 28, 7, 13, 21, NA, NA, NA, NA
)

# Neymanの目標で次の1人の割り当て確率を計算
nextAlloc(alloc, sleep_min,
  k = 3, response = "continuous",
  target.alloc = "Neyman"
)$prob
treatment A treatment B treatment C
  0.6091584   0.1685191   0.2223225

Bai・Hu・Shenの壺モデルのシミュレーション：Bai.Hu.Shen.Urnコマンド

# 3群の改善率を想定し、150人の試験を200回シミュレーション
bhs <- Bai.Hu.Shen.Urn(
  k = 3, p = c(0.55, 0.70, 0.60), ssn = 150,
  nsim = 200, seed = 1
)

# 結果を表示
bhs
Bai, Hu and Shen's Urn
Sample size: 150   Simulations: 200

               p allocation    sd
treatment A 0.55      0.280 0.094
treatment B 0.70      0.411 0.117
treatment C 0.60      0.309 0.096

      failure rate sd of failure rate              power
             0.372              0.043              0.340

出生死滅型の壺モデルのシミュレーション：BirthDeathUrnコマンド

# 壺に各群の玉を2個ずつ入れて始める
bd <- BirthDeathUrn(
  k = 3, p = c(0.55, 0.70, 0.60), ssn = 150,
  Y0 = c(2, 2, 2), nsim = 200, seed = 1
)

# 結果を表示
bd
Birth and Death Urn
Sample size: 150   Simulations: 200

               p allocation    sd
treatment A 0.55      0.226 0.128
treatment B 0.70      0.489 0.173
treatment C 0.60      0.285 0.135

      failure rate sd of failure rate              power
             0.357              0.047              0.320

＜おすすめのRに関する書籍です＞

みんなのR 第2版 | Jared P. Lander, 高柳慎一, 津田真樹, 牧山幸史, 松村杏子, 簑田高志

Rの基礎からモダンな分析までを一冊で。第2版ではdplyr・purrr・tidyverse・caret・R Markdown・Shinyの章が加わった。

二値の反応でDBCDをシミュレーション：DBCD_Binコマンド

# 宇治茶と麦茶の2群。判定の1割が欠損すると想定
uji <- DBCD_Bin(
  p = c(0.55, 0.72), k = 2, ssn = 160, target.alloc = "RSIHR",
  nsim = 200, mRate = 0.1, seed = 1
)

# 結果を表示
uji
Hu and Zhang's DBCD
Sample size: 160   Simulations: 200

               p allocation    sd
treatment A 0.55      0.466 0.024
treatment B 0.72      0.534 0.024

      failure rate sd of failure rate              power
             0.361              0.036              0.550

# 帰無仮説の下での第1種の過誤を確認
uji_null <- DBCD_Bin(
  p = c(0.55, 0.72), k = 2, ssn = 160, target.alloc = "RSIHR",
  nsim = 200, typeI = TRUE, seed = 1
)
uji_null[["type I error"]]
[1] 0.025

# 症例数の半分で中間解析を1回入れる
uji_sq <- DBCD_Bin(
  p = c(0.55, 0.72), k = 2, ssn = 160, target.alloc = "RSIHR",
  nsim = 200, monitor = sqMonitor(c(0.5)), seed = 1
)

# 各解析時点で有効中止となった割合
uji_sq[["stopping probability"]]
look 1 look 2
 0.105  0.895

# 症例数の期待値
uji_sq[["expected sample size"]]
[1] 151.6

連続の反応でDBCDをシミュレーション：DBCD_Contコマンド

# 収縮期血圧の低下量（mmHg）。宇治茶は平均8・分散25、麦茶は平均5・分散16
bp <- DBCD_Cont(
  theta = c(8, 25, 5, 16), k = 2, ssn = 120,
  nsim = 200, seed = 2
)

# 結果を表示
bp
Hu and Zhang's DBCD (Gaussian Response)
Sample size: 120   Simulations: 200

            mu sigma2 allocation    sd
treatment A  8     25      0.559 0.042
treatment B  5     16      0.441 0.042

      mean response sd of mean response               power
              6.720               0.448               0.975

# 黒大豆茶（平均6・分散9）を加えた3群をZRの目標とERADEで
bp3 <- DBCD_Cont(
  n0 = 30, theta = c(8, 25, 5, 16, 6, 9), k = 3, ssn = 150,
  target.alloc = "ZR", lower.bound = 0.2,
  allocation = "ERADE", nsim = 200, seed = 2
)

# 結果の要約を表示
summary(bp3)
ERADE (Gaussian Response)
Sample size: 150   Simulations: 200

            mu sigma2 allocation    sd
treatment A  8     25      0.404 0.100
treatment B  5     16      0.373 0.106
treatment C  6      9      0.223 0.038

      mean response sd of mean response               power
              6.429               0.501               0.885

＜おすすめのRに関する書籍です＞

Rで楽しむベイズ統計入門［しくみから理解するベイズ推定の基礎］ Data Science Library

ベイズ推定のしくみを、仕組みの理解から丁寧に追う入門書。

www.amazon.co.jp

一般化ドロップ・ザ・ルーザー則のシミュレーション：GDLRuleコマンド

# 移入玉を引いたら各群の玉を1個ずつ足す
gdl <- GDLRule(
  k = 3, p = c(0.55, 0.70, 0.60), ssn = 150, aK = c(1, 1, 1),
  nsim = 200, seed = 1
)

# 結果を表示
gdl
Generalized Drop-the-loser Rule
Sample size: 150   Simulations: 200

               p allocation    sd
treatment A 0.55      0.288 0.038
treatment B 0.70      0.403 0.045
treatment C 0.60      0.310 0.043

      failure rate sd of failure rate              power
             0.374              0.039              0.335

二値の反応でグループ単位のDBCDをシミュレーション：Group.DBCD_Binコマンド

# 平均4人ずつまとまって来院し、判定の15%が欠損すると想定
gb <- Group.DBCD_Bin(
  p = c(0.55, 0.72), k = 2, gsize.param = 4, ssn = 160,
  target.alloc = "RSIHR", nsim = 200, mRate = 0.15, seed = 3
)

# 結果を表示
gb
Group DBCD with Binary Response (Random Missing)
Sample size: Total 160, Effective Size 136   Simulations: 200

               p allocation    sd
treatment A 0.55      0.467 0.024
treatment B 0.72      0.533 0.024

      failure rate sd of failure rate              power
             0.354              0.041              0.460

＜おすすめのRに関する書籍です＞

Rではじめるデータサイエンス 第2版 | Hadley Wickham, Mine Çetinkaya-Rundel, Garrett Grolemund, 大橋 真也

tidyverseの事実上の公式教科書。データの取り込み・整形・可視化・モデリング・レポート化までを一本の流れとして学べる。

連続の反応でグループ単位のDBCDをシミュレーション：Group.DBCD_Contコマンド

# 血圧の低下量を、平均4人ずつのまとまりで割り当てる場合
gc <- Group.DBCD_Cont(
  theta = c(8, 25, 5, 16), k = 2, gsize.param = 4,
  ssn = 120, nsim = 200, seed = 4
)

# 結果を表示
gc
Group DBCD with Continuous Response (No Missing)
Sample size: Total 120, Effective Size 120   Simulations: 200

            mu sigma2 allocation    sd
treatment A  8     25      0.556 0.045
treatment B  5     16      0.444 0.045

      mean response sd of mean response               power
              6.637               0.421               0.955

この記事が誰かの役に立ちますように。