院内がん登録から後ろ向きコホートを組んだとします。2015年1月1日に在院していた患者を登録し、そこから死亡まで追跡して、診断からの5年生存率を出す。データはきれいに揃い、生存率は文献値よりいくぶん良く出ました。査読者から返ってきたのは、こういうコメントです。「Patients who died before the index date could not have been included. Please account for left truncation.」
登録日より前に亡くなった患者は、そもそもデータに入っていません。これは打ち切りではなく、集団の作られ方の問題です。
打ち切りと切断は、どちらも「全部は見えていない」状況を指すため混同されやすいのですが、統計的には別の場所に効いてきます。ここでは両者の違いを尤度の形まで下ろして確認し、遅延参入を無視したときにどれだけずれるかをRで計算します。あわせて、打ち切りの扱いで最も誤りが多い競合イベントの問題も見ていきます。
観察が途切れるのか、集団に入らないのか
打ち切りは、観察対象として登録された個体について、イベント時刻が分からないまま追跡が終わることです。右打ち切りなら「少なくとも$t_i$までは起きなかった」という情報が残ります。生存時間解析入門で見たとおり、この情報は尤度に組み込めます。
切断は違います。ある条件を満たさない個体が、そもそも観察対象に入らない。データを見ても、その人たちが何人いたのかは分かりません。
尤度の形で書くと違いがはっきりします。右打ち切りの個体は$S(t_i)$を寄与させます。一方、左切断された（正確には、左切断を生き延びて観察集団に入った）個体は、$L_i$まで生存していたという条件のもとで観察されているので、条件付き確率になります。
$$\frac{f(t_i)}{S(L_i)} \quad(\text{イベント観察}), \qquad \frac{S(t_i)}{S(L_i)} \quad(\text{右打ち切り}) \tag{1}$$
式(1)の分母$S(L_i)$が補正です。この個体は$L_i$まで生き残った人の中から選ばれているので、その条件付き確率で割り戻さないと、集団全体の分布を推定したことになりません。分母を書き忘れると、長生きした人だけを集めた標本を、無選択の標本のように扱うことになります。
遅延参入を無視するとどれだけずれるか
数式の補正がどれくらい効くのかは、実際に数字を出したほうが早いです。発症からの真の生存時間の中央値を24か月とし、発症から登録までの待ち時間を0から24か月の一様分布で与えます。登録日より前に死亡した人はデータに入りません。以下は説明のための架空データです。
library(survival)
set.seed(42)
n <- 800
onset_to_death <- rexp(n, rate = log(2) / 24) # 発症から死亡まで
onset_to_entry <- runif(n, min = 0, max = 24) # 発症から登録まで
enrolled <- onset_to_death > onset_to_entry # 登録日まで生きていた人だけ
naive <- survfit(Surv(onset_to_death[enrolled]) ~ 1)
delayed <- survfit(Surv(onset_to_entry[enrolled], onset_to_death[enrolled],
rep(1, sum(enrolled))) ~ 1)
summary(naive, times = 24)
summary(delayed, times = 24)#> enrolled: 605 of 800
#>
#> time n.risk n.event survival std.err lower 95% CI upper 95% CI
#> 24 421 184 0.696 0.0187 0.66 0.734
#>
#> time n.risk n.event survival std.err lower 95% CI upper 95% CI
#> 24 421 184 0.499 0.0335 0.437 0.56924か月時点の真の生存率は0.500です。参入時点を無視した推定は0.696、遅延参入を入れた推定は0.499でした。生存例と死亡例の数はどちらも同じで、違うのは参入時点をSurv関数の第1引数に渡したかどうかだけです。
この20ポイントの差は標本を増やしても縮みません。800例のうち195例が登録前に死亡して消えており、残った605例は定義上「登録日まで生き延びた人」だからです。
plot(naive, conf.int = FALSE, xlab = "Months since onset", ylab = "S(t)")
lines(delayed, conf.int = FALSE)
curve(exp(-log(2) * x / 24), from = 0, to = 60, lty = 2, add = TRUE)実務でこの補正が要るかどうかの見分け方は単純です。時間の原点（診断日、発症日、手術日）と、コホートに入った日が一致しないなら、左切断が入っています。原点が登録日そのものなら不要です。院内がん登録や保険レセプトのデータベース研究、既存コホートの二次利用では、原点が登録日より前にあることがほとんどなので、まず疑ってよい状況になります。
Rでは参入時点、脱落時点、イベント指示変数の3つをSurv関数に渡す形で対応します。Cox回帰でも同じ書き方が使えるので、対応そのものは重くありません。私なら、後ろ向きデータを受け取った時点で「時間の原点はいつか」を最初に確認します。ここを決めないまま解析を始めると、あとで全部やり直しになります。
他病死を打ち切りにしてはいけない
打ち切りの扱いで、臨床研究の論文に最も多く残っている誤りがここです。がん死亡の累積発生率を出すときに、他病死を右打ち切りとして処理する。手順としては自然に見えますし、Rも警告を出しません。
ところが右打ち切りは「その後もイベントが起こりうる」という前提を含んでいます。心疾患で亡くなった患者に、そのあとがん死亡が起きることはありません。前提が事実と食い違っています。
set.seed(42)
m <- 600
t_cancer <- rexp(m, rate = 0.020) # がん死亡までの時間
t_other <- rexp(m, rate = 0.035) # 他病死までの時間
obs <- pmin(t_cancer, t_other, 60)
cause <- ifelse(obs == 60, 0, ifelse(t_cancer < t_other, 1, 2))
cr <- survfit(Surv(obs, factor(cause)) ~ 1) # Aalen-Johansen推定
km <- survfit(Surv(obs, cause == 1) ~ 1) # 他病死を打ち切り扱い
summary(cr, times = 36)
summary(km, times = 36)#> time n.risk n.event Pr((s0)) Pr(1) Pr(2)
#> 36 96 504 0.16 0.297 0.543
#>
#> time n.risk n.event survival std.err lower 95% CI upper 95% CI
#> 36 96 178 0.527 0.029 0.473 0.58736か月時点のがん死亡の累積発生率は、競合リスクを織り込んだAalen-Johansen推定で0.297です。他病死を打ち切りにして1からKaplan-Meier推定を引くと0.473になります。1.6倍です。
差が生まれる理由は、打ち切り扱いにした患者がリスク集合から抜けたあとも、生存関数の計算上は「これからがんで死ぬ可能性のある人」として扱われ続けることにあります。高齢者コホートのように他病死が多い集団ほど、この乖離は大きくなります。
原因別の累積発生率を報告するなら、Aalen-Johansen推定を使うことを勧めます。survivalパッケージでは、イベント種別をfactorにしてSurv関数に渡すだけで計算できるので、手間はほとんど変わりません。
ただし、他病死を打ち切りにした解析が常に間違いというわけではありません。原因別ハザードを推定したい場合、つまり「この治療はがん死亡のハザードを下げるか」を問う場合は、他病死を打ち切りとして扱うのが正しい定式化です。累積発生率という「実際に何割の患者に起きるか」を答えたいのか、ハザードという「その時点でのリスクの強さ」を答えたいのか。どちらを報告するのかを先に決めれば、扱いは自動的に決まります。Latouche et al.（J Clin Epidemiol 2013;66:648-653）は、この2つを混同した報告が多いことを指摘しています。
検診で見つかるイベントは区間打ち切りになる
糖尿病網膜症の発症、C型肝炎後の肝硬変への進展、関節リウマチの骨びらん。定期検査で判定するアウトカムは、発生時点が正確には分かりません。前回の検査で陰性、今回で陽性なら、その間のどこかで起きたことだけが確定します。
実務では、この区間の中点を発症日として通常の生存時間解析にかける処理がよく行われます。検査間隔が1か月なら実害はほぼありませんが、年1回の検診で12か月幅の区間を中点で埋めると、発症時点の誤差が最大6か月になります。群間で検査間隔が違う研究では、この誤差が群間差として現れることがあります。
区間$[L_i, R_i]$でイベントが起きた確率は、生存関数の差として書けます。
$$\Pr(L_i \le T < R_i) = S(L_i) – S(R_i) \tag{2}$$
式(2)は、$L_i$まで生き延びた確率から$R_i$まで生き延びた確率を引けば、その間にイベントが起きた確率になる、というだけの関係です。中点補完はこの区間を1点に潰す近似で、区間が広いほど近似の粗さが効いてきます。
survivalパッケージでは、区間の下限と上限を渡してTurnbull推定量を求められます。検査間隔が長い研究や、群間で間隔が揃っていない研究では、中点補完の結果と併記して感度分析にすることを勧めます。両者が一致すれば近似で問題ないと書けますし、ずれれば区間打ち切りを正面から扱う根拠になります。
無情報打ち切りは検証できない
ここまでの補正はすべて、打ち切りがイベントの起きやすさと無関係に生じることを前提にしています。共変量$X$を条件として、生存時間$T$と打ち切り時間$C$が条件付き独立である、という仮定です。
問題は、この仮定がデータから確かめられないことにあります。脱落した患者のその後の生存時間は観察されていないので、脱落と予後の関係を推定する材料がありません。有害事象が重い患者ほど来院しなくなる、という状況を否定する統計量は作れません。
できるのは、仮定が崩れる方向を具体的に想定して結果がどこまで動くかを見ることです。ICH E9(R1)の推定目標（estimand）の枠組みでは、脱落を含む中間事象をどう扱うかを解析計画の段階で書くことが求められています。脱落理由を「有害事象」「効果不十分」「転居」「その他」に分けて集計し、予後と関連しそうな理由が偏っていないかを示す。そのうえで、脱落例が最悪の経過をたどったと仮定した場合の推定値を併記する。ここまで書いてあれば、査読で仮定の妥当性を問われても答えられます。
逆に、「打ち切りは無情報と仮定した」とだけ書いて理由を示さない原稿は、指摘を受けやすくなります。仮定そのものが悪いのではなく、検証できない仮定を検証したかのように書くことが問題になります。
解析に入る前に決めておくこと
データを受け取ったら、時間の原点をいつに置くかを最初に決めます。診断日か、治療開始日か、コホート登録日か。原点が登録日より前にあれば左切断が入るので、参入時点を解析に渡す必要があります。次に、イベントの定義と、それを妨げる別のイベントがあるかどうか。他病死や移植のように主要イベントの発生を止める事象があれば、累積発生率とハザードのどちらを報告するのかを決めます。最後に、脱落の理由が記録されているかを確認します。
この3点が決まっていれば、Kaplan-Meier法に進んでよい状態になります。推定量そのものの性質と、信頼区間の作り方はKaplan-Meier推定量の回で扱います。ハザード関数の側から見た定式化は生存関数とハザード関数、共変量調整に進む場合はCox比例ハザードモデルが続きになります。
後ろ向きのデータベース研究で、時間の原点をどこに置くべきか、脱落理由をどう集計して本文に書くか。このあたりは研究デザインとデータの成り立ちの両方を見ないと決まりません。Dr.データサイエンスでは、こうした解析方針のご相談を承っています。
参考文献
Klein JP, Moeschberger ML. Survival Analysis: Techniques for Censored and Truncated Data. 2nd ed. Springer; 2003.
Turnbull BW. The empirical distribution function with arbitrarily grouped, censored and truncated data. Journal of the Royal Statistical Society: Series B. 1976;38(3):290-295.
Latouche A, Allignol A, Beyersmann J, et al. A competing risks analysis should report results on all cause-specific hazards and cumulative incidence functions. Journal of Clinical Epidemiology. 2013;66(6):648-653.
Andersen PK, Geskus RB, de Witte T, Putter H. Competing risks in epidemiology: possibilities and pitfalls. International Journal of Epidemiology. 2012;41(3):861-870.
ICH E9(R1). Addendum on Estimands and Sensitivity Analysis in Clinical Trials. International Council for Harmonisation; 2019.
