数学いらずの医科統計学PART6 CHAPTER30をRでやってみる。
生存曲線については、ただひたすらRだけをやる特別講義でも扱っている。
Kirk AP (1980)らは、プレドニゾロンまたはプラセボを投与された、慢性活動性肝炎患者の生存を比較した（データはAltman DG & Bland LMによる）。
pred <- c(2, 6, 12, 54, 56, 68, 89, 96, 96, 125, 128, 131, 140, 141, 143, 145, 146, 148, 162, 168, 173, 181) pred01 <- c(1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0) ctrl <- c(2, 3, 4, 7, 10, 22, 28, 29, 32, 37, 40, 41, 54, 61, 63, 71, 127, 140, 146, 158, 167, 182) ctrl01 <- c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0) g<- cbind(c(pred01, ctrl01), c(pred, ctrl ), c(rep(1, length(ctrl)), rep(2, length(ctrl))) ) colnames(g) <- c("outcome", "time", "group")
観察結果outcomeは、0が研究打ち切りまで生存、1が死亡。
（プレドニゾロン56ヶ月は、研究から脱落したらしいが、R上でどういれておけばいいか謎だったので56ヶ月で打ち切りということにここではしている。）
生存時間timeは月単位。
群はプレドニゾロンが1、プラセボが2。
カプランマイヤー法をやってみる。
library(survival) Surv(g[, 2], g[, 1]) result <- survfit(Surv(time=g[, 2], event=g[, 1]) ~ g[, 3], type="kaplan-meier") plot(result, col=1:2) legend(120, 1, c("predonisolone", "placebo"), col=1:2, lwd=3) abline(h=0.5, lty=2) for(i in 1:2){ segments(summary(result)$table[i, "median"], 0.5, summary(result)$table[i, "median"], -1, lty=2, col=i) } survdiff(Surv(g[, 2], g[, 1]) ~ g[, 3]) Call: survdiff(formula = Surv(g[, 2], g[, 1]) ~ g[, 3]) N Observed Expected (O-E)^2/E (O-E)^2/V g[, 3]=1 22 11 16.4 1.77 4.66 g[, 3]=2 22 16 10.6 2.73 4.66 Chisq= 4.7 on 1 degrees of freedom, p= 0.0309
ログランクテストによる検定では、p=0.031である。
帰無仮説は、二者の治療法が生存全体を変化させず、観察された差は単に偶然の結果である。
これが確率0.03なので、まあ、棄却されるのだろう。
ハザード比
2つの生存曲線はハザード比として要約され、これは本質的に相対危険度と同じらしい。
このデータでは、ハザード比は0.42、95%CIは0.19~0.92である。
プレドニゾロン治療の患者は、どの時間をとってもプラセボ治療患者に比べて、死亡する確率は42%である、ということらしい。
中央生存期間の比
50%生存に水平線を引いた。プレドニゾロンで146、プラセボで40である。これは3.61倍違う。
これも95%CIを計算できるらしく、3.14~4.07らしい。プレドニゾロンによって50%生存が3~4倍延長することは95%確実である。
p値による生存曲線の比較（参考1、参考2）
ログランク検定法ではすべての時点に等しい重みを置く。
これに対して、Gehan-Breslow-Wilcoxon法は早期の時点における死亡に重みを置く。患者の大部分が早期の時点で打ち切られる場合には、誤った結果をもたらしうる。
この方法は比例ハザード比の前提を必要としない。
ひとつの群が他方の群より一貫して高いリスクを持つことを必要とする。
この方法ではp=0.011らしい。Rでまだ見つからない。
おまけ
Rにはデータセットがたくさんある。
heartというデータセットがhelpで使われていた。
heart heartdata <- Surv(heart$start, heart$stop, heart$event) plot(survfit(heartdata ~ heart$surgery), col=1:2) title("Surgery") legend(1200, 1, c("non surgery", "surgery"), col=1:2, lwd=3) plot(survfit(heartdata ~ heart$transplant), col=1:2) title("Transplant") legend(1200, 1, c("non transplant", "transplant"), col=1:2, lwd=3)
（たぶん）普通の心臓手術ありなしと移植ありなしで比較してみた。
