霞と側杖を食らう

欲望はこちらに[https://www.amazon.jp/hz/wishlist/ls/2EIEFTV4IKSIJ?ref_=wl_share]

好きなモデルは?と問われたならば

大政絢だと私は答える。こんな超絶真面目な冗談は置いておいて、モデル(数理モデル、とくに統計モデル※1)について少し書いていく。ただし、モデル観というものは、人によって違っていると思っていて、ここに書くものが一概に正しいというわけではないであろうことは注意されたい。モデルは分析対象や分析対象を扱う分野によって異なり、見えるところで同じような表現をされていても、暗黙のところに存在する見方のようなものがあるように思う。自分の知らないモデル観はきっと存在して、それについて知らないので書けない(違うモデル観を言語化して文章として公開してくださると私は嬉しく思う)。また、偉そうに色々書いているが、自分の理解もけっして深く広くあるわけではないので注意されたい。

さて、モデルを作って使っていく過程について、ざっくりと説明する。分析対象をまず決める。分析対象は複雑すぎるので(そうでなければモデルで考える必要がない)、モデルに必要な要素を残して、それ以外は捨象する。モデルが数理モデルであれば、必要な要素が仮定となって、数理構造となる。その仮定、構造に従って、結果を導いていく。その結果を分析対象に照らし合わせ、解釈して、理解を促し、必要ならば意思決定に繋げていく。

モデルというものに不慣れな人もいるかもしれないので※2、例を使って、モデルで考えてみる。とある隣町に、とある時刻までに到着したいとする。このとき、移動時間がどのくらいになるのかということが分析対象となる。最も簡単なモデルを設定するとして、隣町まで直線距離が8kmということが分かっているとする。時速4kmで歩けるとしたら、2時間で着けると解析的に導き出せる。それならば、2時間前に出発すれば良いという結論が出るが、これで良いのだろうか。最も簡単なモデルには二つの大きな要素があった。隣町まで8kmというものと時速4kmで歩けるというものだ。直線距離は8kmかもしれないが、歩ける道があるのは8kmより長いかもしれないし、信号や踏切、坂道、天候、疲労の関係で時速4kmでは歩けないかもしれない。2時間ぴったりで着けなくても大体の時間が分かればいいというのならば、先の最も簡単なモデルによる結論で十分なのだが、遅れると困る場合や、早く着きすぎるのも嫌だという場合は、モデルの要素をもう少し複雑にして分析する必要がある。複雑化した結果、解析的に解けなければ、数値シミュレーションするなどして、結果を導くことになる。疲れているならば2時間20分、雨の場合は2時間10分になるとか導けて、意思決定ができる。

簡単にモデルを作ってみたが、持っている情報や知識、分析の目的で、モデルの構築の仕方は変わってくるだろう。データを持っているならば統計モデルを組んで、統計的推測をして確率論的に評価することもできる。隣町の例で言えば、2時間5分だとか1時間55分だとか、かかった時間の記録があれば、確率分布を推測して考えることもできる。付随して、天候だとか時間帯だとか曜日とか疲弊具合のデータがあれば、それを用いた構造のモデルを組むことも可能だ。かかる時間の予測分布を出せば、リスク評価して、自分の損得に応じた出発時刻の決定ができる。

今では統計ソフトも充実しているので、PCでポチポチ触れば、何かしら結果を出して、それに基づいた結論が出せるわけだが、中身が分からなくとも、ただ手順をたどっているだけの人もいる。モデルの目的が予測であって、それで利益が出るからいいのだという人もいるだろうが、モデルの目的が説明である場合、果たしてそれでいいのかと思うことがある。少しでも理解を促すために、以下では統計モデルについて少し書いていくことにする。

統計モデリングの一連の手続きの中には、モデルと推定手法と計算アルゴリズムがある。モデルでいえば、GLM(Generalized Linear Model)なら、指数分布族の確率分布とリンク関数を指定して数理構造を表す※3。選んだ数理構造に応じて、表現できる範囲が決まる。この数理構造の範囲の中から、何かしらの指標を最適化したり、揃えたりすることで、パラメータを推定して、モデルの構造の中から特定する。推定手法は最尤法であれば、尤度関数を最大化してモデルのパラメータを推定する。推定手法が解析的な解が得られるのならばそれでいいが、そうでない場合はコンピュータに計算してもらう必要がある。最尤法でよく使われる最適化の方法はBFGS法だとか、Nelder-Mead法だとかがある。これらのアルゴリズムに従って、最適解を得て、モデルを決定する。

これらの理解がごちゃごちゃの人がいる。たとえば、線形モデルのことを推定手法である最小二乗法だと言っている人はたまにいる。モデルを選ぶと、よく使う手法やアルゴリズムがセットになって理解するから一緒のものだと思ってしまいやすいのかもしれないが、OLS(Ordinary Least Squares)とGLM(Generalized Linear Model)が同カテゴリーのものとして並べられているのを見ると、違和を感じてしまう。

雑然とした理解にしないようにするには、モデルと推定手法と計算アルゴリズムのそれぞれの理論や体系を知るべきなのだと思う。モデルについては、仮定があって数理に従って結果が出るので、仮定が妥当でなければ適切な結果が得られるとは限らないという話は冒頭の方で書いた(ただし、統計モデルの場合、結果が出るまでの数理に推定手法が入っていて、モデルは表現できる範囲なので、注意が必要かもしれない)。推定手法については、手法はいくつもある(おぼろげに数字が浮かんでくる、おぼろげ推定法だって、推定法のひとつだ)ということが大前提で、ということは、推定手法の良し悪しは比較したくなるわけで、推定する上で望ましい性質(不偏性や一致性、十分性、精度などなど)で以て、推定手法を評価して、良い推定手法が使われているということを知っておくべきかもしれない。そして、そういった望ましい性質が成立するためには条件があるので、その条件は何かを知っておくと良いだろう(回帰分析のガウス=マルコフ定理あたりを丁寧に追うと、何が条件で導かれているかの簡単な例になるだろう※4)。数値計算アルゴリズムについては、まず単純なアルゴリズムでいいから、解析解が分かっている問題を自分で実装してみるのがよいかもしれない。許容誤差の話とか、オーバーフローの話とか、実行速度の話とか、上手くいったり上手くいかなかったりと経験できる。私は全然詳しくないので、このことについては詳しい人に聞いてほしい※5。

諸々の仮定だとか条件が妥当なのかどうかは、ドメイン知識から判断するしかない場合が多い。得られた結果が妥当でないとか、測定やサンプリングの問題だとか、その分野の知識から仮定に反していることと考えられるだとかから判断する※6。仮定や条件がズレていたとして、そのズレが本来の目的にとってどのくらい問題になるのかだとか、現状の制約(経済的理由だとか倫理的理由だとか技術的理由だとか)から実行可能かどうかだとか、そういったことも気にしないといけない。

色々とつらつらと書いてきたが、結局何が書きたかったんだ?と問われたならば、こういうことだと私は答える。モデルを構築するときの仮定、数理モデルの数理の性質、実装があるならその実装の性質、ドメイン知識などなどが分析、推論、意思決定に影響する。これらのすべてを完全に理解するのは極めて困難だが、そういうものの上に成り立っていることを知らないのでは、認識できていない危うさを抱えてしまう。「予測が外れた、ああ無情」で終わるなら、その危うさは無視してもいいかもしれないが、モデル分析、とくに統計分析は簡単に使えてしまえて、重要な意思決定につながりうる。ときとして、人の命を奪ってしまう判断にもつながってしまうことがある。公衆衛生や経済政策・金融政策の決定といった多大な影響を与えるものは当然だが、もっと身近な決定も命につながりうる。言ってしまえば、統計モデルは便利だが危うさもある刃物なのだ。その刃物の特性を事前に把握しておくためには、完全な理解はできないなりにも、知らないことを知ろうとする努力を続けるべきだろうし、果てしなく広大な関連知識が背景にあるということを知っておいてほしいということが書きたかった※7。完璧なモデルなどといったものは存在しないが、モデルで考えることは有益となりうる。そのためには、モデルで考えることへの理解を深く広くしていくしかないのだ。

 

※1) ワイスバーグ『科学とモデル シミュレーションの哲学入門』では、モデルを構造と解釈の組み合わせから成り立つとし、モデルの構造を具象モデル、数理モデル数値計算モデルの3種類に分類している。作者がシミュレーション寄りの人ではあるが、モデルとは何かを理解する一歩目になる。ペイジ『多モデル思考:データを知恵に変える24の数理モデル』では、数理モデルの用途をREDCAPEと7つに分類している(Reason, Explain, Design, Communication, Act, Predict, Explore)。

