- パッケージとデータセットの準備
- transMat( )を使って遷移行列（transition matrix）を指定する
- msprep( )でデータを縦長に変形する
- coxph( )とmisfit( )を使ってモデルを当てはめる
- probtrans( )を使って各状態を取る確率を推定する
- 遷移確率をプロットする
- おわりに
- 参考資料
生存時間データを、イベントの種類と1人あたり許容される回数をもとに4つに分けることは述べた。
ここでは、「複数の種類のイベントを合計2回以上経験しうる状況（多状態型生存時間データ, multi-state type survival data）」について、（自分の理解のため）まとめたい。
例えば「対象者が再発後に死亡する場合もあるし、再発せずに死亡する場合もある」という状況では、イベントの種類は「再発」・「死亡」の2種類あり、人によっては「再発→死亡」のように2回イベントを経験しうる。
多状態モデル（multi-state model）では、各時点において「ある状態にとどまる or 移る確率（遷移行列） 」を推定し、それを目的の時点まで掛け合わせることで、「ある時点において各状態に収まっている確率
」を求める（Aalen-Johansen推定量）。
ここで は初期時点における各状態の割合。遷移行列
の要素は「ある状態にいる条件の下で、次の瞬間に他の状態へ移行する瞬間確率」で、遷移ハザード（transition hazard）または遷移強度（transition intensity）と呼ばれるもの。Nelson-Aalen推定量で推定される（詳細は割愛）。
パッケージとデータセットの準備
mstateパッケージはsurvivalパッケージも要求されるのでインストールしておく。
library(mstate) data(ebmt3)
European Society for Blood and Marrow Transplantation（EBMT）が1995年から1998年に収集した2204人の骨髄移植患者のデータ。 データセットの概要は以下のとおり。1人1行の横長データからスタート。
> head(ebmt3) id prtime prstat rfstime rfsstat dissub age drmatch tcd 1 1 23 1 744 0 CML >40 Gender mismatch No TCD 2 2 35 1 360 1 CML >40 No gender mismatch No TCD 3 3 26 1 135 1 CML >40 No gender mismatch No TCD 4 4 22 1 995 0 AML 20-40 No gender mismatch No TCD 5 5 29 1 422 1 AML 20-40 No gender mismatch No TCD 6 6 38 1 119 1 ALL >40 No gender mismatch No TCD
id：患者識別番号prtime：移植後、血小板回復までの時間（日）prstat：血小板回復状態（1=回復、0=未回復）。1つ目のイベント。rfstime：移植後、再発または死亡（= relapse-free survival）までの時間（日）rfsstat：無再発生存状態（1=再発または死亡、0=無再発生存）。2つ目のイベント。dissub：原疾患サブカテゴリー（AML=急性骨髄芽球性白血病、ALL=急性リンパ芽球性白血病、CML=慢性リンパ芽球性白血病）age：年齢drmatch：ドナーとレシピエントの性別がマッチしていたかtcd：T細胞除去（移植後のGVHD予防のため移植骨髄からT細胞を除去する処理のこと？）
想定する多状態モデル
骨髄移植後の血小板回復と無再発生存をイベントとして、以下のような多状態モデルを考える。
3つの状態：
- status 1：BMT直後。血小板未回復かつ無再発生存
- status 2：血小板回復
- status 3：無生存再発
状態間の移行：
- trans 1：BMT後（status 1）→ 血小板回復（status 2）
- trans 2：BMT後（status 1）→ 再発または死亡（status 3）（つまり、血小板が回復しないまま再発 or 死亡）
- trans 3：血小板回復（status 2）→ 再発または死亡（status 3）
再発または死亡（status 3）からは他の状態に移行することがない。このように「入ったら出てこない状態」をabsorbing statusと呼ぶ。
transMat( )を使って遷移行列（transition matrix）を指定する
遷移行列は、起こりうる状態間の移行を指定して、その移行に対して番号を振ったもの。
引数xにリスト形式で指定する。
リストの 番目の要素は、「status
から、どのstatusへ移行しうるか」をベクトルで与える。
下の指定方法だと、リストの1番目の要素が
c(2,3)なので、「status 1からはstatus 2とstatus 3に移行しうる」ことを指定している。
namesは分かりやすいように各statusに名前をつけただけ（Tx = treatment; PR = platelet recovery; RelDeath = relapse or death）
tmat <- transMat(x = list(c(2, 3), c(3), c()), names = c("Tx", "PR", "RelDeath"))
> tmat
to
from Tx PR RelDeath
Tx NA 1 2
PR NA NA 3
RelDeath NA NA NA
各移行には勝手に識別番号を振ってくれている。
上記のように「健康」→「病気」→「死亡」と不可逆的に推移していくモデルに対しては、trans.illdeath( )という関数が準備されている。
tmat <- trans.illdeath(c("Healthy", "Illness", "Death"))
path( )で遷移経路を確認する
あまり必要ないかもしれないが、辿りうる経路はpath( )で表示できる。
> paths(tmat) [,1] [,2] [,3] [1,] 1 NA NA [2,] 1 2 NA [3,] 1 2 3 [4,] 1 3 NA
この場合は4通りの経路がある。 すなわち、
- status 1から動かずに打ち切りを迎える
- status 1 → status 2と移行して、その後は動かずに打ち切りを迎える
- status 1 → status 2 → status 3
- status 1 → status 3
msprep( )でデータを縦長に変形する
msprep( )*1を使って、1行1患者の横長データを「1行1移行」の縦長データへ変換する。 指定する主な引数は以下のとおり。
time,status：モデルで扱う状態（status）と、その状態に至るまでの時間（time）を表している変数を指定する。初期状態は通常変数になっていないので、NAを明示的に書く。data：元のデータtrans：先程作成した遷移行列keep：変換後も保持したい変数。ここでは最終観察期間が分かるようにprtimeを残した。
d_msbmt <- msprep(time = c(NA, "prtime", "rfstime"), status = c(NA, "prstat", "rfsstat"), data = ebmt3, trans = tmat, keep = c("dissub", "age", "drmatch", "tcd", "prtime"))
変形後はこんな感じ。
> head(d_msbmt) An object of class 'msdata' Data: id from to trans Tstart Tstop time status dissub age drmatch tcd prtime 1 1 1 2 1 0 23 23 1 CML >40 Gender mismatch No TCD 23 2 1 1 3 2 0 23 23 0 CML >40 Gender mismatch No TCD 23 3 1 2 3 3 23 744 721 0 CML >40 Gender mismatch No TCD 23 4 2 1 2 1 0 35 35 1 CML >40 No gender mismatch No TCD 35 5 2 1 3 2 0 35 35 0 CML >40 No gender mismatch No TCD 35 6 2 2 3 3 35 360 325 1 CML >40 No gender mismatch No TCD 35
from：移行前のstatusto：移行先のstatustrans：移行の種類を識別する番号Tstart：観察開始 or 直前の状態から移行した時点Tstop：観察終了時点（移行した時点）status：その移行が観察されたかどうか。
keepで指定した共変量については説明省略。
例えば1行目は、「患者1ではtime=0から観察して、time=23でstatus 1からstatus 2への移行が観察された（status=1）」ことを示している。このため、「同期間でstatus 1からstatus 3への移行が確認されなかった（status=0）」ことが2行目で書かれている。
events( )で移行の様子を確認する
event( )を使うと下のように移行の様子が確認できる。
> events(d_msbmt) $Frequencies to from Tx PR RelDeath no event total entering Tx 0 1169 458 577 2204 PR 0 0 383 786 1169 RelDeath 0 0 0 841 841 $Proportions to from Tx PR RelDeath no event Tx 0.0000000 0.5303993 0.2078040 0.2617967 PR 0.0000000 0.0000000 0.3276305 0.6723695 RelDeath 0.0000000 0.0000000 0.0000000 1.0000000
$Frequenciesで頻度が集計されている。
Txからスタートした人（total entering）が2204人いる。PRに移行した人が1169人、RelDeathに移行した人が458人、何もイベント発生せず打ち切りになった人（no event）が577人。PRに移行した1169人のうち、RelDeathに移行した人が383人、打ち切りが786人。RelDeathからは移行する人がいない
$Proportionsは割合にしてくれている。説明は割愛。
coxph( )とmisfit( )を使ってモデルを当てはめる
移行の種類ごとに異なったベースラインハザードを仮定することになるので、strata(trans)と指定して層別Coxハザードモデルを当てはめる。
ここでは引数method = "breslow"としてBreslow近似法を用いることを指定した。他にEfron近似法（"efron"）も可。直接法（"exact"）は計算にすごく時間がかかってしまう（というかフリーズした）ので、やめた方が良さそう。*2
coxph( )の結果と遷移行列をmsfit( )へ渡す。引数vartypeで共分散の推定方法（"greenwood"または"aalen"）を指定できる（点推定値には影響なし。どちらもそんなに変わらない？）。
fit <- coxph(Surv(Tstart, Tstop, status) ~ strata(trans), data = d_msbmt, method = "breslow") simdatfit <- msfit(fit, trans = tmat, vartype = "greenwood")
ここで得られたsimdatfitは、以下の結果のリスト。
Haz：移行の種類別の累積ハザード推定値。time,Haz,transを変数として含んだデータフレームvarHaz：累積ハザード推定値の分散共分散。time,Haz,trans1,trans2を変数として含んだデータフレーム（trans1=trans2なら分散、異なるなら共分散）trans：使用された遷移行列
probtrans( )を使って各状態を取る確率を推定する
misfit( )で作成したオブジェクトをprobtrans( )に渡す。指定する主な引数は
direction：予測の方向？="forward"ならpredtから開始して最後まで。=fixedhorizonならpredtまで。predt：予測の開始時点または終了時点（directionの指定による）。method：分散の推定方法。misfit( )と同じくGreenwood法かAalen法を選択できる。misfit( )で選んだものに合わせておく。
simdatprob <- probtrans(simdatfit, predt = 0, method = "greenwood")
結果の中身は次のようになっている。
> str(simdatprob) List of 7 $ :'data.frame': 506 obs. of 7 variables: ..$ time : num [1:506] 0 1 3 4 6 7 8 9 10 11 ... ..$ pstate1: num [1:506] 1 0.999 0.999 0.998 0.997 ... ..$ pstate2: num [1:506] 0 0.000454 0.000908 0.000908 0.000908 ... ..$ pstate3: num [1:506] 0 0.000454 0.000454 0.000908 0.001817 ... ..$ se1 : num [1:506] 0 0.000642 0.000786 0.000907 0.001111 ... ..$ se2 : num [1:506] 0 0.000454 0.000642 0.000642 0.000642 ... ..$ se3 : num [1:506] 0 0.000454 0.000454 0.000642 0.000908 ... $ :'data.frame': 506 obs. of 7 variables: ..$ time : num [1:506] 0 1 3 4 6 7 8 9 10 11 ... ..$ pstate1: num [1:506] 0 0 0 0 0 0 0 0 0 0 ... ..$ pstate2: num [1:506] 1 1 1 1 1 1 1 1 1 1 ... ..$ pstate3: num [1:506] 0 0 0 0 0 0 0 0 0 0 ... ..$ se1 : num [1:506] 0 0 0 0 0 0 0 0 0 0 ... ..$ se2 : num [1:506] 0 0 0 0 0 0 0 0 0 0 ... ..$ se3 : num [1:506] 0 0 0 0 0 0 0 0 0 0 ... $ :'data.frame': 506 obs. of 7 variables: ..$ time : num [1:506] 0 1 3 4 6 7 8 9 10 11 ... ..$ pstate1: num [1:506] 0 0 0 0 0 0 0 0 0 0 ... ..$ pstate2: num [1:506] 0 0 0 0 0 0 0 0 0 0 ... ..$ pstate3: num [1:506] 1 1 1 1 1 1 1 1 1 1 ... ..$ se1 : num [1:506] 0 0 0 0 0 0 0 0 0 0 ... ..$ se2 : num [1:506] 0 0 0 0 0 0 0 0 0 0 ... ..$ se3 : num [1:506] 0 0 0 0 0 0 0 0 0 0 ... $ trans : int [1:3, 1:3] NA NA NA 1 NA NA 2 3 NA ..- attr(*, "dimnames")=List of 2 .. ..$ from: chr [1:3] "Tx" "PR" "RelDeath" .. ..$ to : chr [1:3] "Tx" "PR" "RelDeath" $ method : chr "greenwood" $ predt : num 0 $ direction: chr "forward" - attr(*, "class")= chr "probtrans"
リストのうち最初の3つの項目は、それぞれ、状態1,2,3からスタートして移行していく様子のデータ。この3つのリスト項目の中には、pstate1, pstate2, pstate3とそれに対応する標準誤差（se1など）が含まれている。例えば、状態1の人が状態2に移行していく遷移確率（transition probability）はsimdatprob[[1]]$pstate2に入っている。
遷移確率をプロットする
probtrans( )で作成したオブジェクトをplot( )に渡す。指定する主な引数は次のとおり。
from：どの初期状態からスタートした遷移確率を表示するか指定する。ord：積み重ねる順序を下から指定するtype：プロットの形式を指定。- "filled"：積み重ねプロット（合計が100%で一定）で、中を塗ったもの
- "stacked"：積み重ねプロットで中を塗らないもの
- "separate"：それぞれの状態について別のプロットを作成（状態の数だけプロットが作成される）
- "single"：上記の複数のプロットを1枚に重ねる
plot(simdatprob, from = 1, ord = c(1,2,3), type = "filled")
あまり色が見やすいとは言えないし、ggplotを使ったらもっとキレイに描けそう（な気がするだけ）。
おわりに
- 「生存時間データ4分類シリーズ」もようやく終わりました。
- 次は時間依存性共変量の扱いをまとめる予定です。
参考資料
- Multi-stateモデルとmstateパッケージについての論文。
- Aalen-Johnson推定について日本語で説明された記事。詳しく説明されていて参考になりました。
- オスロ大学の講義資料（lecture7_16.pdfがmstateパッケージの解説スライド）
