Jackson, C. (2019) Multi-state modelling with R: the msm pachage. Version 1.6.8.

　パネルデータに多状態マルコフモデルをあてはめるRパッケージ msm の解説。実戦投入のための儀式として読んだ。もとはJ. Statistical Software の2011年の論文だそうだ。
　おっと、いまみたら今年9月に1.6.9が出ている… なんてこった…

　いわく、
　連続時間上の時点で、ある個人がとりうる状態が何個かあるとしよう。時点\(t\)における状態を\(S(t)\)とする。状態\(r\)から\(s\)への遷移の即時リスクを遷移強度 $$ q_{rs}(t, z(t)) = \lim_{\delta t \rightarrow 0} \frac{P(S(t+\delta t) = s | S(t) = r) }{ \delta t } $$ とする。\(z(t)\)は時変する説明変数。[右辺から\(z(t)\)が抜けているけど、気持ちとしては、分子は\( P(S(t + \delta t) = s | S(t) = r, Z(t) = z(t)) \) という感じであろう]
　遷移強度の行列を\(Q\)とする。対角要素を \(q_{rr} = – \sum_{s \neq r} q_{rs} \)と定義する。各行の和は0になる。
　\(Q\)を推定したい。
　[なるほどね。連続時間マルコフモデルってまったく触れたことがなかったけど、推定すべきパラメータは一定期間における遷移確率じゃなくて遷移の即時リスクなわけだ。要は連続時間生存時間モデルで一定期間における生存率じゃなくてハザードを推定するようなものね]

　以下ではマルコフモデルについて考える。つまり、遷移強度\(q_{rs}(t, z(t))\)は\(t\)までの遷移過程の履歴から独立であるとする。
　ある状態\(r\)が維持される時間を滞在時間と呼ぶ。時間等質な連続時間マルコフモデルでは、滞在時間は率\(-q_{rr}\)で指数分布に従う。\(r\)の次の状態が\(s\)となる確率は\(-q_{rs} / q_{rr} \)である。

　msmパッケージをつくったのは疾患のモデリングに適用したかったからである。この分野ではillness-deathモデルがよく使われる。人は健康状態から疾患状態に遷移する。逆向きに遷移することもある。どちらの状態からも、吸収状態absorbing stateと呼ばれる状態へと遷移することがあり、ここに遷移したらもう戻れない。典型的には死のことである。

　観察時点は遷移時点と異なる。それは固定されていたりランダムだったりするし、informativeなこともある(自分自身によって選ばれている)。msmパッケージはinformativeな観察時点のデータを扱わない。[←なるほどね、離散時間データでの状態依存欠損と同じように、連続時間データでは状態依存な観察タイミングという問題が起きるのか。超めんどくせえな]

　遷移確率行列\(P(t)\)を考える。過程が時間等質だとして、要素\(p_{rs}(t)\)は、時点\(u\)において状態が\(r\)であるという条件の下で、時点\(t + u\)で状態が\(s\)である確率を指す。遷移がいつ起きるかについて語っているわけではない点に注意。[おっとっと。ここまでとは\(t\)の意味が変わってきたぞ。もはやカレンダー時間ではなくて観察時点間の経過時間である]
　遷移確率行列と遷移強度行列の間には $$ P(t) = Exp(tQ) $$ という関係がある。ただし、行列\(X\)の\(k\)回の積を\(X^k\)として、\(Exp(X) = 1 + X^2 /2! + X^3 / 3! + \cdots\)である。信頼性の高い計算が難しいことで有名である。[ははは]
　もっと単純なモデルで、解析的に表現すると… [事例。メモ省略]
　msmパッケージでは、5状態モデルまでは解析的に求める。6状態以上では… [難しくてわからんので省略]

　ある個人\(i\)の、\(S(t_j)\)から\(S(t_j+1)\)への状態遷移の、尤度への貢献を考えよう。