※2) 浜田『その問題、数理モデルが解決します』は、社会科学や統計学のモデルを平易な言葉で説明していて、これまで数理モデルというものに触れてこなかった人が、その考え方を身に着けるのに役立つと思う。ほかに、江崎『データ分析のための数理モデル入門 本質をとらえた分析のために』あたりも参考になる。アヒル本こと、松浦『StanとRでベイズ統計モデリング』は、統計モデルについて、丁寧に言語化されていて理解を促進してくれる。ベイズ統計については色々な話があるみたいだが、自分で分かっていないことが多いので、ここでは書かない。

※3) Dobson『一般化線形モデル入門』あたりが詳しい。LMMに進みたかったら、McCulloch『線形モデルとその拡張』あたりが良いかもしれない。

※4) こういった話は浅野・中村『計量経済学』などの、数理を解説している計量経済学のテキストを読むと分かるようになると思う。数式を介さないで理解に到達するのは極めて難しいと思われる。ちなみに、Hansenが"A Modern Gauss-Markov Theorem"という論文でガウス=マルコフ定理を拡張しているという話が昨年か今年に話題になったが、私は読んでいないので詳細は分からない。

※5) 数値計算の話を知りたいというときに、最初に紹介してもらったのが、東大の公開している講義動画(https://ocwx.ocw.u-tokyo.ac.jp/course_11404/)だったが、松尾先生の「僕らは脆弱な数の体系の上に暮らしている」という言葉が印象的だったし、講義内容も分かりやすくて良かった。もっと奥深くに行くならば、コンピュータの中身へ掘り下げていく必要があるように思う。開発環境だとかバージョン管理とかも、私は全く詳しくないが重要だと思う。

※6) 江崎『分析者のためのデータ解釈学入門 データの本質をとらえる技術』が、こういったことに関して散在している知識を一冊にまとめていて、参考になるかもしれない。サンプリングについては社会調査法など、測定については心理学や工学などの分野を漁るといいかもしれない。

