StatModeling Memorandum

StatModeling Memorandum

StanとRとPythonでベイズ統計モデリングします. たまに書評.

統計モデリングで癌の5年生存率データから良い病院を探す

現実のデータ Stan

概要

2017年8月9日に国立がん研究センターは、がん治療拠点の約半数にあたる全国188の病院について、癌患者の5年後の生存率データを初めて公表しました（毎日新聞の記事）。報告書は国立がん研究センターが運営するウェブサイトからダウンロードできます（ここ）。報告書をダウンロードしようとすると注意点を記したポップアップが表示されます。大切な部分を抜粋すると以下です。

本報告書には、施設別の生存率を表示していますが、進行がんの多い少ない、高齢者の多い少ないなど、施設毎に治療している患者さんの構成が異なります。そのため、単純に生存率を比較して、その施設の治療成績の良し悪しを論ずることはできません。

一般に高齢者が多い病院ほど、進行癌（ステージが進んだ癌）が多い病院ほど、その病院の生存率は下がるわけです。それならば、統計モデリングで年齢と進行度（ステージ）の影響を取り除いて（専門的な言葉で言えば「調整して」）病院の良し悪しを論じてみようというのがこの記事の内容になります。しかし、あくまでも「手法」の節で書いた仮定のもとでの推定結果であり、真実として断定するものではありませんのでその点はご理解ください。

なお、影響を取り除く前の（カプランマイヤー法で算出した）実測生存率のヒストグラムは以下になります。病院によってかなりばらつきがありそうに見えます。何もしないで比べると「国立研究開発法人国立がん研究センター中央病院」や「がん研有明病院」など大きな病院の生存率が高いです。患者の平均年齢と平均的な進行度が低いためと思われます。

f:id:StatModeling:20201106161444p:plain

結果

それでは結果から先に述べます。

以下では病院ごとの生存率や手術率を比較するために、癌種t・病院hにおける男性比率を0.5・平均年齢を60・平均進行度を2.5（おおよそステージIIに相当）に仮に固定して議論をすすめます。

=== 2018.2.6 追記 ===

この部分が分かりにくかったようなので補足します。今知りたいことは「あなたが60歳でステージIIだったらどこへ行くのが5年生存率が高いか」です。しかし、病院によっては若い人が多かったり、重症な人が多かったりして、そのままの生存率を単純に比較できません。そこで、統計モデリングによってすべての病院の平均年齢と平均進行度を仮想的に揃えることで、病院にだけ依存する生存率の差が残るわけです。ここでは一例として平均年齢を60・平均進行度を2.5に固定して算出しています。とにかく揃えればよいので、仮に平均年齢を50・平均進行度を1.5に固定してもよいです（この場合、生存率は全体的に高い方へ少しずれますがランキングは変わりません）。

癌種ごとの生存率のヒストグラム

生存率の推定値（MCMCサンプルの中央値）のヒストグラムは以下になります。

f:id:StatModeling:20201106161448p:plain

驚くべきことに胃・大腸・肺においてはほとんど病院の差がありません。がん診療連携拠点病院、さすがです。肝癌と乳癌は少し病院によって差がありますので、それぞれbest10を紹介します。ただし、病院による差は提供されたデータ以外の要因（他の病気をもつ患者が多い少ないなど）も含まれます。検討によってそういう情報が重要と分かれば、今後データとして収集して公表する必要があると思います。

　肝癌における生存率best10の病院

f:id:StatModeling:20201106161456p:plain

黒点は（MCMCサンプルの）中央値、横に伸びる線は80%ベイズ信頼区間です。中央値の高い順に並び替えています。ベイズ信頼区間を見ると分かる通り、病院による差はそこまで明らかではありません。

「大垣市民病院」があまりに良いので病院のコメントを読むと「4.肝（C22）は切除症例に限る。」と書いてありました。これだと切除していない患者さんとその死亡数が考慮されていないので、フェアではなく正しいランクとは言えません。ただし、他の癌種における結果を見ると、全てを考慮してもかなり良い生存率である可能性はあります。

仮に「大垣市民病院」を除くと「信州大学医学部附属病院」が1位となります。

　乳癌における生存率best10の病院

