2019 医薬生物学 （統計検定1級応用 過去問の解説です）
- 生存時間解析
- 仮想的コホート研究
- 感度
- ２標本t検定
2019医薬生物学① 生存時間解析の基礎
［参考サイト］ 平松雄司 様 「生存時間分析の基礎２」
https://note.com/maxwell/n/nc3ab385128d8?magazine_key=mc84c90fb0e35
［昇順に並べた生存時間の表］
イベントや打ち切りが発生した時点を昇順に並べ、各時点での生存数・イベント発生数・打ち切り数を表形式にしたもの。
time : イベントや打ち切りが発生した時点（生存時間）
n.risk : その時点での生存数（リスクセットとも呼ばれます）
n.event : イベント発生数（gehan データの場合は死亡数）
n.censor : 打ち切り数
はじめの３つを、\(t_{(j)}\:\:,\:\:n_j\:\:,\:\:d_j\:\:\:\:(j=1,2,\cdots r)\) で表す。
(1)は、生存時間データから［昇順に並べた生存時間の表］を作成する問題である。
生存時間データが、カウンティングプロセス（Counting Process : CP）形式で与えられている。
（「生存時間分析の基礎1」https://note.com/maxwell/n/nc78c55afe944 参照 ）
観測された生存時間は５回であるから、\(n=5\)
［昇順に並べた生存時間の表］は、イベントのあった時間のみを表示するので３段である。
| 時間 | \(n_j\) | \(d_j\) | 死亡確率 | \(\widehat{S}(t)\) |
| 10 | 5 | 1 | \(\dfrac{1}{5}\) | \(1\cdot \dfrac{4}{5}=\dfrac{4}{5}=0.8\) |
| 13 | 4 | 1 | \(\dfrac{1}{4}\) | \(\dfrac{4}{5}\cdot \left(1-\dfrac{1}{4}\right)=\dfrac{3}{5}=0.6\) |
| 19 | 2 | 1 | \(\dfrac{1}{2}\) | \(\dfrac{3}{5}\cdot \left(1-\dfrac{1}{2}\right)=\dfrac{3}{10}=0.3\) |
一般に、離散の場合の生存確率 \(S(i)\) は、瞬間死亡率を \(h(i)\) として、
直前の生存確率 \(S(i-1)\) に瞬間生存率 \(1-h(i)\) を乗じて得られる。
即ち、\(S(i)=S(i-1)P(X>i|X\geq i)=S(i-1)\left(1-h(i))\right)\)
このように順次条件確率を乗じていくが、上の例のように
打ち切りのない箇所では分母が約分できる。
打ち切りのある箇所では分母の約分ができなくなる。
Rによるカプランマイヤー曲線では、打ち切りの点ではグラフに ”|” が入る。
(2) ネルソン・アーレン推定は累積ハザード関数から生存関数を求める方法である。
［参考サイト］ Jin’s page 様 「Ｒと生存時間分析(1)」 https://www.cis.doshisha.ac.jp/mjin/R/36/36.html 参照
［累積ハザード関数］
連続の場合 \(\displaystyle H(t)=\int_0^th(t)dt=-\log(1-F(t))=-\log{S(t)}\)
離散の場合 \(\displaystyle H(j)=\sum_{k=1}^jh_k=-\log S(j)\) （これは、小サンプルに有効と言われている）
ネルソン・アーレン推定量 \(\displaystyle \widetilde{S}(j)=\exp\left({-\sum_{k=1}^jh_k}\right)= \prod_{k=1}^j e^{-h_k}\)
カプランマイヤー推定量 \(\displaystyle \widehat{S}(j)=S(j-1)(1-h_j)= \prod_{k=1}^j (1-h_k)\)
一般に、\(e^x\geq 1+x\) （ \(x\fallingdotseq 0\) のとき近似的にイコールが成立）なので
\(e^{-h_k}\geq 1-h_k\) \(\displaystyle \therefore \prod_{k=1}^j e^{-h_k}\geq \prod_{k=1}^j (1-h_k)\)
よって、 \(\widetilde{S}(j)\geq\widehat{S}(j)\) 題意は示されたといえる。
(3) 境界内生存時間
［基本］
［生存数の曲線下面積が平均寿命］
\(T\) は生存時間（死亡時刻・イベントが起きる時刻）を表す確率変数、 \(f(t)\) は \(T\) の確率密度関数
・終了時刻 があるとき
\(G\)は最長の生存時間で \(S(G)=0\) である。
\(S(t)=1-F(t)\:\:\:(0\leq t\leq G)\) であるので、
\(\begin{split}
E(T)&=\int _0^{+\infty}tf(t)dt=\int _0^{G}tf(t)dt\\&=\left[tF(t)\right]_0^{G}-\int_0^{G}F(t)dt=G-\int_0^{G}F(t)dt
\\&=G-\int_0^{G}(1-S(t))dt=\int_0^{G}S(t)dt
\end{split}\)
\(\begin{split}
E(T^2)&=\int _0^{+\infty}t^2f(t)dt=\int _0^{G}t^2f(t)dt\\
&=\left[t^2F(t)\right]_0^{G}-2\int_0^{G}tF(t)dt=G^2-2\int_0^{G}t(1-S(t))dt
\\&=2\int_0^{G}tS(t)dt+G^2-2\int_0^{G}tdt=2\int_0^{G}tS(t)dt
\end{split}\)
\(\displaystyle V(T)=2\int_0^{G}tS(t)dt-\left(\int_0^{G}S(t)dt\right)^2\)
第１式は生存曲線の下部の面積が平均寿命になることを意味する。
・終了時刻 G が定まらないとき
境界内平均生存時間（RMST）で平均寿命を推定する。
境界時間 \(\tau\) とすると境界内生存時間は、 \(X(\tau)=\min(T,\tau)\)
境界内平均生存時間(RMST)は、 \(E(X(\tau))=E(\min(T,\tau))\) である。
\(\begin{split}
E(X(\tau))&=\int _0^{+\infty}tf(t)dt
=\int _0^{\tau}tf(t)dt+\int _{\tau}^{+\infty}\tau f(t)dt\\
&=\left[tF(t)\right]_0^{\tau}-\int_0^{\tau}F(t)dt+\tau(1-F(\tau))\\
&=\tau F(\tau)-\int_0^{\tau}(1-S(t))dt+\tau(1-F(\tau))=\int_0^{\tau}S(t)dt
\end{split}\)
これは、\(0\leq t\leq\tau\) での \(S(t)\) の曲線下の面積を表す。
(4) 境界内生存時間の分散については、
\(\displaystyle V(X(\tau))=2\int_0^{\tau}tS(t)dt-\left(\int_0^{\tau}S(t)dt\right)^2\)
が成立する。(計算は解答参照)
(5) のとき \(S(t)=\exp(-\lambda t)\)
ハザード関数は、 \(\displaystyle h(t)=\frac{-S'(t)}{S(t)}=\dfrac{\lambda\exp(-\lambda t)}{\exp(-\lambda t)}=\lambda\) （瞬間死亡率は一定）
分布関数は、\(F(t)=1-S(t)=1-\exp(-\lambda t)\)
確率密度関数は、\(f(t)=F'(t)=\lambda \exp(-\lambda t)\)
これは指数分布であり、\(\lambda\) は単位時間のイベント発生回数である。
以下は、モデルを指数分布とした場合のパラメーターの最尤推定である。
ただし、打ち切りがあるので注意が必要である。
［打ち切りのある場合の最尤推定］
尤度関数は、イベント発生時には \(f(t)\) を乗じ、打ち切り発生時には \(S(t)\) を乗ずる。
［参考資料］福山平成大学 福井正康 様
「社会システム分析のための統合化プログラム31 -生存時間分析-」（p5）
https://www.heisei-u.ac.jp/ba/fukui/pdf/analysis31.pdf
尤度関数 \(L(\lambda)=\prod_{i=1}^n f(t_i)^{\delta_i}\cdot S(t_i)^{1-\delta_i}\) （打ち切り・イベント発生を \(\delta_i=0,1\) で表す）
観測された生存時間を \(t_1,t_2,\cdots ,t_n\) 、イベントの観察時刻を \(t_{(1)},t_{(2)},\cdots ,t_{(r)}\) とする。
\(t_{(i)}\) にイベントが発生する確率は、\(f(t_{(i)})=\lambda \exp(-\lambda t_{(i)})\)
\(t_{j}\) に打ち切りが発生する確率を、 \(S(t_j)=\exp(-\lambda t_j)\) と考えると、
尤度関数は、\(\displaystyle L=\prod_{i=1}^r \lambda \exp(-\lambda t_{(i)})
\cdot \prod_{j=1}^{n-r}\exp(-\lambda t_{j})=\lambda^r\exp\left(-\lambda\sum_{i=1}^n t_i\right)\)
対数尤度関数は、 \(\displaystyle \log L=r\log\lambda-\lambda\sum_{i=1}^n t_i\)
これを最大とする \(\lambda\) が最尤推定量である。
\(\displaystyle \frac{\partial L}{\partial \lambda}=r\left(\dfrac{1}{\lambda}\right)-\sum_{i=1}^n t_i=0\) とすると、最尤推定量は \(\widehat{\lambda}=\dfrac{r}{\sum_{i=1}^n t_i}\)
表から \(r=3\) であるから、
\(\widehat{\lambda}=\dfrac{3}{10+13+18+19+23}=0.036\)
境界内生存時間については、
\(\displaystyle
E(X(\tau))=\int_0^{\tau}S(t)dt=\int_0^{\tau}\exp(-\lambda t)dt
=\dfrac{1-\exp(-\lambda \tau)}{\lambda}\)
\(V(X(\tau))=\dfrac{1-2\lambda\tau \exp(-\lambda \tau)-\exp(-2\lambda \tau)}{\lambda}\) （計算は解答参照）
\(\lambda=0.036\:\:,\:\:\tau=20\) を代入すると最尤推定値が求まる。
\(\exp(-\widehat{\lambda} \tau)=\exp(-0.72)=\exp(0.28-1)=\dfrac{e^{0.28}}{e}=\dfrac{1.3231}{2.7183}=0.486738\) を用いると、
\(\displaystyle
\widehat{E(X(\tau))}=\dfrac{1-\exp(-\widehat{\lambda} \tau)}{\widehat{\lambda}}
=\dfrac{1-\exp(-0.72)}{0.036}=\dfrac{1-0.486738}{0.036}=\dfrac{0.573262}{0.036}=14.257\)
\(\widehat{V(X(\tau))}=\dfrac{1-2\cdot(0.72)(0.486738)-(0.486738)^2}{(0.036)^2}=\dfrac{0.062183}{(0.036)^2}=47.98\)
この境界内生存時間の期待値は、モデルを指数分布としてパラメーターを最尤推定し、そこから算出したものである。解答には、カプランマイヤー法による推定値が16.3 と求められており、最尤推定による推定値より大きくなっている。
［参考サイト］
理数アラカルト様 「指数分布に対する最尤法」
https://risalc.info/src/st-maximum-likelihood-exponential.html
大阪大学大学院基礎工学研究科 狩野 裕 様 「不完全データの分析」
http://seisan.server-shared.com/611/611-71.pdf