※7) 何も考えずにテキトーに予測して上手くいったという人と、問題が生じないか熟慮してモデリングして予測して上手くいっている人の評価において、予測して利益を出すことが至上とされるならば、前者の方がスピーディーで正しいのかもしれないが、前者は後者に比べて大きなリスクを抱えている可能性があることを評価者は知っていてほしい。知らないこと、思いつかないことに対して、我々は意を注ぐことはできないのだ。

 

生存時間解析の一般化ガンマ分布のモデルの推定の学習記録

【学習動機】

一般化ガンマ分布について、前回色々調べたが、

一般化ガンマ分布の数式の学習記録 - 霞と側杖を食らう

実際に、生存時間解析の、一般化ガンマ分布のパラメトリックモデルで推定すると、上手く推定できないケースがちょこちょこある。一般化ガンマ分布でtime-to-event dataの人工データを生成して、それを一般化ガンマ分布のモデルで推定することにする。

【学習内容】

パッケージの読み込みと乱数の種の設定。

library(tidyverse)
library(flexsurv)
set.seed(2021)

シミュレーションのコードはこのように。100回シミュレーション。サンプルサイズ100でデータ生成。\((\mu,\sigma,Q)\)はテキトーに乱数生成。一般化ガンマ分布でモデル指定して、BFGS法で推定。

nsim <- 100    # シミュレーション回数の設定
# 推定値を保存する
m_par_true <- matrix(NA_real_, nrow=nsim, ncol=3)
m_par_est_BFGS <- matrix(NA_real_, nrow=nsim, ncol=3)
m_est_opt <- matrix(NA_real_, nrow=nsim, ncol=2)

for(i in 1:nsim){
  n <- 100    # サンプルサイズ
  # 真のパラメータ生成
  mu <- rnorm(1,mean = 0, sd=10)
  sigma <- rgamma(1, shape = 2, rate = 1)
  Q <- rnorm(1, mean = 0, sd=2)
  #Q <- 4 # Qが4あたりより大きいとn=100じゃ推定が厳しいように見える(後ほど実行)
  #Q <- -4 # Qが負値がある程度以上だとBFGS法では収束しないことが増えるように見える(後ほど実行)

  par_true <- c(mu, sigma, Q)
  names(par_true) <- c("mu", "sigma", "Q")
  
  # データ生成
  df <- data.frame(
    time = rgengamma(n = n, mu = mu, sigma = sigma, Q = Q),
    event = sample(c(0,1), size = n, replace = TRUE, prob = c(0.2, 0.8))  # 2割がランダムで打ち切り
  )
  m_par_true[i,] <- par_true
  
  tryCatch(
    {
      fit_BFGS <- flexsurvreg(Surv(time, event)~1, data = df,
                              dist = "gengamma", method = "BFGS")
      m_par_est_BFGS[i,] <- fit_BFGS$res[,1]
    },
    error=function(e){
      m_par_est_BFGS[i,] <- numeric(3)
    },
    warning=function(e){
      m_par_est_BFGS[i,] <- numeric(3)
    },
    finally = {},
    silent = TRUE
  )
}

df_res <- 
  cbind(m_par_true, m_par_est_BFGS) %>% 
    as.data.frame() %>% 
    set_names(c("mu_true","sigma_true","Q_true",
                "mu_est","sigma_est","Q_est"))

df_res2 <- df_res %>% 
  mutate(
    mu_diff = abs(mu_true-mu_est),
    sigma_diff = abs(sigma_true-sigma_est),
    Q_diff = abs(Q_true-Q_est)
  )

結果を可視化してどうなっているか確認。赤色は真値を表し、青色は推定値を表す。青色がない部分は推定ができていない部分。赤色と青色が近くにあればその回の推定は上手くいっている一方で、離れていれば上手くいっていない。軸タイトルはちゃんと変更していない。