f:id:StatModeling:20201106161452p:plain

グラフの見方は肝癌の場合と同じです。

1位は「愛知県がんセンター中央病院」です。圧倒的な手術率もさることながら、手術以外の影響（病院による効果）も1位です。謎の民間療法に頼らず、こういうところに入院したいものです。

=== 2018.12.8 追記 ===

識者の方から、がんセンターは一般的に合併症のないピュアな癌患者を受け入れる傾向が強いので生存率が高いのでは、とコメントをいただきました。合併症のデータも重要で調整する必要があるかもしれません。今後の課題と思います。

生存率への年齢・進行度・手術率の影響（オッズ比）

f:id:StatModeling:20201106161441p:plain

点は中央値、横に伸びる線は80%ベイズ信頼区間です。例えば、胃癌の年齢_10のOdds ratioの中央値は約0.5ですが、これは「ある病院における患者の平均年齢が10上がると『生存率/死亡率』が約0.5倍になる」ことを意味します。性別_0.1は男性比率が0.1上がることを意味します。他も同様です。

結果を見ると、予想通り性別の影響はあまりなくて年齢や進行度の影響が強いです。特に胃癌においては年齢の影響が強く、乳癌においては進行度の影響が強いことが示唆されています。

手術率への年齢・進行度の影響（オッズ比）

f:id:StatModeling:20201106161437p:plain

グラフの見方は上の場合と同様です。肝癌において進行度の影響が低そうなのが意外に思いました。このあたり、病院によって生存率がばらついていることの原因があるのかもしれません。

その他

ここでは記しませんが、他にも色々知ることができました。例を挙げます。

  • 手術の好きな病院・嫌いな病院
  • 手術率が比較的低いわりに生存率がまあまあ高い病院
  • 生存率worst10の病院

手法

データの抽出

このような報告書を公表するのは英断であり、 “がん診療連携拠点病院ががん患者さんの治療に透明性を確保し、拠点病院全体として責任をもって取り組んでいる意気込み” （報告書からの抜粋）をまさに感じることができました（なかには難癖つけて公表を拒否している病院があり最後にまとめました）。一点、ちょっと残念なのは報告書がpdf（580ページほど）であり、利活用が大変しにくいことでした。今回は素敵なRによるスクレイピング入門を出している株式会社ホクソエムの市川さんと牧山代表取締役に担当してもらい、Rのtabulizerパッケージを使って抽出しました。

データの難所

この報告書には188の病院が含まれていて、各病院について5つの癌（胃、大腸、肝、肺、女性乳房）のデータが含まれています。1つの病院の1つの癌のデータ例を以下に示します（数値は架空のものです）。

f:id:StatModeling:20201106161605p:plain

個人情報保護の観点から、「1人以上10人以下」の場合に-（ハイフン）に置き換えられています。これがこのデータの解析の難易度を大幅に引き上げています。

ハイフンへの対策は以下のようにしました。

  • 合計から数値が計算できる場合は計算して置き換える。
    • 例: 「観血的治療の実施」（＝外科的手術の実施）の有と無の合計が「対象数」（＝患者数）となるため、上記の例では「観血的治療の実施_無」は100-91=9と算出できる。
    • 例: 「(100.0 - 生存状況把握割合(%))*対象数」を四捨五入して「打ち切り数」が算出できる。
  • 場合の数を列挙し、統計モデリングに組み込む。

生存率はどの値を使うか

打ち切りを考慮してカプランマイヤー法で算出した生存率を使う方法がまず考えられますが、対象数の大小に由来する生存率の推定幅をうまく組み込むのが容易ではありません。また、がん診療連携拠点病院の意気込みが凄くて、現時点で生存状況把握割合は95%以上が多く、今後は100%に近づきそうです。これらを踏まえて、対象数と5年後生存数を用いた二項ロジスティック回帰にしました。*1

対象数が少ない場合への対応

元の報告書では、対象数が50人未満の場合、定された生存率の信頼性が低くなるため公表しないとしてハイフンになっています。本解析においても対象数が少ない場合、年齢や進行度のデータがほぼすべてハイフンとなって、列挙する場合の数が非常に多くなります。そこに推定の時間を取られるのは本意ではないため、同様に50人未満の場合は解析対象から除外しました。

