Top 5 This Week

関連記事

【医療統計版】論文用の図表と報告:Kaplan–Meier図・リスク集合・フォレストプロット

- 本サイト運営者のサービスの紹介 -

生存時間解析の結果を論文にまとめる段階で、原稿に載る図は生存曲線1枚とハザード比1つ、という形になりがちです。ところが査読で追加を求められるのは、たいていその周辺にある情報のほうで、各時点で何人が追跡下に残っていたのか、打ち切りはいつ起きたのか、部分集団の違いをどの推定量で判断したのか、といった点が問われます。曲線そのものは正しく計算されていても、読者が誤読しない図になっているかは別の問題です。この記事では、生存時間解析の結果を報告するときに図表へ何を配置するかを、survivalパッケージのveteranデータを使って具体的に決めていきます。

扱うのは、Kaplan–Meier図に添えるリスク集合の人数と打ち切り印と信頼区間、群間比較で書く効果量と区間推定、事前に定めた部分集団を図にするときの交互作用です。読み終えたときに、自分の解析結果のどこを図にしてどこを表や本文に回すか、その配分を判断できる状態を目指します。Kaplan–Meier推定量そのものの計算や比例ハザード性の診断はKaplan–Meier推定量と比例ハザード仮定の検証と診断で扱っているので、ここでは推定が済んだ後の報告に絞ります。

使うデータは肺がん患者137例の無作為化試験記録で、標準治療群69例、試験治療群68例、観察されたイベントは128件です。打ち切りは9例(6.6%)にとどまります。以下では標準治療群をStandard、試験治療群をTestと表記し、$S(t)$を時点$t$までイベントが起きない確率、$n_j$を時点$t_j$の直前に追跡下にある人数、$d_j$をその時点のイベント数とします。

生存曲線だけでは伝わらない情報

Kaplan–Meier推定量は、イベントが起きた各時点での条件付き生存確率を掛け合わせた形をしています。

$$\hat{S}(t)=\prod_{t_j \le t}\left(1-\frac{d_j}{n_j}\right) \tag{1}$$

この式は、時点$t$までに観測されたイベント時点を順にたどり、それぞれの時点で「直前に追跡下にいた$n_j$人のうち$d_j$人にイベントが起きなかった割合」を掛け合わせていく、という計算を表しています。曲線の高さは$n_j$に依存して決まりますが、描かれた線を見ても$n_j$がいくつだったかは分かりません。veteranデータでは360日時点で追跡下にいるのがStandard群4例、Test群6例まで減っていますが、曲線の右端はそのことを何も語りません。線の太さも滑らかさも、137例の部分と10例の部分で変わらないからです。

そのため、生存曲線を1枚だけ載せる報告では、読者は末端の上下差を前半と同じ確からしさで受け取ります。

リスク集合の人数が果たす役割

veteranデータのStandard群では、90日時点の生存確率が54.7%(95%信頼区間42.2〜65.6%)、そのときのリスク集合が37例です。同じ群の360日時点では生存確率7.1%(2.3〜15.5%)で、リスク集合は4例まで落ちています。区間の幅が23ポイントから13ポイントへ狭まって見えるのは推定が安定したからではなく、推定値そのものが0に近づいたためで、相対的な不確実性はむしろ大きくなっています。人数を併記しておけば、読者はこの違いを自分で判断できます。

Morrisら(2019)のKMunicate研究は、Kaplan–Meier図の表示方法について研究者と統計家に意見を集め、曲線の下にリスク集合の人数を並べた形式が支持されることを報告しています。リスク集合の行を加えると、曲線の各部分を支える人数が明示されます。表示時点は追跡期間とイベントの分布に合わせて定め、人数の減少過程が読み取れる間隔にします。veteranデータでは30日刻みから始め、180日以降の間隔を広げました。

人数を載せない図を勧めない理由は、末端の解釈が読者任せになるからです。末端のリスク集合が小さくなるほど、曲線の段差は少数のイベントに敏感になります。したがって、交差の見た目だけで群差が安定しているとは判断できません。

標準治療群と試験治療群のKaplan–Meier曲線に半透明の95%信頼区間と打ち切り印を重ね、下に各時点のリスク集合の人数を並べた図
図1 veteranデータの治療群別Kaplan–Meier曲線。半透明の帯が対数対数変換に基づく95%信頼区間、+印が打ち切り時点で、下段の数字が各時点で追跡下にいる人数。180日を過ぎると両群とも15例前後まで減り、その後の上下差は少数例で決まっている。