p_mu <- 
ggplot(df_res2) +
  theme_bw() +
  geom_line(aes(x=1:nsim,y=mu_true),colour="red") + 
  geom_line(aes(x=1:nsim,y=mu_est),colour="blue")

p_sigma <- 
ggplot(df_res2) +
  theme_bw() +
  geom_line(aes(x=1:nsim,y=sigma_true),colour="red") + 
  geom_line(aes(x=1:nsim,y=sigma_est),colour="blue")

p_Q <-
ggplot(df_res2) +
  theme_bw() +
  geom_line(aes(x=1:nsim,y=Q_true),colour="red") + 
  geom_line(aes(x=1:nsim,y=Q_est),colour="blue")

gridExtra::grid.arrange(p_mu, p_sigma, p_Q, ncol = 1)

f:id:moratoriamuo271:20210405224737p:plain

Qが真値と推定値が大きく離れているケース確認。Qが大きいと起こりやすい?

# Qが真値と推定値が大きく離れているケース確認
df_res2 %>% filter(Q_diff>1)
##      mu_true sigma_true    Q_true     mu_est  sigma_est     Q_est    mu_diff
## 1  -7.256602  2.3907504  3.309848  -5.269766 0.22579896 33.316234 1.98683683
## 2  11.170970  3.5850739  5.138105  13.496568 0.40728452 32.798580 2.32559811
## 3   5.498213  2.6399826  4.675271   6.547039 0.30846360 32.267718 1.04882608
## 4 -14.714114  0.6803831 -4.738894 -14.805919 0.49633897 -7.011095 0.09180474
## 5   2.708861  1.5552280  2.783587   3.619433 1.14671749  3.842243 0.91057210
## 6  -3.969919  1.2063532  2.131750  -2.580868 0.05956122 46.632574 1.38905148
## 7 -19.410441  1.6320494  1.888442 -18.352138 0.79317497  3.818769 1.05830329
##   sigma_diff    Q_diff
## 1  2.1649514 30.006386
## 2  3.1777894 27.660475
## 3  2.3315190 27.592447
## 4  0.1840442  2.272201
## 5  0.4085105  1.058656
## 6  1.1467919 44.500824
## 7  0.8388745  1.930327

Qが推定できていないケース確認。Qが負だと起きやすい?

# Qが推定できていないケース確認
df_res2 %>% filter(is.na(Q_diff))
##       mu_true sigma_true    Q_true mu_est sigma_est Q_est mu_diff sigma_diff
## 1 10.75032179  0.9671487 -3.786987     NA        NA    NA      NA         NA
## 2 -0.05690626  2.2481857 -3.910603     NA        NA    NA      NA         NA
##   Q_diff
## 1     NA
## 2     NA

\(Q=4\)で同じことをやって、結果確認。青線はあるが、大きく外れていることが見て取れる。

f:id:moratoriamuo271:20210405224950p:plain



サンプルサイズを増やせば解決するかもしれないので、\(Q=4\)のまま、サンプルサイズを100から10000にしてみる。一見、離れているように見えるが、縦軸をよく見ると、真値の4付近に来ている。Qがある程度以上大きいときは、サンプルサイズが大きくないと、真値付近で尤度関数が膨れないということなのかもしれない。尤度関数の形状によっては、そういうことが起きうるんだろうなと考えられる。これは推察に過ぎないが、Qが大きいときというのは、ガンマ分布の形状母数が1より小さい場合に近いかもしれない(ガンマ分布について言うと、逆に形状母数が大きい時は、つまり指数分布をたくさん足し合わせていることになるので、正規分布に近づいく。)。生存時間が高確率でかなり短く終わる分布だから上手く推定できないのかもしれない。

f:id:moratoriamuo271:20210405225134p:plain

 

\(Q=-4\)で同じことをやって、結果確認。青線が途切れることが多いのが見て取れる。

f:id:moratoriamuo271:20210405225221p:plain

 

同様に、サンプルサイズを増やせば解決するかもしれないので、\(Q=-4\)のまま、サンプルサイズを100から10000にしてみる。推定できているときは\(Q=4\)のときと同じで、近いところに来ている。しかし、推定がそもそもできていないのは多い。尤度関数が変な形をしているのかもしれない。どうすべきか分からない。

f:id:moratoriamuo271:20210405225300p:plain



【学習予定】

生存時間解析のパラメトリックモデルのざっくりしたまとめを書こうとしてまだやっていないからそのうち。

 

一般化ガンマ分布の数式の学習記録

  

【学習動機】

生存時間解析のパラメトリックモデルを見ていると、一般化ガンマ分布というものが現れるのだが、生存時間解析の和書や日本語資料を読んでもよく分からない。式の形も分かりにくいし、ソフトウェアを見てると2種類の確率密度関数があって、何が何か分からない。今回は、数式関連で英語で調べて分かったことをまとめておく。次回はシミュレーションして推定の様子を見る。