統計モデル

使用したモデルは以下になります。

f:id:StatModeling:20201106161511p:plain

f:id:StatModeling:20201106161515p:plain

f:id:StatModeling:20201106161518p:plain

f:id:StatModeling:20201106161523p:plain

f:id:StatModeling:20201106161526p:plain

f:id:StatModeling:20201106161530p:plain

f:id:StatModeling:20201106161534p:plain

f:id:StatModeling:20201106161537p:plain

f:id:StatModeling:20201106161541p:plain

f:id:StatModeling:20201106161545p:plain

f:id:StatModeling:20201106161549p:plain

f:id:StatModeling:20201106161554p:plain

f:id:StatModeling:20201106161557p:plain

f:id:StatModeling:20201106161602p:plain

青色の項は女性乳房以外の癌において存在し、赤色の項は肺癌のみに存在します。inv_logitはロジスティック関数です。ここでデータは以下です。

  • t: 癌種のインデックス（tissue）。 Tが癌種の数で、ここでは T=5。
  • h: 病院のインデックス（hospital）。 Hが病院の数で、ここでは報告書にならって患者数が50以上のみを解析対象としたため、 H=183。
  • Surv: 生存数
  • N: 対象数（患者数）
  • Ope: 観血的治療の実施有の数（外科的手術をした人数）
  • Age: 各年代に属する人数
  • Stage: 各進行度に属する人数
  • Male: 性別が男性の人数
  • SC: 肺癌において、部位が小細胞である人数
  • Cutoff_{age}: 平均 \mu_{age}の正規分布を0.5, 0.6, 0.7, 0.8で切って5つの領域の面積をそれぞれ p_{age}に割り当てます。その閾値です。
  • Cutoff_{stage}: 平均 \mu_{stage}の正規分布を0.2, 0.3, 0.4で切って4つの領域の面積をそれぞれ p_{stage}に割り当てます。その閾値です。

ちなみに、発見経緯のデータは進行度に反映されると仮定して今回は使用していません。

推定するパラメータは以下です。推定にはStanとRを用いました。

  • p_{surv}: 生存率
  • p_{ope}: 手術率
  • q_{male}: 男性比率
  • q_{SC}: 肺癌における小細胞率
  • \mu_{age}: 癌 t, 病院 hにおける「平均年齢」。*2
  • \mu_{stage}: 癌 t, 病院 hにおける「平均進行度」。*3
  • p_{age}: 各年代の比率
  • p_{stage}: 各進行度の比率
  • r_{surv}: 生存率に対する説明変数以外の影響。病院差。
  • r_{ope}: 手術率に対する説明変数以外の影響。病院差。
  • b, a: 回帰係数
  • \sigma_{r_{surv}}: 生存率の病院差の標準偏差。データから推定します（階層モデル）。
  • \sigma_{r_{ope}}: 手術率の病院差の標準偏差。データから推定します（階層モデル）。

年齢と癌の進行度の影響を取り除くため、これらを説明変数とした二項ロジスティック回帰にするのがポイントです。さらに手術率も年齢や進行度の影響をうけると考えられるので上記のモデルとなっています。 r_{surv}と r_{ope}がそれぞれ別の多変量正規分布から生成されるモデルも試してみましたが、現段階では少しデータが足りないようで、MCMCが収束しませんでした。また探索的な解析からは、年齢や進行度によって生存率や手術率が指数関数的に落ちるのではなく線形に近いことが示唆されたため、今回は非線形な関数を含まないモデルとしました。実際のコードは上記の数式に加えてハイフンの処理が入るのでかなり複雑になります。今回はコードの解説を省略します。

最後にモデルの事後予測チェックの図を以下に載せておきます。生存数の実測値と予測値（点は中央値、縦に伸びる線は80%ベイズ信頼区間）は以下です。

f:id:StatModeling:20201106161503p:plain

生存率のカプランマイヤー法による実測値と予測値（点は中央値、縦に伸びる線は80%ベイズ信頼区間）は以下です。肝癌は少しあてはまりが悪いですが、全体的に問題なく推定できていると思いました。

f:id:StatModeling:20201106161508p:plain