• 通常は$$ L_{i,j} = p_{S(t_j), S(t_j+1)} (t_{j+1} – t_j) $$
• 慢性疾患の観察研究では、死亡時点は正確にわかるけどその直前の状態がわからないということもある。死亡状態を\(D\)とし、\(S(t_{j+1}) = D\)としよう。この場合、\(S(t_j)\)から\(S(t_{j+1})\)のあいだにとりえた全状態\(m\)について足し上げればよい。$$ L_{i,j} = \sum_{m \neq D} p_{S(t_j), m} (t_{j+1} – t_j) q_{m,D} $$[時間等質ならば、ということであろう。なるほどね、これも観察時点が状態依存である例だが、特定の状態に遷移した瞬間に観察が生じると決まっているなら話は簡単なわけだ]
• 状態遷移の時点を常に正確に観察できるならば、$$L_{i,j} = \exp(q_{S(t_j), S(t_j)} (t_{t+1} – t_j)) q_{S(t_j)S(t_j+1)}$$ [添字が長くなりすぎてわけがわかんないですね。前の状態を\(x = S(t_j)\), 後の状態を\(y = S(t_j+1)\)、観察間の間隔を\(t = t_{t+1} – t_j\)と書けば$$L_{i,j} = \exp(q_{xx} \times t) q_{xy}$$\(\exp(q_{xx} \times t)\)は2時点目まで状態\(x\)が持続すること、\(q_{xy}\)は2時点目で遷移した行き先が\(y\)であること、を表しているのであろう]
• 観察の右側打ち切りは、打ち切り時点の状態がわかっているなら気にしなくて良い。いっぽう、たとえば患者が未知の状態で脱落しちゃうときみたいに、打ち切り時点の状態が未知ということもある。この場合は、打ち切り時点の状態\(S(t_{j+1})\)のありうる集合を\(C\)として、$$ L_{i,j} = \sum_{m \in C} p_{s(t), m} (t_{j+1}-t_j)$$
• 極端な話、すべての対象者において打ち切りがあったとして…[わけわかんなくなってきたのでパス]

まあいずれにせよ、完全尤度\(L(Q)\)はすべての\(L_{i,j}\)の積である。

　共変量について。
　遷移強度\(q_{rs}\)について次の比例ハザードモデルを考える: $$ q_{rs}(z(t)) = q_{rs}^{(0)} \exp(\beta^T_{rs} z(t)) $$ 共変量\(z(t)\)の観察時点が状態の観察時点とは異なるという場合もある。そういう場合は、共変量は観察時点間では変わらないと仮定する。[←そうするしかないのか… 他になんとかやりようがあるような気もするけどな…]

　隠れマルコフモデル(HMM)について。
　HMMとは、観察が観察されていない状態によって条件付けられている確率分布(emission分布)に従っていると考えるモデルである。HMM推定の理論はたいてい離散時間を考えており、医学ではあまり用いられない。
　たとえば誤分類モデル。対象者\(i\)の真の状態を\(S_i(t)\)、観察を\(O_i(t)\)とし、$$ r_{rs} = Pr(O(t_{ij}) = s | S(t_{ij} = r))$$を要素とする誤分類行列を\(E\)とする。ここに説明変数\(w(t)\)を入れたい場合は、多項ロジスティック回帰モデル $$ \log \frac{e_{rs}(t)}{e_{r,s_0}(t)} = \gamma^\top_{rs} w(t) $$を考える。ベースライン状態\(s_0\)はもっとも確率の高い状態にしておくと数値的に安定する。
　一般化すると…[パス]
　尤度は…[パス]
　一般HMMの例…[パス]

　…ここからはmsmパッケージの使い方。メモしてもしょうがないので省略する。わかりやすい説明であった。
　なお、入力データはある対象者のある観察が行になる。\(Q\)は非対角要素のうちどこがゼロでないかを指定できる。面白いのは、たとえばデータのうえで、隣り合う2時点で状態1から2に遷移しているようにみえるケースがあっても、\(q_{1,2}\)を0と指定できるという点。2時点のあいだで別の状態に遷移しているかもしれないからである。なるほど。
　\(Q\)の推定だけでなく、対象者ごとにもっともありそうなパスを復元したりできます。一般HMMも推定できます。などなど、いろいろ書いてあったけど、必要になったらちゃんと読む、ということで…