【学習内容】

対数正規分布とワイブル分布を包含する一般化として、Stacy(1962)でGeneralized Gamma Distribution(GGD)が出てきた。しかし、最尤法が上手く収束しない問題があった。とくに、対数正規分布に近いパラメータのところ(\(k \rightarrow \infty\)あたり)で上手くいかなかった。

GGDに従う確率変数の対数の分布を調べたのがPrentice(1974)で、こちら側からパラメータ化すると、推定も安定したようだ。こういう経緯があるため、こちらのパラメータ化したものGGDはGeneralized log-gamma distributionとも呼ばれる。

Rの{flexsurv}パッケージで言うと、Stacy(1962)のパラメータ化したGGDがgengamma.origに該当して、Prentice(1974)のものがgengammaに該当する。後者の方が、より一般化されていて、広い分布がとれる({flexsurv}のgengammaの\(Q<0\)の部分。後で説明する式で言う\(\lambda\))。

Stacy(1962)の定式化では、確率変数\(T\)が一般化ガンマ分布に従う(\(T\)~GG( \(\alpha\) , τ , \(k\) ))とき、確率密度関数は、

f:id:moratoriamuo271:20210405224218p:plain

と表される。一般化ガンマ分布は特殊形として、指数分布、ワイブル分布、ガンマ分布、対数正規分布を含む。

  • \(k= τ =1のときは指数分布\)

\[ f = \frac{1}{\alpha}(\frac{t}{\alpha})^0 exp[-(\frac{t}{\alpha})] =\frac{1}{\alpha} exp[-(\frac{t}{\alpha})] \]

  • \(k=1のときはワイブル分布\)

\[ f = \frac{τ}{\alpha}(\frac{t}{\alpha})^{τ -1} exp[-(\frac{t}{\alpha})^τ] \]

  • τ =1のときはガンマ分布

f:id:moratoriamuo271:20210405224327p:plain



  • \(k \rightarrow \infty\)のときは対数正規分布となる

この式は未確認。

パラメータを\((\alpha,τ,k)\)を以下のように、\((μ,σ,\lambda)\)と書き換える。

\[ μ = log(\alpha) + \frac{1}{τ} log(k) \] \[ σ = \frac{1}{τ \sqrt{k}} \] \[ \lambda = \frac{1}{\sqrt{k}} \]

として、\(T\)の対数をとって変数変換やらなんやらで操作して、元に戻すとPrentice(1974)の定式化した確率密度関数が得られる。この\(\lambda\)はRの{flexsurv}パッケージの\(Q\)に該当する。

