Rで解析：入院や死亡をまとめて治療効果を比べるwin time法「wintime」パッケージ

「wintime」パッケージは、臨床試験の2群を、入院や死亡のように重さの違う複数のイベントをまとめて比べ、治療群が対照群より良い状態で過ごせた時間（win time）として治療効果を求めるパッケージです。渡すのは、イベントを軽い順に行、参加者を列に並べた発生日数の行列と、同じ形のイベント有無の行列、それに治療群か対照群かを示すベクトルです。年齢などの共変量を加えた解析や、ブートストラップによるp値の計算もできます。生存時間解析に慣れた人が、複合エンドポイントの結果を日数という分かりやすい単位で示したいときに役立つのではないかと考えます。

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

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

StanとRでベイズ統計モデリング Wonderful R 2 | 松浦健太郎, 石田基広

Stanによるベイズ統計モデリングの実践書。現実のデータ解析を念頭に、モデリングの考え方から丁寧に解説する。

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

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

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

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

コマンド例

詳細はコメント、パッケージのヘルプを確認してください。例は、京都府内の診療所に通う高齢者20人を、宇治茶を毎日飲む群（治療群）と飲まない群（対照群）に10人ずつ分け、1年間（365日）入院と死亡を追った架空の試験です。

win time法では、参加者の状態を「イベント無し」「入院した」「亡くなった」のように重さの順に並べ、各時点でどちらの群が軽い状態にいるかを比べます。その差を追跡期間で積み上げた日数が治療効果です。そのため、Timeとその有無を示すDeltaは、行を臨床的に軽いイベントから重いイベントの順に並べます。イベントが起きなかった人の日数には、観察を終えた日を入れます。

# 入院までの日数（入院しなかった人は観察を終えた日）
# 1〜10人目は対照群、11〜20人目は宇治茶を飲む群
hosp_day <- c(
  60, 120, 365, 45, 200, 90, 365, 30, 180, 365,
  210, 365, 95, 365, 300, 365, 150, 365, 365, 240
)

# 死亡までの日数（亡くならなかった人は観察を終えた日）
death_day <- c(
  250, 365, 365, 160, 365, 365, 365, 120, 365, 365,
  365, 365, 365, 365, 365, 365, 330, 365, 365, 365
)

# 入院したら1、しなかったら0
hosp_event <- c(
  1, 1, 0, 1, 1, 1, 0, 1, 1, 0,
  1, 0, 1, 0, 1, 0, 1, 0, 0, 1
)

# 亡くなったら1、亡くならなかったら0
death_event <- c(
  1, 0, 0, 1, 0, 0, 0, 1, 0, 0,
  0, 0, 0, 0, 0, 0, 1, 0, 0, 0
)

# 軽い入院を1行目、重い死亡を2行目にした日数の行列
Time <- rbind(hosp_day, death_day)

# 同じ並びのイベント有無の行列
Delta <- rbind(hosp_event, death_event)

# 群の指標（対照群は0、宇治茶を飲む群は1）
trt <- rep(c(0, 1), each = 10)

# 共変量：年齢と握力（kg）
age <- c(
  74, 79, 67, 85, 71, 76, 70, 83, 78, 65,
  72, 68, 81, 75, 70, 66, 84, 73, 69, 77
)
grip <- c(
  24, 20, 33, 16, 29, 23, 30, 18, 21, 35,
  28, 32, 19, 25, 30, 34, 17, 27, 31, 22
)

# 参加者を行、共変量を列にした行列
cov <- cbind(age, grip)

# 対照群の人数、治療群の人数、イベントの数
n0 <- sum(trt == 0)
n1 <- sum(trt == 1)
m <- nrow(Time)

カプラン・マイヤー法で状態の確率を推定：kmコマンド

kmコマンドとmarkovコマンドは群の指標を受け取らず、TimeとDeltaの先頭n0列を対照群、残りを治療群として扱います。対照群を先に並べてから実行してください。

# カプラン・マイヤー法で群ごとの状態の確率を推定
z <- km(n0, n1, m, Time, Delta)