手術件数の実測値と予測値（点は中央値、縦に伸びる線は80%ベイズ信頼区間）は以下です。

f:id:StatModeling:20201106161459p:plain

その他にも癌 t, 病院 hにおける「平均年齢」と「平均進行度」の推定などがうまくいっているかもチェック済みです。

後記

謝辞

がん診療連携拠点病院と国立がん研究センターがん対策情報センターがん登録センター院内がん登録室に感謝します。相談に乗ってくれた@Med_KUさんと@happyningenさんに感謝します。

展望

全ての癌種の最新かつ時系列のデータで再解析がしたいです。もし患者レベルのデータが利用できれば、さらに良い解析ができます。いつか解析できることを楽しみにしています。

難癖つけて公表しなかった残念な病院一覧

治療のレベルはお察し。毎回更新したいです。

*1:打ち切りがある場合は5年後生存数を正確に得ることはできませんが、生存状況把握割合が非常に高いので、打ち切りに含まれる生存数を「打ち切り数*既知の生存数/対象数」で近似的に算出しました。

*2:本来は患者の年齢と生存か否かが結び付けれればさらに良い解析ができますが、現段階の情報ではこうやって多少強引に平均年齢を算出する他なさそうです。年代の人数は割と山型の分布なので、これで良さそうです。

*3:年代と同様です。しかし、癌種によっては「ステージIとIVは多いけどIIとIIIは少ない」ということがあり、進行度を要約するために平均＋正規分布を使うのではなく、他の指標と分布を使った方がよいかもしれません。検討の余地があります。

逆温度1の事後分布のサンプルからWBICを計算する

Stan R 理論

この記事は以下のツイートを拝見してやってみようと思いました。

ツイートで言及されている渡辺先生の論文は以下です。

  • S Watanabe (2013) "A widely applicable Bayesian information criterion" Journal of Machine Learning Research 14 (Mar), 867-897 (pdf file)

この記事では、以前WAICとLOOCVの比較をした時に使った3つのモデル（重回帰、ロジスティック回帰、非線形回帰）において、「定義通りに算出したオリジナルのWBIC」と「近似式（上記論文の(20)式）で求めたWBIC」を比較してみました。

手法

case 1 重回帰

真のモデルは以下です。

　 Y \sim Normal(1.3 - 3.1 X_1 + 0.7 X_2, 2.5)

あてはめたモデルは以下です。

　 Y \sim Normal(b_1 + b_2 X_1 + b_3 X_2, \sigma)

データ点の数Nについては20,100を試しました。例としてN = 20の場合を説明します。まず乱数でデータX（すなわち X_1, X_2）を生成します。次にそのXの値を使ってYを生成しますが、以下の二つの場合について計算しました。

  • 1) MCMCサンプルの出方の違いによる影響: Yを1つだけ生成して固定し、MCMCの乱数の種を変えて200回推定を行い、 WBIC_{original}・ WBIC_{approx}のそれぞれの分布を確認した
  • 2) データの出方の違いによる影響: Yを1000通り生成し、 WBIC_{approx} - WBIC_{original}および相対的な差である (WBIC_{approx} - WBIC_{original})/WBIC_{original}の分布を確認した

事後分布の推定はStanで行いました。iter=11000, warmup=1000, chains=4で実行して合計40000個のMCMCサンプルを得ています。

case 2 ロジスティック回帰

手順は重回帰の場合と同じです。使用したモデルだけが異なります。真のモデルは以下です。

　 Y \sim Bernoulli(inv\_logit(0.8 - 1.1 X_1 + 0.1 X_2))

あてはめたモデルは以下です。

　 Y \sim Bernoulli(inv\_logit(b_1 + b_2 X_1 + b_3 X_2))

　 b_1,b_2,b_3 \sim Student\_t(4,0,1)

case 3 非線形回帰 ミカエリス・メンテン型

手順は重回帰の場合と同じです。使用したモデルだけが異なります。真のモデルは以下です。

　 Y \sim Normal(10.0 X / (2.0 + X), 0.8)

あてはめたモデルは以下です。

　 Y \sim Normal(m X / (k + X), \sigma)

　 k \sim Uniform(0, 12)