打ち切り印と信頼区間の置き方

打ち切り印を省いても、Kaplan–Meier推定値そのものは1つも動きません。式(1)の計算には打ち切りが$n_j$の減少として織り込まれているからです。それでも印を消すと、曲線の平坦な区間が「その間イベントが起きなかった」ためなのか「追跡から外れた人がいて情報が薄い」ためなのかを区別できなくなり、後半の情報量を実際より多く見せる方向に働きます。打ち切りが図の読み方へ与える影響は割合だけで決まらず、どの時点で起きたかにも依存します。

曲線を2本並べただけでは、その差が推定の揺れに収まる程度なのかどうかは判断できません。veteranデータの90日時点はStandard群54.7%、Test群38.0%(26.6〜49.4%)で、点推定の差は約17パーセントポイントですが、両群の区間は重なっています。

もう一つ注意したいのが縦軸の範囲です。生存確率の軸を0から1にせず、差が出ている帯だけを拡大すると、同じデータが別の主張を持った図に変わります。

同じveteranデータのKaplan–Meier曲線を2通りに表示した図。左は縦軸を0から0.35に切り打ち切り印と信頼区間を省いたもの、右は縦軸を0から1にとり打ち切り印と半透明の95%信頼区間を表示したもの
図2 同一の推定結果を2通りに表示したもの。パネルAは縦軸を0.35までに切り、打ち切りと信頼区間を省いたため、200日以降の群差だけが目立つ。パネルBでは95%信頼区間を半透明の帯で示しており、同じ差が全体の中では小さく、360日を過ぎた領域に残るのが4例と6例であることも分かる。

パネルAとパネルBの数値は完全に同じで、違うのは軸の範囲と追跡状況の表示の有無だけです。そのため、縦軸の範囲は群間差を強調する目的で決めるのではなく、生存確率の全体像と差の大きさを同時に読めるように定めます。軸の切断が必要な場合は、切断した範囲を明示し、全範囲の図または主要時点の数値を併記すると表示による印象の偏りを抑えられます。

中央値とハザード比が食い違う場面

veteranデータの生存期間中央値は、Standard群103日(95%信頼区間54〜126日)、Test群52.5日(43〜90日)です。Test群のほうが半分程度に短く、群間差があるように読めます。ところが全期間を1つの係数でまとめたCoxモデルのハザード比は1.02(0.71〜1.45、p=0.92)、ログランク検定でもp=0.93で、群差を支持する結果になりません。

この食い違いは計算の誤りではなく、2つの要約が別の情報を使っているために起きます。中央値は生存確率が0.5を切る1点だけで決まり、図1のように曲線が130日付近で入れ替わる場合、交差の後ろで起きたことを一切反映しません。一方の全期間のハザード比は、どの時点でどれだけイベントが起き、その時点のリスク集合がどう構成されていたかを反映した、追跡期間全体の要約値になります。前半のハザード比と後半のハザード比を単純に平均した値ではないので、「中央値では差があるがハザード比では差がない」という状態は矛盾ではありません。

ここで書いてはいけないのが、p=0.92を根拠に「両群の生存に差はなかった」と結論することです。信頼区間0.71〜1.45は、Test群のハザードが低い向きと高い向きの両方を含み、効果の向きと大きさに広い不確実性を残しています。有意でないことは差がないことの証明にならないので、報告では効果量と区間を並べて書き、検定結果だけで締めないほうが読者に伝わります。CONSORT 2025声明も、群ごとの結果に加えて効果量と精度を示すことを求めています。曲線が交差している研究では、期間を分けた推定や制限付き平均生存時間の差を併記する方法もあり、この判断にはCox比例ハザードモデルの前提が保たれているかどうかが関わります。

部分集団の図と交互作用の推定

Test群とStandard群の比較を、既治療歴の有無で分けてみます。この部分集団は解析を見てから選んだのではなく、データを開く前にpriorという変数1つに固定したものです。既治療歴なしの97例(イベント91件)ではハザード比1.25(0.82〜1.90、p=0.29)、ありの40例(イベント37件)では0.67(0.34〜1.33、p=0.25)となり、点推定の向きが逆になります。