# 宇治茶を飲む群のイベント時点
z$unique_event_times1
[1]  95 150 210 240 300 330 365

# 宇治茶を飲む群の状態の確率（行はイベント無し・入院・死亡、列はイベント時点）
round(z$dist_state1, 3)
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,]  0.9  0.8  0.7  0.6  0.5  0.5  0.5    0
[2,]  0.1  0.2  0.3  0.4  0.5  0.4  0.4    0
[3,]  0.0  0.0  0.0  0.0  0.0  0.1  0.1    0

マルコフモデルで状態の確率を推定：markovコマンド

# マルコフモデルで群ごとの状態の確率を推定
z2 <- markov(n0, n1, m, Time, Delta)

# 宇治茶を飲む群の状態の確率をkmコマンドの結果と比べる
round(z2$dist_state1, 3)
     [,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,]  0.9  0.8  0.7  0.6  0.5  0.5  0.5
[2,]  0.1  0.2  0.3  0.4  0.5  0.4  0.4
[3,]  0.0  0.0  0.0  0.0  0.0  0.1  0.1

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

Rで学ぶ統計的データ解析 (データサイエンス入門シリーズ)

データサイエンス入門シリーズの一冊。統計的な考え方の数理を押さえつつ、Rで手を動かして確認できる。

www.amazon.co.jp

win time法で治療効果を推定：wintimeコマンド

主なtypeは次のとおりです。”ewt”は期待win time、”ewtr”は対照群を基準にした期待win time、”rmt”は打ち切り日数までで治療群が有利だった平均日数、”wtr”はwin time比、”pwt”は参加者の組ごとに比べるwin timeです。結果のtrt_effectが治療効果の推定値です。

# 期待win time（EWT）で宇治茶を飲む群が良い状態で過ごせた日数を推定
r_ewt <- wintime("ewt", Time, Delta, trt)

# 治療効果の推定値と、イベントごとの内訳
r_ewt$trt_effect
[1] 118.4

r_ewt$components
[1] 65.4 53.0

# win time比（WTR）をブートストラップ100回でp値まで計算
r_wtr <- wintime("wtr", Time, Delta, trt,
  resample = "boot", resample_num = 100, seed = 2026
)

# 治療効果の推定値、p値、勝ちと負けの数
r_wtr[c("trt_effect", "p", "wins", "losses")]
$trt_effect
[1] 0.328125

$p
[1] 0.01175476

$wins
[1] 64

$losses
[1] 21

# 最初の180日に限った、治療群が有利だった平均日数（RMT）
r_rmt <- wintime("rmt", Time, Delta, trt, time_restriction = 180)
r_rmt$trt_effect
[1] 45.3

# 年齢と握力を共変量に加え、対照群を基準にした期待win time（EWTR）を推定
r_ewtr <- wintime("ewtr", Time, Delta, trt, cov = cov)
r_ewtr$trt_effect
     trt
94.94223

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

基礎から学ぶ統計学実践編〜PCを使った計算とカテゴリカルデータ解析

PCで実際に計算しながら進める実践編。クロス集計などカテゴリカルデータの解析を扱う。

www.amazon.co.jp

Cox比例ハザードモデルで治療効果のz統計量を計算：COMPコマンド

wintimeコマンドのtype = “max”の中で使われるコマンドです。結果は名前の無いリストで、z統計量、治療効果の推定値、その分散、p値の順に入っています。

# 年齢と握力で調整したCox比例ハザードモデルを当てはめる
r_comp <- COMP(length(trt), Time, Delta, cov, trt)
r_comp
[[1]]
     trt
2.160131

[[2]]
      trt
0.1825737

[[3]]
[1] 0.6197907

[[4]]
       trt
0.03076255

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

多変量解析の基礎 II　主成分分析（改訂版）: 理論とRによる演習

主成分分析に絞って理論を丁寧に追い、Rの演習で確認する構成。数式の意味を腹落ちさせたい人向け。

www.amazon.co.jp

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