　 m \sim Uniform(0, 20)

case 3b 真のモデルが含まれない場合

あてはめたモデルが以下の場合も試しました。

　 Y \sim Normal(a + b X, \sigma)

結果

計算速度

Stanを使う場合、近似式のモデルの方がサンプリングが速いので、計算速度は近似式の方が少し速いです。どれだけ速くなるかはモデル依存で場合によりけりです。

MCMCサンプルの出方の違いによる影響

f:id:StatModeling:20201106161809p:plain

近似式の方はMCMCサンプルの出方によってかなりばらつくようです。また、少し値が低くなっています。

データの出方の違いによる影響

差

f:id:StatModeling:20201106161728p:plain

横軸は WBIC_{approx} - WBIC_{original}です。少しマイナスに偏った分布になりました。これが近似で捨てた項の影響なのか、Stanによるサンプリングの影響なのかは分かりません。

相対的な差

f:id:StatModeling:20201106161732p:plain

横軸は相対的な差である (WBIC_{approx} - WBIC_{original})/WBIC_{original} * 100です。ロジスティック回帰のN = 20の場合は、しばしば WBIC_{original}が0に近くなるので尾を引いています。Nが増えるに従って相対的な差は小さくなり、N = 100では±5%ぐらいに収まりそうです。

まとめ

Nが大きいときは近似式でスピードを重視しても大丈夫そう。でもNが小さいときは定義通り計算した方が無難に思えます。

ソースコード

case 1のソースコードを以下に載せます。

オリジナルのWBICを算出するためのStanコード

model/model1-ori.stanというファイル名とします。

data {
  int D;
  int N;
  matrix[N,D] X;
  vector[N] Y;
}

parameters {
  vector[D] b;
  real<lower=0> sigma;
}

transformed parameters {
  vector[N] mu;
  mu = X*b;
}

model {
  target += 1/log(N) * normal_lpdf(Y | mu, sigma);
}

generated quantities {
  vector[N] log_lik;
  for (n in 1:N)
    log_lik[n] = normal_lpdf(Y[n] | mu[n], sigma);
}

近似式でWBICを算出するためのStanコード

model/model1-apx.stanというファイル名とします。

data {
  int D;
  int N;
  matrix[N,D] X;
  vector[N] Y;
}

parameters {
  vector[D] b;
  real<lower=0> sigma;
}

transformed parameters {
  vector[N] mu;
  mu = X*b;
}

model {
  Y ~ normal(mu, sigma);
}

generated quantities {
  vector[N] log_lik;
  for (n in 1:N)
    log_lik[n] = normal_lpdf(Y[n] | mu[n], sigma);
}

各WBICを算出するRコード

library(rstan)

wbic_original <- function(log_lik) {
  wbic <- - mean(rowSums(log_lik))
  return(wbic)
}

wbic_approx <- function(log_lik) {
  b2 <- 1.0/log(ncol(log_lik))
  b1 <- 1.0
  log_denominator <- statnet.common::log_sum_exp(-(b2-b1)*(rowSums(-log_lik)))
  log_numerator   <- statnet.common::log_sum_exp(-(b2-b1)*(rowSums(-log_lik)) + log(rowSums(-log_lik)))
  wbic <- exp(log_numerator - log_denominator)
  return(wbic)
}

set.seed(123)
D <- 3
b <- c(1.3, -3.1, 0.7)
SD <- 2.5
N <- 100

X <- cbind(1, matrix(runif(N*(D-1), -3, 3), N, (D-1)))
Mu <- X %*% b
Y <- rnorm(N, Mu, SD)
data <- list(N=N, D=D, X=X, Y=Y)

sm_ori <- stan_model(file='model/model1-ori.stan')
sm_apx <- stan_model(file='model/model1-apx.stan')
fit_ori <- sampling(sm_ori, pars='log_lik', data=data, iter=11000, warmup=1000, seed=123)
fit_apx <- sampling(sm_apx, pars='log_lik', data=data, iter=11000, warmup=1000, seed=123)
wbic_ori <- wbic_original(rstan::extract(fit_ori)$log_lik)
wbic_apx <- wbic_approx(rstan::extract(fit_apx)$log_lik)
c(wbic_ori=wbic_ori, wbic_apx=wbic_apx)
  • 3～6行目: WBIC_{original}を算出します。この記事参照。
  • 8～15行目: WBIC_{approx}を算出します。途中で{statnet.common}パッケージのlog_sum_exp関数を使っています。前の記事のようにStanに含まれるlog_sum_exp関数を使っても構いません（全く同じ数値になります）。
  • 9行目: log_likはN_mcmc×N（データの数）のmatrix型ですのでncol(log_lik)でNを取得しています。