向きが逆でも、どちらのp値も0.05を超えています。ここで「片方だけ効いている」と読むのはよくある誤りで、Altmanら(2003)が指摘したとおり、一方が有意で他方が非有意という状況は2つの効果が異なることの検定になっていません。効果の差を知りたいなら、差そのものを推定します。既治療歴で層別したCoxモデルに交互作用項を入れると、ハザード比の比が推定できます。

$$\theta=\frac{\mathrm{HR}_{\text{prior=yes}}}{\mathrm{HR}_{\text{prior=no}}}=\exp(\beta_{\mathrm{int}}) \tag{2}$$

$\theta$は既治療歴ありの群で得られる治療のハザード比が、なしの群のハザード比の何倍にあたるかを表し、$\beta_{\mathrm{int}}$は交互作用項の係数です。$\theta$が1なら、治療効果は既治療歴によらず同じ大きさになります。veteranデータでは$\theta$が0.54(0.24〜1.19、p=0.13)で、点推定は既治療歴ありのほうで治療が有利な向きを示しているものの、区間は1を大きくまたいでいます。部分集団ごとのp値を見比べたときの印象より、判断できることはずっと少ないと分かります。

Pocockら(2002)は主要な医学雑誌に載った無作為化試験の報告を調べ、部分集団解析の結果が交互作用の検定を伴わずに示されている例が多いことを報告しました。フォレストプロットは行を増やすほど情報が多いように見えますが、Cuzick(2005)が述べるように、各行を独立した結論として読ませる図ではありません。全体効果と交互作用を同じ図に置くのは、行ごとの読み方を防ぐためです。

全体、既治療歴なし、既治療歴ありのハザード比と95%信頼区間、および交互作用としてのハザード比の比を対数軸に並べたフォレストプロット
図3 治療効果の全体推定と、事前に固定した部分集団別の推定、および両者の差にあたるハザード比の比。部分集団の点推定は1を挟んで反対側にあるが、最下行の交互作用は0.54で95%信頼区間が0.24〜1.19となり、効果が異なると判断できる状態にはない。

報告ガイドラインが求める記載

報告に必要な要素は、図の慣習ではなく文書化された基準から決めるほうが確実です。CONSORT 2025声明は無作為化試験の報告項目を定めていて、主要アウトカムについて群ごとの結果と、群間差の効果量および精度を示すことを求めています。生存時間解析なら、ハザード比や制限付き平均生存時間の差と、その95%信頼区間がこれに当たります。部分集団については、事前に規定したものと探索的なものを区別して報告することが求められます。

追跡期間の記載も抜けやすい項目です。最終イベントの時点を追跡期間として書いてしまう原稿を見かけますが、これは追跡の長さではなくイベントの起きた範囲を表しているだけになります。veteranデータでは、打ち切りをイベントとして扱う逆Kaplan–Meier曲線が0.5へ到達しないため、追跡期間中央値は未到達です。この結果は未到達と報告し、観察時間の中央値で置き換えることは適切ではありません。補足が必要なら、推定可能な追跡期間分布の分位点、追跡期間の定義に沿った範囲、データカット日などから目的に合う量を選び、その定義を明示します。

必要イベント数の考え方は必要イベント数:検出力は症例数だけでは決まらないで扱っていて、報告時に例数だけでなくイベント数を書く理由もそこにつながります。図3の各行に例数とイベント数を併記したのも、既治療歴ありの40例でイベントが37件と、情報量が例数の印象ほど小さくないためです。

Rでの作図と出力の読み方

次の実行済みコードは、本文に残したKaplan–Meier推定値、中央値、Coxモデル、ログランク検定、部分集団別推定、交互作用、逆Kaplan–Meier結果を同じデータから生成します。準備処理を共有し、その後に図1から図3を作ります。

