生存時間データを、イベントの種類と1人あたり許容される回数をもとに4つに分けることは述べた。
ここでは、1種類のイベントを2回以上経験しうる状況（反復型生存時間データ, repeated type survival data）について掘り下げたい。 例えば、増悪・寛解や入退院などがこの種のイベントに該当する。
最終的に必要なデータ構造
反復イベント型の生存時間データを解析する際には、1人のデータを複数に区切ったカタチにする必要がある。 この構造は、説明変数がベースライン時点から一定ではなく、時間経過によって変化する「時間依存性共変量（time-dependent covariates）」を含んだデータを扱うときも用いられる。
下に例を示す。各変数は以下のとおり。
- ID：患者識別番号
- time1：観察開始時点（ここでは日単位としておく）
- time2：観察終了時点
- status：観察終了時点でのイベント発生有無（1=イベント発生, 0=打ち切り）
- covariate：観察開始時点の共変量（=イベント発生リスク因子）
この例では1人の対象者の経過を以下の4つの区間に分けている。
- 1行目：観察開始時は共変量が1.2で、15日目まで観察した時点でイベントは発生していなかった。
- 2行目：15日目に共変量が1.5に変化した。46日目時点ではイベント発生なし。
- 3行目：46日目に共変量が1.6に変化した。73日時点でイベントが発生した。
- 4行目：共変量は1.6のままだが、91日目に2回目のイベントが発生した。
繰り返しになるが、 「共変量は区間開始時点での値、イベント発生状況は区間終了時点での値」である。
tmerge( )でデータ整形する
tmerge( )を使って2つのデータフレームを結合（merge）させることで、step by stepに目的の形にしていく。 コードは次の要領で書く。
newdata <- tmerge(data1, data2, id, newvar=tdc(time, value))
data1：元となるデータ。このdata1に、data2の情報を新しい変数として追加していく。data2：追加する情報を含んだデータid：データ同士を紐つけるキーを指定するnewvar=tdc(time, value)：data2に含まれるtimeとvalueを指定する。新しいデータセット内での変数名をnewvarとして指定する。
3つのデータを用意する
以下の3つのデータフレームを用意する。
d_base：ベースラインの患者背景データid：患者IDz：時間非依存性変数（time-constant variable）futime：フォローアップ期間
d_tdcov：時間依存性共変量のデータid：患者IDctime：変数を測定した時点value：そのときの値
d_event：イベント発生データid：患者IDetime：イベント発生時点event_type：イベントの種類（1種類しかなければ不要）*1
それぞれのデータの中身は以下のとおり（いずれも手入力で作成）。
> d_base id z futime 1 1 60 100 2 2 70 120
> d_tdcov id ctime value 1 1 0 1.2 2 1 15 1.5 3 1 46 1.6 4 2 0 0.9 5 2 30 1.0 6 2 80 1.2
患者1は1.2 → 1.5 → 1.6と変化し、患者2は0.9 → 1.0 → 1.2と変化したことを示している。
> d_event id etime event_type 1 1 73 1 2 1 91 1 3 2 20 1 4 2 120 2
患者1は73日と91日にtype1のイベントを発生し、患者2は20日にtype1のイベント、120日にtype2のイベントを発生したことを示している。
tstart, tstopを追加する
d_baseにd_baseを結合させる際に、futimeをtstop（観察終了時点）に指定する。tstart（観察開始時点）も自動的に作成される。
> data1 <- tmerge(d_base, d_base, id=id, tstop=futime) > data1 id z futime tstart tstop 1 1 60 100 0 100 2 2 70 120 0 120
tmerge( )は使わずにtstartとtstopを作成してもよい。
時間依存性共変量データを追加する
上で作成したdata1に、共変量の時間変化を示すデータd_tdcovを追加する。tdc( )（time-dependent covariates）を使うと追加するものが共変量であることが示される。こうするとtdc( )の中のtimeが区間開始時点での値になる。
> data2 <- tmerge(data1, d_tdcov, id=id, cov=tdc(ctime, value)) > data2 id z futime tstart tstop cov 1 1 60 100 0 15 1.2 2 1 60 100 15 46 1.5 3 1 60 100 46 100 1.6 4 2 70 120 0 30 0.9 5 2 70 120 30 80 1.0 6 2 70 120 80 120 1.2
イベント発生データを追加する
上で作成したdata2に、イベント発生を示すデータd_eventを追加する。event( )を使うと追加するものがイベントであることが示される。こうするとevent( )の中のtimeが区間終了時点での発生状況と認識される。
> data3 <- tmerge(data2, d_event, id=id, status=event(etime, event_type)) > data3 id z futime tstart tstop cov status 1 1 60 100 0 15 1.2 0 2 1 60 100 15 46 1.5 0 3 1 60 100 46 73 1.6 1 4 1 60 100 73 91 1.6 1 5 1 60 100 91 100 1.6 0 6 2 70 120 0 20 0.9 1 7 2 70 120 20 30 0.9 0 8 2 70 120 30 80 1.0 0 9 2 70 120 80 120 1.2 2
累積回数を追加する
必要なら各患者における累積回数を追加できる。
それぞれの患者での観察期間数を追加したければ、cumtdc(tstart)とすればよい（cumtdc(tstop)とすると、「それまでに観察期間がいくつあったか」になる）。
> data4 <- tmerge(data3, data3, id=id, cumobs= cumtdc(tstart)) > data4 id z futime tstart tstop cov status cumobs 1 1 60 100 0 15 1.2 0 1 2 1 60 100 15 46 1.5 0 2 3 1 60 100 46 73 1.6 1 3 4 1 60 100 73 91 1.6 1 4 5 1 60 100 91 100 1.6 0 5 6 2 70 120 0 20 0.9 1 1 7 2 70 120 20 30 0.9 0 2 8 2 70 120 30 80 1.0 0 3 9 2 70 120 80 120 1.2 2 4
時間依存性共変量が変化した累積回数はcumtdc(time)で追加する。
（cumtdc(time, value)とすると値の累積値が計算されてしまうので注意。）
> data5 <- tmerge(data3, d_tdcov, id=id, cumcov=cumtdc(ctime)) > data5 id z futime tstart tstop cov status cumcov 1 1 60 100 0 15 1.2 0 1 2 1 60 100 15 46 1.5 0 2 3 1 60 100 46 73 1.6 1 3 4 1 60 100 73 91 1.6 1 3 5 1 60 100 91 100 1.6 0 3 6 2 70 120 0 20 0.9 1 1 7 2 70 120 20 30 0.9 0 1 8 2 70 120 30 80 1.0 0 2 9 2 70 120 80 120 1.2 2 3
「発生したイベントが、その人にとって何回目のイベントか」を追加したい場合はcumevent(time)で追加する。
> data6 <- tmerge(data3, d_event, id=id, cumevent=cumevent(etime)) > data6 id z futime tstart tstop cov status cumevent 1 1 60 100 0 15 1.2 0 0 2 1 60 100 15 46 1.5 0 0 3 1 60 100 46 73 1.6 1 1 4 1 60 100 73 91 1.6 1 2 5 1 60 100 91 100 1.6 0 0 6 2 70 120 0 20 0.9 1 1 7 2 70 120 20 30 0.9 0 0 8 2 70 120 30 80 1.0 0 0 9 2 70 120 80 120 1.2 2 2
tdc( ), cumtdc( )は開始時点に、event( ), cumevent( )は終了時点に対して累積数をくっつける。
例えばcumevent(tstop)は、前述のcumtdc(tstart)と同じものを出力する。
あと、それぞれにおいてvalueやstatusは省略できても、timeは省略できない。ちょっと直感的に分かりにくいときもあるので、実際はやってみて確認することになりそう。
Cox比例ハザードモデルで解析する
慢性肉芽腫症（chronic granulomatous disease, CGD）における反復感染に対して、インターフェロン治療とプラセボを比較したデータを使う。
データcgd0は感染イベントが起こるごとにetimeを変数として横に追加していった横長データ（=1患者1行）。出力すると横長で見にくいのでhead( )で確認してみてください。
これを以下のコードで反復型生存時間データに加工する。
test <- tmerge(cgd0[, 1:13], cgd0, id=id, tstop=futime, infect = event(etime1), infect= event(etime2), infect = event(etime3), infect= event(etime4), infect = event(etime5), infect= event(etime6), infect = event(etime7)) test <- tmerge(test, test, id= id, enum = cumtdc(tstart))
一部の変数を抜いて表示すると下のような感じ。
> head(test[-c(2,3,6:11)]) id treat sex hos.cat futime tstart tstop infect enum 1 1 1 2 2 414 0 219 1 1 2 1 1 2 2 414 219 373 1 2 3 1 1 2 2 414 373 414 0 3 4 2 0 1 2 439 0 8 1 1 5 2 0 1 2 439 8 26 1 2 6 2 0 1 2 439 26 152 1 3
これをに対してCoxモデルを当てはめる。普通と違うところは、
- Survオブジェクトを開始時点、終了時点、イベント有無の3つを使って指定した
cluster=idを指定してロバスト分散推定が行われるようにする（GEEみたいなもの）。
fit <- coxph(Surv(tstart, tstop, infect) ~ treat, data = test, cluster = id)
結果は以下のとおり。
> summary(fit) Call: coxph(formula = Surv(tstart, tstop, infect) ~ treat, data = test, cluster = id) n= 203, number of events= 76 coef exp(coef) se(coef) robust se z Pr(>|z|) treat -1.0953 0.3344 0.2610 0.3119 -3.511 0.000446 *** ~~~~~（以下省略）~~~~~
おわりに
- 共変量やその効果が時間によって変化することを許容した解析はまた今度まとめたい。
- そんなにモフモフを着てたら暑いだろうと思って冷房をかけてあげてるのに、朝ごはんが終わると、とっとと冷房のかかってない寝室の出窓に移動するのはなぜ？
- 誤植修正しました（2022-10-05）
参考資料
*1:ここではイベント種類が1つの場合を想定して説明しているが、複数種類にも対応できるため明示した。