渡辺先生の論文の(20)式の通りに計算しようとすると、expの内側が50ぐらい以上の数値になるため、計算が不安定になります。そのため(20)式の両辺の対数をとって計算して、最後にexpをかませて戻します。MCMCを使っている場合、(20)式の左辺の対数は以下のように式変形できます。

f:id:StatModeling:20201106161736p:plainf:id:StatModeling:20201106161739p:plain

なので、まず易しい分母の方から計算すると、

f:id:StatModeling:20201106161743p:plain

f:id:StatModeling:20201106161746p:plain

f:id:StatModeling:20201106161750p:plain

f:id:StatModeling:20201106161753p:plain

f:id:StatModeling:20201106161758p:plain

分子の方も同様に、

f:id:StatModeling:20201106161801p:plain

f:id:StatModeling:20201106161805p:plain

これをそのまま実装しています。

Enjoy!

Stanの関数を使ってRを拡張して高速化する

R Stan

（C++に自動で変換される）Stanの関数を使ってRを拡張できる機能が、Stan/RStanの2.16で実装開始されて2.17でほぼ完成しました。Rを高速化するためにC++（とRcpp）はあまり書きたくないけれど、Stanの関数なら書いてもいいよという僕得な機能です。この記事ではその方法を簡単に紹介します。

元にした資料はRStanの開発者であるBenさんがStanCon2018で発表したこちらの資料です。

ここでは例として、以下の2つの関数をRで使えるようにしましょう。

  • 1) 機械学習分野でおなじみのlog_sum_exp関数
    • 引数はN個の正の実数
  • 2) データにemax modelという曲線をあてはめた場合の対数尤度を返す関数
    • 引数はデータ（N個のXとYのペア）とパラメータの値

手順は簡単で以下だけです。

  • functionsブロックだけ書いたstanファイルを用意する
  • R側でrstan::expose_stan_functions()を呼び出す

具体的には、以下の内容を例えばmy_functions.stanという名前で保存します。

functions {
  real stan_log_sum_exp(vector x) {
    return(log_sum_exp(x));
  }

  real emax_model_ll(vector Y, vector X, real e0, real emax, real ed50, real sigma) {
    int N = num_elements(Y);
    vector[N] mu = e0 + emax*X ./ (ed50 + X);
    real ll = normal_lpdf(Y | mu, sigma);
    return(ll);
  }
}
  • 3行目：log_sum_exp関数の方はすでにStanに用意されているものを使います。
  • 6行目：接尾の_llはlog likelihoodの意味です。
  • 7行目：たしかStan 2.14ぐらいから変数の宣言と定義を同時にできるようになりました。ここではデータポイントの数NをYから動的に取得しています。

次にR側でrstan::expose_stan_functions()を呼び出します。

rstan::expose_stan_functions('my_functions.stan')

するとfunctionsブロックに入っていた関数は全てR側で使えるようになります。

stan_log_sum_exp(c(0.1, 0.05, 0.8))
emax_model_ll(Y=c(1,10,12,24), X=c(0,5,10,10), e0=0.8, emax=30, ed50=5.1, sigma=3)

StanのコードはC++のコードに変換されるのでfor文はRと比べて非常に高速です。また、Stanのマニュアル（特にBuilt-In Functionsの章以降）を見ていただくとよいのですが、Stanはマニアックな分布の対数尤度を求める関数や、行列演算の関数（Eigen由来）や、log1m_expなどの指数対数を安定に計算する関数が豊富です。

これらのStanの長所を非常に簡単にRのコーディングに取り込めるのは最高です。まあ、少しばかりStanの関数やデータ型に慣れている必要がありますが。

Enjoy!