\[ f(t;\lambda,σ,μ)= \left\{ \begin{array}{ll} \frac{|\lambda|}{σ t \Gamma(\lambda^{-2})} (\lambda^{-2})^{\lambda^{-2}} exp[\lambda^{-2}[(\frac{log(t)-μ}{σ})\lambda - exp[(\frac{log(t)-μ}{σ})\lambda]]] & (\lambda \neq 0) \\ \frac{1}{\sqrt{2\pi}σ t}exp[-\frac{(log(t)-μ)^2}{2σ^2}] & (\lambda=0) \end{array} \right. \]

指数分布となるのは\(\lambda=σ=1\)のとき、ワイブル分布となるのは\(\lambda=1\)のとき、ガンマ分布となるのは\(\lambda=σ\)のとき、対数正規分布となるのは\(\lambda=0\)のとき。

【参考文献】

【学習予定】

次回はシミュレーションして推定の様子を見る。

生存時間解析の一般化ガンマ分布のモデルの推定の学習記録 - 霞と側杖を食らう

 

生存時間解析(生存関数とハザード関数)の学習記録

 

【学習動機】

つい先日, 生まれてから1万日が経過した. そんなわけなので, 記念に生存時間解析について書く(生存時間解析に関する日本語の良い感じの記事少ないから書く).

生存時間解析をやるとなると, モデルの種類が, ノンパラメトリック(カプランマイヤー法やログランク検定の話), セミパラメトリック(Cox比例ハザードモデルの話), パラメトリック(ワイブル分布や比例ハザードモデル, 加速故障モデルの話)と出てくるが, ここではパラメトリックモデルについてのみ扱うこととする.

パラメトリックモデルに絞っても分量が多いので, 何回かに分けて書く.

今回は, 生存時間の分布関数と生存関数とハザード関数の関係について扱う.

【学習内容】

今回は, パラメトリックモデルに限らず, 生存時間解析一般について共通する話を書く.

生存時間解析では, 起点(time origin)からイベント(死亡, 失敗, 故障, 再発, 寛解など)の発生までの時間のデータ(time to event data)をモデリングする. 以降, 言葉の統一のため, とくに何か書かない場合は, イベント発生までの時間については生存時間と書き, イベントについては死亡と書き, 各データの対象は患者として書くこととする. 生存時間解析という風に, 通常とは別に区切られて扱われているのは, この生存時間のデータは非対称であることと打ち切りが存在することが特徴として挙げられるためである.

さて, ここから生存時間の分布関数と生存関数とハザード関数の関係について書く. この関係性は, 数式展開で度々登場するので, 理解しておく必要がある.

生存時間\(T\)の分布関数(distribution function)は\(t\)以下で死亡する確率であり, \[ F(t)=P(T<t)=\int_0^t {f(u)}du \] 生存関数\(S(t)\)は生存時間が\(t\)以上になる確率であり, \[ S(t)=P(T \geq t)=1-F(t) \] ハザード関数は, 患者が時点\(t\)まで生存したという条件のもとでその時間に死亡する確率である. ハザード率(hazard rate), 瞬間死亡率(instantaneous death rate), 強度率(intensity rate)や死力(force of mortality)などと呼ばれることもある.

\[ h(t) = \lim_{\delta \to 0} {\frac{P(t \leq T < t+\delta | T \geq t)}{\delta}} \] ハザード関数\(h(t)\)は生存時間の確率密度関数\(f(t)\)と生存関数\(S(t)\)で表せる.

f:id:moratoriamuo271:20210119190045p:plain

 また, \(h(t)=-\frac{d}{dt} log S(t)\)と表せる.

これは

f:id:moratoriamuo271:20210119190218p:plain

となるため. したがって, \(H(t)=\int_0^t{h(u)}du\)とすると(累積ハザード関数(cumulative hazard function)という), \[ S(t)=exp(-H(t)) \] と表せる.

以上見てきたように, 生存時間の分布は, 生存関数, ハザード関数, 累積ハザード関数のどれでも表すことができる. どれか一つを定めれば残りも決まる, 数学的に等価な関係となっている. ただし, 分析の面では等価でない. 分析のしやすさや誤差の出やすさが異なってくる.

 

比例ハザードと生存関数の関係もここで書いておく.

任意の時点でもハザードが\(a\)倍になるハザードに対応する生存関数は, \(h_a(t)=a*h(t)\) とハザード関数を表すとすると,

f:id:moratoriamuo271:20210119190315p:plain

となり, ハザードの\(a\)倍は生存関数が\(a\)乗となって現れる.

【学習予定】

次回, よく使われるパラメトリックな確率分布について書く.

次々回は, 比例ハザードモデル(PH)と加速故障時間モデル(AFT)について書くつもり.

【参考文献】

  • Collett本

www.amazon.co.jp

基本的にはこれを中心に学習している。数式が追いやすくて分かりやすいと個人的には感じている。

 

  • 大橋本

www.amazon.co.jp

 SASでの解析方法とともに解説している本。文章部分に結構、理解を深めてくれるところが多く、併読している。別に出ている応用編もあるが、あっちはSASの使い方がメインなのと、もう少し難しいことも解説されている。

  • クラインバウム本

www.amazon.co.jp

 デカい。スライドとその説明という形で解説される。重くてあんまり読んでいなかったが、デカいだけあって詳しく書いてあることも多いので、これも読まねば。

Maxwellさんがnoteで生存時間解析の解説をしているが、この解説はクラインバウム本に準拠しているようだ。

note.com

 

  • ホスマー本

www.amazon.co.jp

 個人的には、あんまり読み心地が良くないからほとんど読めていないが、この本もよく使われているみたい。

 

 

2020年の後日譚々

 2021年になったので、2020年を振り返る。どうやったって目を背けることができないのが、新型コロナであった。行く予定だったライブも中止となり、気に入っていた居酒屋も閉店に追い込まれ、友人たちと会う機会も減少し、東京五輪も延期となった。日常がどんどん侵食されている気分だった。一方で、在宅勤務ができるようになったりと、従来の生活様式で、不便だが継続されていたものが見直され、変化してきたのは良いことではあった。今後、どのように世界が変わっていくのかは見通しが全く立たない。

 この状況に対する感覚で少し似ていたのは、2011年の東日本大震災のときの感覚だった。当時高校二年生から高校三年生になる年だった私は、翌年の大学受験を控えていた。震災後数か月は自分の受験勉強に集中していたつもりだったが、秋ごろになって、今後どうなっていくのかよく分からない不安からか、少し鬱気味になってしまった。そのときは二ヶ月ほど、苦しい時期を過ごしてなんとか乗り越えることができた。2020年に感じたものは、そのときの、不安感に少し似ていた。過去の経験があったからこそ、今年は鬱な状態に深くは陥らずに済んだ(とはいえ、夏頃は少し色々なモチベーションは落ちていた)。

 昨年から、自分の滞った感性や価値観を見つめ直すと決めていた。その取っ掛かりの一つとして、映画を最低50本は見るということにしていた。実際に、51本は見ることができた(51本の中では『ミッドナイトスワン』が一番良かった)。他にも漫画や小説など、人にオススメされたものを多く読んだし、音楽も昔より幅広く聞くようにした。来年も最低映画50本、できれば100本見ていきたいし、同様に感性を鍛えるようにしていきたい。

 勉強面に関しては、医療統計・生物統計、因果推論あたりをターゲットとして深く学んでいきたいとしていた。実際、医療統計分野のいくつかの手法を学ぶことができたし、分析する前の準備諸々を経験することができた(詳しくはここでは書けないが)。因果推論は、What Ifを読みたかったが、その余裕はなかったので、2021年に読めたらいいなと思っている(勉強会が開かれているのは知っていて参加したかったが、その時期大変だった)。対外的な発表をする予定ではあったが、オンラインの発表がまだ苦手意識があってできなかった。発表用のスライドは2本作ったので、そのリンクを貼っておく。

 

医学論文を探そう ~{easyPubMed}パッケージの紹介~ / Introduction to easyPubMed - Speaker Deck
確率分布の紹介 - Speaker Deck


周辺知識として、python機械学習、英語もやっておきたかったが、あまりできなかった。機械学習は今後近いうちに使うことは少なそうなので、python、英語、あとはガウス過程あたりを少しずつ継続して学習しようと思う。

 これからどうなっていくのか見えないものが多いが、先行き見えぬときこそ、着実に堅実に自分の力を蓄えて、来るべきチャンスに備えること。そして、身体的・精神的に不安定な時期となりうるが、健やかに過ごしていくこと。この二つを軸に2021年を過ごしていこう。

喫茶店ルノアールでつかまえて

こうして文章を書き始めるとなると、結婚の予定はいつだとか、自分の健康について気にしていることはあるかだとか、年収はいくらだとか、そういうしょうもないことを知りたがるかもしれない。でもはっきり言ってね、その手の話を書く気になれないんだよ。一度ウケたもんだからって繰り返しがちなワンパターンなネタだとか、誰でも分かる下世話な話や人の悪口だとか、そういうことを書きたいわけでもないんだ、けっしてね。今から書くのは、ここ1、2年間に起こった、マジな話、気に入らないことをやっている奴らのお話だ。

まず最初に、奴らに出会ったのは、趣味でやっているバドミントンのサークルだ。社会に出てからバドミントンを続けるには、どこか社会人サークルに入らなくてはならなかった。サークルがメンバー募集している掲示板をネットで見つけて、いくつか回ってみた。その中の一つに、奴らはいた。何度かそのサークルに遊びに行ったところ、バドミントンのレベル的にも、会話の内容的にも、ある程度、自分と合っていた。サークルの代表の人と仲良くなってきたところで、奴らの一人は、話を持ち掛けてきたんだ。「今度やるイベントに、めちゃくちゃすごい人が来るのだけど、君も来ないか?」と。あ、こいつらは黒だ。真っ黒なんだと、このとき確信した。このとき確信したというのは、それ以前にグレーだと思う点がいくつかあったからだ。お酒を飲んでは、やたら夢を語りたがるし、何か夢はないのかと聞いて、語らせたがる。夢の話も、不労所得がどうだかとかって話をしてくる。この時点でグレーで判断を止めていたのは、これまで真っ黒な奴に実際に会ったことがなかったからだ。仲良くなってから、わけのわからないイベントに誘って、人のお金をむしり取ろうとしてくる行為、その魂胆に酷く失望し、少しの間、人間不信になった。

また別の機会でも、奴らに遭遇した。バドミントン以外の愉しみの一つで、ボードゲームを好んでいて、知り合いの知り合い経由で、5、6人でボードゲームをやりながらお酒を飲もうということで、一度集まった。場所は、新宿の歌舞伎町あたりの居酒屋だった。このときは、ふつうに、ふつうのボードゲームをやってお終いだった。「ふつうに」と言っても、居酒屋は値段より多少不釣り合いなサービスだった(歌舞伎町周りだと、もっと酷いところがたくさんあるから、まだマシとその場では納得してはみたが、おそらく幹事の奴は店からキックバックを受け取っていると思う)。まぁ、知らないボードゲームを楽しめたので、次に呼ばれたときも参加してみた。東京タワー付近のレンタルスペースを借りたので、そこでやろうとのことだった。とくに中身も聞かずに現地に行ってみたら、3、40人くらいいただろうか(今になって思えば、集めた数に応じてキックバックが入る仕組みなのだろう)。想像していたより人数は多いが、テーブルに分かれて、それぞれのテーブルでゲームをやることになっていた。さて、どんなゲームをやるんだろうかと見ていたら、どうやら、対戦ゲームや協力ゲームといった類じゃないゲームだった。ふつうのボードゲームは勝ち負けがある。ふつうなら、対戦している人との戦いだったり、協力してゲームマスターを倒すなり、共通課題をクリアするかのようなものである。この場で始まったのは、モノポリーに近いが、勝ちや負けの存在しないゲームで、キャッシュフローゲームというものだった。簡単にいうと、すごろくのようにマスを進めていって、みんなが不労所得を得て終わりというものだった。ルール説明をした奴らの一人は、一見、ノリが良いように見えるが、どこかコミュニケーションに、ぎこちなさが感じられる奴だった。どうも無理をしているというか、不自然というか。ゲームの内容も、無駄に夢を語らせたがる。やれやれ、真っ黒だ。こいつらも奴らなんだと、ガッカリした。ゲーム名をググってみても、同様の報告が見受けられた。気持ち悪くなって、途中で抜け出して、近くのお店でトランプ大統領が食べたとかいうハンバーガーを食べて帰った。

奴らは、夢をやたら語らせたがるし、コミュニケーションがどこかおかしい。後ろめたさがあるのだろうか。バドミントンのサークルにいた奴は普段見かけないエナジードリンクを飲んでいた。奴らの仲間から買ったのだろうか。目元や顔のハリが、連日あまり眠れてない人のそれにとても似ていた。これはボードゲームにいた奴の一人もそんな感じがしていた。おかしさの雰囲気はだんだん掴めてきた。不思議に思うのは、奴らは人を集めるのが上手い気がするってことさ。人の集め方に関するマニュアルでも出回っているのだろうか。人を損させて自分が得する仕組みを作ろうとする奴らの、ふざけたシステムがぶっ壊れしまえばいいのに、まじな話でね。

僕の話はこれでおしまい。これを書く少し前、とある喫茶店にいたんだけど、近くのテーブルで、悪い大人と悪い友人の紹介のせいでFXの情報教材に49万円を捨ててしまう人がいた。ハメられているのは、大学生くらいだろうか。自分の会計をした後、そいつらのテーブルの近くで大きな溜め息を吐いて、外に出た。やれやれ。これを読んでいる君たちは、こんないかれトンチキな奴らに出会わないことを願うし、もし会ってしまったら、距離をとってほしいし、できることなら、必勝法ってやつを思いつくものなら、そういうい仕組みを滅ぼしてほしいんだ。それが僕の夢ってやつなのさ。

 

 

 

自動手記人形サービス

劇場版ヴァイオレット・エヴァーガーデンを映画館で観賞してきた。綺麗な画と感動的な物語で、たしかに見る価値はあった。しかし、多くの知り合いが打たれている感動の大きさほどのものは感じることができなかった。その理由について、考え付いたことを少し記すことにする。

本作は、元々テレビシリーズのアニメであった。テレビシリーズは、主人公ヴァイオレットが、手紙代筆の依頼主との関わりの中で、孤児かつ戦争の道具として生まれ育ったゆえに、それまで持ち合わせていなかった思考や感情、人の心に気付きながら、人々に評価される手紙を書けるようになっていくというものであったと思う。一方で、劇場版では、お客様のために手紙を書くシーンももちろんあるのだが、メインはヴァイオレット自身の物語であると感じた。愛する少佐との再会を望むヴァイオレットの話だったのだが、それはどうも、テレビシリーズ後の物語を整った形で終わらせるための、物語であるように見えた。少佐の生存を描くことで、ヴァイオレットというキャラクターを幸せにするために作られた物語だと思われたのだ。王道で、たしかに力強いのだけれど、自分が、ヴァイオレット・エヴァーガーデンというタイトルに求めていたものは、テレビシリーズで見ていた、手紙を代筆することを通じたヴァイオレットの自己への気付きだったのだ。他者とのやり取りを通じて自己を見るという話は、先日書いた記事で少し触れた(https://moratoriamuo.hatenablog.com/entry/2020/09/22/135107)。そういうところに良さを見出していたのだろうと、劇場版を見て、テレビシリーズと比較して気付いた。映画冒頭の家のシーンを見て、この家はもしかしてあの家なんじゃないかとすぐ気付いたくらいには、テレビシリーズが結構印象に残っていたんだろう。強い印象に残るくらいのアニメだったため、観賞前の期待がかなり高まっていたというのも、良くなかったのかもしれない。

劇場版では、電信が普及し始めていた時代設定になっていた。産業革命期の時代の潮流の大きな変化が起きている、あのような世界観は個人的に、とても好みなのだとも再認識した。プレステ2のゲームで『ポンコツ浪漫大活劇バンピートロット』というゲームに昔、熱中した。このゲームの世界観も同様で、それが原因かもしれない。私も、誰が読んでいるのか全く分からないこのブログで文章を書きながら、自分を少しずつ発見していると、いつも感じている。