library(survival); d <- survival::veteran
d$trt2 <- factor(d$trt, 1:2, c("Standard", "Test")); d$prior2 <- factor(d$prior, c(0,10), c("No","Yes"))
d$test <- as.integer(d$trt2 == "Test"); d$prior_yes <- as.integer(d$prior2 == "Yes")
km <- survfit(Surv(time,status) ~ trt2, d, conf.type="log-log"); tt <- c(0,30,60,90,120,180,360)
ko <- summary(km, times=tt, extend=TRUE); med <- summary(km)$table[,c("median","0.95LCL","0.95UCL")]
overall <- coxph(Surv(time,status) ~ trt2, d, ties="efron"); lr <- survdiff(Surv(time,status) ~ trt2, d)
sub <- lapply(levels(d$prior2), \(z) coxph(Surv(time,status) ~ trt2, d[d$prior2==z,], ties="efron"))
inter <- coxph(Surv(time,status) ~ test + test:prior_yes + strata(prior2), d, ties="efron")
reverse_km <- survfit(Surv(time,1-status) ~ 1, d, conf.type="log-log")
effect <- \(fit,term) {s <- summary(fit)$coefficients; c(HR=exp(s[term,"coef"]),L=exp(s[term,"coef"]-1.96*s[term,"se(coef)"]),U=exp(s[term,"coef"]+1.96*s[term,"se(coef)"]),p=s[term,"Pr(>|z|)"])}
c(n=nrow(d),events=sum(d$status),censored=sum(!d$status)); med
summary(km, times=c(30,90,180,360), extend=TRUE)[c("time","n.risk","surv","lower","upper")]
effect(overall,"trt2Test"); c(logrank_p=pchisq(lr$chisq,1,lower.tail=FALSE))
e_no <- effect(sub[[1]],"trt2Test"); e_yes <- effect(sub[[2]],"trt2Test")
e_int <- effect(inter,"test:prior_yes"); rbind(No=e_no,Yes=e_yes,Interaction=e_int)
summary(reverse_km)$table[c("median","0.95LCL","0.95UCL")]
rows <- rbind(effect(overall,"trt2Test")[1:3],e_no[1:3],e_yes[1:3],e_int[1:3])
#> Overall HR 1.02 (0.71–1.45), p=0.92; log-rank p=0.93
#> No prior HR 1.25 (0.82–1.90), p=0.29
#> Prior therapy HR 0.67 (0.34–1.33), p=0.25
#> Interaction ratio 0.54 (0.24–1.19), p=0.13
#> Reverse Kaplan–Meier median: NA (not reached)

逆Kaplan–Meier出力のNAは、曲線が0.5へ到達せず中央値を推定できないことを示します。図1の数値はkoのn.risk、曲線と外側の線はsurv、lower、upperに対応します。

png("fig1-km-risk-table.png",1100,760,res=200); layout(matrix(1:2,2),heights=c(3,1))
par(mar=c(3,7,1,1)); plot(km,conf.int=FALSE,mark.time=FALSE,col=NA,xlim=c(0,400),ylim=0:1,xaxt="n",xlab="Time since randomisation (days)",ylab="Overall survival probability")
stepxy <- function(t,y) list(x=c(t[1],rep(t[-1],each=2)),y=c(rep(y[-length(y)],each=2),tail(y,1)))
e <- cumsum(km$strata); b <- c(1,head(e,-1)+1)
for(i in 1:2){j <- b[i]:e[i]; j <- j[km$time[j]<=400]; z <- c(0,km$time[j],400)
lo <- stepxy(z,c(1,km$lower[j],tail(km$lower[j],1))); hi <- stepxy(z,c(1,km$upper[j],tail(km$upper[j],1)))
polygon(c(lo$x,rev(hi$x)),c(lo$y,rev(hi$y)),border=NA,col=adjustcolor(c("#1f4e79","#b45309")[i],alpha.f=.18))}
lines(km,conf.int=FALSE,mark.time=TRUE,col=c("#1f4e79","#b45309"),lwd=2); axis(1,tt)
legend("topright",c("Standard","Test"),col=c("#1f4e79","#b45309"),lwd=2,bty="n")
par(mar=c(1,7,0,1)); plot.new(); plot.window(c(0,400),c(0,3)); text(-90,2.5,"No. at risk",adj=0,xpd=NA)
for(i in 1:2){j <- ko$strata==levels(ko$strata)[i]; text(-90,2-i*.6,levels(d$trt2)[i],adj=0,xpd=NA); text(ko$time[j],2-i*.6,ko$n.risk[j])}
dev.off()
png("fig2-axis-and-censoring.png",1100,620,res=200); par(mfrow=1:2)
plot(km,conf.int=FALSE,mark.time=FALSE,col=c("#1f4e79","#b45309"),lwd=2,xlim=c(0,999),ylim=c(0,.35),xlab="Time (days)",ylab="Survival probability",main="A. Truncated axis, no censoring, no CI")
plot(km,conf.int=FALSE,mark.time=FALSE,col=NA,xlim=c(0,999),ylim=0:1,xlab="Time (days)",ylab="Survival probability",main="B. Full axis, censoring marks, 95% CI")
for(i in 1:2){j <- b[i]:e[i]; j <- j[is.finite(km$lower[j])&is.finite(km$upper[j])]; z <- c(0,km$time[j],999)
lo <- stepxy(z,c(1,km$lower[j],tail(km$lower[j],1))); hi <- stepxy(z,c(1,km$upper[j],tail(km$upper[j],1))); polygon(c(lo$x,rev(hi$x)),c(lo$y,rev(hi$y)),border=NA,col=adjustcolor(c("#1f4e79","#b45309")[i],alpha.f=.18))}
lines(km,conf.int=FALSE,mark.time=TRUE,col=c("#1f4e79","#b45309"),lwd=2); abline(v=360,lty=3); text(370,.85,"Beyond 360 days:\n4 and 6 at risk",adj=0); dev.off()
png("fig3-forest.png",1100,620,res=200); y <- c(4,3,2,.8); plot.new(); plot.window(log(c(.15,9)),c(.2,4.8))
abline(v=0,lty=2); segments(log(rows[,2]),y,log(rows[,3]),y,lwd=2); points(log(rows[,1]),y,pch=c(15,16,16,18),cex=1.5)
axis(1,at=log(c(.25,.5,1,2,4)),labels=c(.25,.5,1,2,4)); title(xlab="Estimate on ratio scale (log scale; null = 1)")
lab <- c("Overall (n=137, 128 events)","No prior therapy (n=97, 91 events)","Prior therapy (n=40, 37 events)","Interaction (ratio of HRs)")
text(log(.15),y+.3,lab,adj=0,xpd=NA); text(log(9),y+.3,sprintf("%.2f (%.2f-%.2f)",rows[,1],rows[,2],rows[,3]),adj=1,xpd=NA); dev.off()

図3の横軸は、上の3行にあるTest対Standardの治療ハザード比と、最下行にある交互作用のハザード比の比をともに含む比尺度です。破線の1が両方の推定量に共通する帰無値になります。

図表の配分を決めるときの判断

どの情報を図に載せ、どれを本文や表に回すかは、読者が図だけを見て誤った結論に至る経路を一つずつ塞ぐ作業になります。veteranデータで確かめたように、リスク集合の人数がなければ末端の差を過大に読み、打ち切り印がなければ平坦部の意味を取り違え、部分集団のp値だけを並べれば効果の違いを推定していないのに違いを語ってしまいます。

主要アウトカムの効果量と95%信頼区間は本文または表に示し、その推定を支える追跡状況をKaplan–Meier図とリスク集合の行で補います。部分集団を示す場合は、事前に定めた区分と交互作用の推定を対応させます。多数の部分集団を並べたフォレストプロットを主要な図に据える方法は推奨しません。行ごとの推定を独立した結論として読ませやすく、Cuzick(2005)が指摘した解釈上の問題につながるためです。

推定対象の設定まで戻って考える必要がある場合は、推定対象と中間事象:治療中止後のイベントをどう扱うかで扱った内容が出発点になります。図の作り込みより先に、その図が何を推定した結果なのかが決まっている必要があります。

図に載せる範囲は、追跡状況、打ち切りの発生時期、事前に規定した部分集団、準拠する報告指針によって変わります。ただし、図だけで推定結果を完結させる必要はありません。図には時間的な経過と不確実性を示し、正確な効果量と95%信頼区間は本文または表で補うと、それぞれの役割が明確になります。

参考文献

  1. Hopewell S, Chan AW, Collins GS, Hróbjartsson A, Moher D, Schulz KF, et al. CONSORT 2025 statement: updated guideline for reporting randomised trials. Lancet. 2025.
  2. Morris TP, Jarvis CI, Cragg W, Phillips PPJ, Choodari-Oskooei B, Sydes MR. Proposals on Kaplan-Meier plots in medical research and a survey of stakeholder views: KMunicate. BMJ Open. 2019;9(9):e030215.
  3. Pocock SJ, Assmann SE, Enos LE, Kasten LE. Subgroup analysis, covariate adjustment and baseline comparisons in clinical trial reporting: current practice and problems. Statistics in Medicine. 2002;21(19):2917-2930.
  4. Altman DG, Bland JM. Interaction revisited: the difference between two estimates. BMJ. 2003;326(7382):219.
  5. Cuzick J. Forest plots and the interpretation of subgroups. Lancet. 2005;365(9467):1308.

前の記事:【医療統計版】40. RMST差の推移:制限時点によって群間差はどう変わるか

Popular Articles