以下のイベントで話しました。
発表資料は以下です。
前にやった解析において、最終的なモデルにいたるまでのプロセスとモデリングのコツを多めにした内容にしました。
自分以外の発表内容も大変面白くて楽しめました。運営と発表者に感謝です。ありがとうございました!
この記事の続きです。
ここではPRMLの10.1.3項の一変数ガウス分布の例題(WikipediaのVariational_Bayesian_methodsのA basic exampleと同じ)をSymPyで解きます。すなわちデータが
に従い*1、とが、
に従うという状況です。ここでデータ()が得られたとして事後分布を変分ベイズで求めます。
まずはじめに、上記の確率モデルから同時分布を書き下しておきます。
なので、
となります。
この問題は単純なので事後分布は厳密に求まるのですが、ここでは変分ベイズで解きます。すなわち、事後分布をで近似します。さらにと因子分解可能と仮定します。そして、前の記事の最後の2つの式を使って、とが収束するまで繰り返し交互に更新して求めるのでした。以下ではこれをSymPyでやります。
from sympy import * from sympy.stats import * init_printing(use_unicode=True) from IPython.core.interactiveshell import InteractiveShell InteractiveShell.ast_node_interactivity = 'all'
y, mu, mu0 = symbols('y mu mu0', real=True) Y_vec = symbols('Y1:4', real=True) tau, lambda0, a0, b0 = symbols('tau lambda0 a0 b0', positive=True)
symbols
関数で作成しておく必要があります。real=True
と指定することで、実数と仮定することができます。何も指定しなければ複素数になります。このように仮定を入れておかないと、のちの式変形や積分がうまくいかない場合があります。positive=True
と指定することで、正の実数だと仮定することができます。なお、SymPyでは要素数やデータ数をN
とするような一般の場合の式変形は基本的に難しいです。しかし具体的な値に決めれば実行できます。そこで、ここでは2行目でデータ数をY1
,Y2
,Y3
の3
個として先に進めます。あとで値を色々変えて試すとN
の場合の見当がつくので、そこから一般の場合を証明することもできます。
p_y = density(Normal('', mu, 1/sqrt(tau)))(y) p_mu = density(Normal('', mu0, 1/sqrt(lambda0*tau)))(mu) p_tau = density(Gamma('', a0, 1/b0))(tau)
sympy.stats
には確率分布の密度関数の式がありますので、それを使っています。ここではデータ1つあたりのy
の分布とmu
とtau
の事前分布を定義しています。
SymPyの正規分布はNormal(平均, 標準偏差)
なので、精度であるtau
を1/sqrt(tau)
を代入しています。また、PRMLやWikipediaのガンマ分布はGamma(shape, rate)
である一方*2、SymPyのガンマ分布はGamma(shape, scale)
なので、1/b0
を代入しています。
integrate(p_mu, (mu, -oo, oo))
simplify(integrate(p_tau, (tau, 0, oo)))
試しにmu
の分布をからまで積分してみましょう。期待通り1が返ります。tau
の分布でも同様に積分すると1にならずに整理されていない式が返ってきますが、simplify
関数で整理すると1になります。
前の記事の最後の2つの式でやっていることを日本語で書くと以下です。
そこでまず同時分布の対数を準備します。
log_p = sum([log(p_y.subs(y, x)) for x in Y_vec]) + log(p_mu) + log(p_tau) log_p = simplify(log_p) log_p
expr.subs(y, x)
は式expr
のy
にx
を代入します。次に積分にすすみます。
を含まない項にを含まない分布を掛けてで積分したところで、やはりに関係がない定数になります。定数は最後に規格化して求めればよいので、途中の計算はなるべく簡単になるように余計な項を取り除きます。これがSymPyで計算をうまくさせるポイントになります。
log_p_for_mu = integrate(diff(log_p, mu), mu) log_p_for_mu = collect(log_p_for_mu, mu) log_p_for_mu
log_p_for_tau = integrate(diff(log_p, tau), tau) log_p_for_tau = collect(log_p_for_tau, tau) log_p_for_tau
で積分してを求める方も同様なのでそうしておきます。
SymPyの練習のため、事前分布から積分を1回実行してとを求めるところをやってみます。
log_q1_mu = integrate(log_p_for_mu * p_tau, (tau, 0, oo))
log_q1_mu
log_q1_mu = simplify(log_q1_mu)
log_q1_mu
log_q1_mu = collect(expand(log_q1_mu), mu)
log_q1_mu
simplify
していますが、解析者が意図しない形になることはよくあります。ここでは、expand
してcollect
することでmu
の多項式にしています。log_q1_mu
の式はの二次関数のマイナスなので、このすぐあとのq1_mu
は正規分布になることが分かります。共役事前分布を使っているからです。規格化定数をもとめて規格化しましょう。
q1_mu = exp(log_q1_mu)
const = simplify(integrate(q1_mu, (mu, -oo, oo)))
const
q1_mu = 1/const * exp(log_q1_mu)
q1_mu
const
が規格化定数になります。以下の部分です。
q1_mu
は規格化された分布のです。以下になります。
同じようにq1_tau
を求めます。変分ベイズの手順としては、上で求めたばかりのを掛けてで積分します。しかしSymPyではその計算は重くて実行できないので、ここではの事前分布p_mu
を使ってq1_tau
を求めてみます。
log_q1_tau = integrate(log_p_for_tau * p_mu, (mu, -oo, oo)) log_q1_tau log_q1_tau = integrate(diff(log_q1_tau, tau), tau) log_q1_tau log_q1_tau = collect(log_q1_tau, tau) log_q1_tau
このすぐあとのq1_tau
はガンマ分布になることが分かります。これも共役事前分布を使っているからです。
q1_tau = logcombine(exp(log_q1_tau)) q1_tau # const = integrate(q1_tau, (tau, 0, oo)) # const # q1_tau = 1/const * q1_tau # q1_tau
logcombine
関数を使うことでをにします。simplify
関数だとこの変形をやってくれないことがあります。q1_tau
は以下です。
このの肩にのっているの係数が負だとSymPyが分からないから積分できないようです。ちなみにこのあたりはMathematicaの方が圧倒的に賢くて、例えば以下の入力できちんと積分できます。
Integrate[tau^(a+1)*Exp[tau * (-1/2*x^2 + x*mu - 1/2* y^2 + y*mu - b - mu^2)], {tau, 0, Infinity}, Assumptions -> {b > 0, a > 0, Element[x, Reals], Element[y, Reals], Element[mu, Reals]} ]
これをうまく積分させるには、の係数が負であることを確認してから変数で置き換えて実行します。
まずの係数が負であることを確認します。
coef = collect(log_q1_tau, tau).coeff(tau) coef sol = solve(diff(coef, Y_vec[0]), Y_vec[0])[0] sol #=> mu0 replacements = [(var, sol) for var in Y_vec] coef.subs(replacements) #=> -b0
coef
を取得しています。coef
の最大値が負であることを示せばOKです。まずはY_vec[0]
についてcoef
が最大になる値を探します。それは微分して0(&2階微分が負)になる点を求めればOKです。Y_vec[0]
と他のY_vec[*]
は区別がある形ではないので、Y_vec[*]
についても同じ点でcoef
が最大となります。-b0
と分かるので、の係数は負であることがわかります。次に変数で置き換えて積分します。
xi = symbols('xi', positive=True) const = simplify(integrate(tau**(a0+1)*exp(-xi*tau), (tau, 0, oo))) const = const.subs(xi, -coef) const q1_tau = 1/const * q1_tau q1_tau
const
が規格化定数になります。以下の部分です。
q1_tau
は正規化された分布のです。以下になります。
このように解析解を求めることはコンセプトの理解に役立ちます。一方で、積分を繰り返して事後分布が収束するか確認するようなことは数値的に求めた方が分かりやすいです。
仮に得られたデータY_vec
は1.1
,1.0
,1.3
とします。また、事前分布はa0 = 1
, b0 = 1
, mu0 = 0
, lambda0 = 1
とします。
replacements = [(a0, 1), (b0, 1), (mu0, 0), (lambda0, 1)] data_vec = [1.1, 1.0, 1.3] replacements.extend([(var, val) for var, val in zip(Y_vec, data_vec)]) log_p_for_mu_subs = log_p_for_mu.subs(replacements) log_p_for_tau_subs = log_p_for_tau.subs(replacements) [log_p_for_mu_subs, log_p_for_tau_subs]
の初期分布をp_tau
として、を求める→を求める→を求める→...と7回ほど繰り返してみます。
q_tau = N(p_tau.subs(replacements)) q_tau for i in range(7): log_q_mu = N(integrate(log_p_for_mu_subs * q_tau, (tau, 0, oo))) const = N(integrate(exp(log_q_mu), (mu, -oo, oo))) q_mu = 1/const * exp(log_q_mu) log_q_tau = N(integrate(log_p_for_tau_subs * q_mu, (mu, -oo, oo))) const = N(integrate(exp(log_q_tau), (tau, 0, oo))) q_tau = 1/const * exp(log_q_tau) [q_mu, q_tau]
N
関数は数値による近似を求める関数です。7回ほどの繰り返しのあとでほぼ収束していそうなことがわかります。
最後に求めた事後分布(の近似)を可視化してみましょう。SymPyにもsympy.plotting
やsympy.plotting.plot
が存在するのですが、ちょっと凝った図を書こうとするとすぐ厳しくなってしまいます。そこで、得られた事後分布をlambdify
関数で関数化し、NumPyとMatplotlibで描くのが拡張性が高くてオススメです。
from sympy.utilities.lambdify import lambdify import numpy as np import matplotlib.pyplot as plt delta = 0.05 x = np.arange(-1.0, 3.0, delta) y = np.arange(0.0, 6.0, delta) X, Y = np.meshgrid(x, y) func = lambdify((mu, tau), q_mu * q_tau, 'numpy') Z = func(X, Y) plt.figure() CS = plt.contour(X, Y, Z) plt.clabel(CS, inline=1, fontsize=10)
Enjoy!
北大電子研の佐藤勝彦氏に感謝します。僕が院生の頃に輪読していた ニコリス プリゴジーヌ『散逸構造』の例題をMathematicaで10分ぐらいで一般解を求めるという衝撃のデモを見せてもらい、その後もたまにMathematicaを教えてもらい、数式処理を学ぶきっかけをもらいました。
2017年8月9日に国立がん研究センターは、がん治療拠点の約半数にあたる全国188の病院について、癌患者の5年後の生存率データを初めて公表しました(毎日新聞の記事)。報告書は国立がん研究センターが運営するウェブサイトからダウンロードできます(ここ)。報告書をダウンロードしようとすると注意点を記したポップアップが表示されます。大切な部分を抜粋すると以下です。
本報告書には、施設別の生存率を表示していますが、進行がんの多い少ない、高齢者の多い少ないなど、施設毎に治療している患者さんの構成が異なります。そのため、単純に生存率を比較して、その施設の治療成績の良し悪しを論ずることはできません。
一般に高齢者が多い病院ほど、進行癌(ステージが進んだ癌)が多い病院ほど、その病院の生存率は下がるわけです。それならば、統計モデリングで年齢と進行度(ステージ)の影響を取り除いて(専門的な言葉で言えば「調整して」)病院の良し悪しを論じてみようというのがこの記事の内容になります。しかし、あくまでも「手法」の節で書いた仮定のもとでの推定結果であり、真実として断定するものではありませんのでその点はご理解ください。
なお、影響を取り除く前の(カプランマイヤー法で算出した)実測生存率のヒストグラムは以下になります。病院によってかなりばらつきがありそうに見えます。何もしないで比べると「国立研究開発法人国立がん研究センター中央病院」や「がん研有明病院」など大きな病院の生存率が高いです。患者の平均年齢と平均的な進行度が低いためと思われます。
それでは結果から先に述べます。
以下では病院ごとの生存率や手術率を比較するために、癌種t
・病院h
における男性比率を0.5
・平均年齢を60
・平均進行度を2.5
(おおよそステージIIに相当)に仮に固定して議論をすすめます。
=== 2018.2.6 追記 ===
この部分が分かりにくかったようなので補足します。今知りたいことは「あなたが60歳でステージIIだったらどこへ行くのが5年生存率が高いか」です。しかし、病院によっては若い人が多かったり、重症な人が多かったりして、そのままの生存率を単純に比較できません。そこで、統計モデリングによってすべての病院の平均年齢と平均進行度を仮想的に揃えることで、病院にだけ依存する生存率の差が残るわけです。ここでは一例として平均年齢を60
・平均進行度を2.5
に固定して算出しています。とにかく揃えればよいので、仮に平均年齢を50
・平均進行度を1.5
に固定してもよいです(この場合、生存率は全体的に高い方へ少しずれますがランキングは変わりません)。
生存率の推定値(MCMCサンプルの中央値)のヒストグラムは以下になります。
驚くべきことに胃・大腸・肺においてはほとんど病院の差がありません。がん診療連携拠点病院、さすがです。肝癌と乳癌は少し病院によって差がありますので、それぞれbest10を紹介します。ただし、病院による差は提供されたデータ以外の要因(他の病気をもつ患者が多い少ないなど)も含まれます。検討によってそういう情報が重要と分かれば、今後データとして収集して公表する必要があると思います。
黒点は(MCMCサンプルの)中央値、横に伸びる線は80%ベイズ信頼区間です。中央値の高い順に並び替えています。ベイズ信頼区間を見ると分かる通り、病院による差はそこまで明らかではありません。
「大垣市民病院」があまりに良いので病院のコメントを読むと「4.肝(C22)は切除症例に限る。」と書いてありました。これだと切除していない患者さんとその死亡数が考慮されていないので、フェアではなく正しいランクとは言えません。ただし、他の癌種における結果を見ると、全てを考慮してもかなり良い生存率である可能性はあります。
仮に「大垣市民病院」を除くと「信州大学医学部附属病院」が1位となります。
グラフの見方は肝癌の場合と同じです。
1位は「愛知県がんセンター中央病院」です。圧倒的な手術率もさることながら、手術以外の影響(病院による効果)も1位です。謎の民間療法に頼らず、こういうところに入院したいものです。
=== 2018.12.8 追記 ===
識者の方から、がんセンターは一般的に合併症のないピュアな癌患者を受け入れる傾向が強いので生存率が高いのでは、とコメントをいただきました。合併症のデータも重要で調整する必要があるかもしれません。今後の課題と思います。
点は中央値、横に伸びる線は80%ベイズ信頼区間です。例えば、胃癌の年齢_10
のOdds ratioの中央値は約0.5ですが、これは「ある病院における患者の平均年齢が10上がると『生存率/死亡率』が約0.5倍になる」ことを意味します。性別_0.1
は男性比率が0.1上がることを意味します。他も同様です。
結果を見ると、予想通り性別の影響はあまりなくて年齢や進行度の影響が強いです。特に胃癌においては年齢の影響が強く、乳癌においては進行度の影響が強いことが示唆されています。
グラフの見方は上の場合と同様です。肝癌において進行度の影響が低そうなのが意外に思いました。このあたり、病院によって生存率がばらついていることの原因があるのかもしれません。
ここでは記しませんが、他にも色々知ることができました。例を挙げます。
このような報告書を公表するのは英断であり、 “がん診療連携拠点病院ががん患者さんの治療に透明性を確保し、拠点病院全体として責任をもって取り組んでいる意気込み” (報告書からの抜粋)をまさに感じることができました(なかには難癖つけて公表を拒否している病院があり最後にまとめました)。一点、ちょっと残念なのは報告書がpdf(580ページほど)であり、利活用が大変しにくいことでした。今回は素敵なRによるスクレイピング入門を出している株式会社ホクソエムの市川さんと牧山代表取締役に担当してもらい、Rのtabulizer
パッケージを使って抽出しました。
この報告書には188の病院が含まれていて、各病院について5つの癌(胃、大腸、肝、肺、女性乳房)のデータが含まれています。1つの病院の1つの癌のデータ例を以下に示します(数値は架空のものです)。
個人情報保護の観点から、「1人以上10人以下」の場合に-
(ハイフン)に置き換えられています。これがこのデータの解析の難易度を大幅に引き上げています。
ハイフンへの対策は以下のようにしました。
打ち切りを考慮してカプランマイヤー法で算出した生存率を使う方法がまず考えられますが、対象数の大小に由来する生存率の推定幅をうまく組み込むのが容易ではありません。また、がん診療連携拠点病院の意気込みが凄くて、現時点で生存状況把握割合は95%以上が多く、今後は100%に近づきそうです。これらを踏まえて、対象数と5年後生存数を用いた二項ロジスティック回帰にしました。*1
元の報告書では、対象数が50人未満の場合、定された生存率の信頼性が低くなるため公表しないとしてハイフンになっています。本解析においても対象数が少ない場合、年齢や進行度のデータがほぼすべてハイフンとなって、列挙する場合の数が非常に多くなります。そこに推定の時間を取られるのは本意ではないため、同様に50人未満の場合は解析対象から除外しました。
使用したモデルは以下になります。
青色の項は女性乳房以外の癌において存在し、赤色の項は肺癌のみに存在します。inv_logit
はロジスティック関数です。ここでデータは以下です。
ちなみに、発見経緯のデータは進行度に反映されると仮定して今回は使用していません。
推定するパラメータは以下です。推定にはStanとRを用いました。
年齢と癌の進行度の影響を取り除くため、これらを説明変数とした二項ロジスティック回帰にするのがポイントです。さらに手術率も年齢や進行度の影響をうけると考えられるので上記のモデルとなっています。とがそれぞれ別の多変量正規分布から生成されるモデルも試してみましたが、現段階では少しデータが足りないようで、MCMCが収束しませんでした。また探索的な解析からは、年齢や進行度によって生存率や手術率が指数関数的に落ちるのではなく線形に近いことが示唆されたため、今回は非線形な関数を含まないモデルとしました。実際のコードは上記の数式に加えてハイフンの処理が入るのでかなり複雑になります。今回はコードの解説を省略します。
最後にモデルの事後予測チェックの図を以下に載せておきます。生存数の実測値と予測値(点は中央値、縦に伸びる線は80%ベイズ信頼区間)は以下です。
生存率のカプランマイヤー法による実測値と予測値(点は中央値、縦に伸びる線は80%ベイズ信頼区間)は以下です。肝癌は少しあてはまりが悪いですが、全体的に問題なく推定できていると思いました。
手術件数の実測値と予測値(点は中央値、縦に伸びる線は80%ベイズ信頼区間)は以下です。
その他にも癌, 病院における「平均年齢」と「平均進行度」の推定などがうまくいっているかもチェック済みです。
がん診療連携拠点病院と国立がん研究センターがん対策情報センターがん登録センター院内がん登録室に感謝します。相談に乗ってくれた@Med_KUさんと@happyningenさんに感謝します。
全ての癌種の最新かつ時系列のデータで再解析がしたいです。もし患者レベルのデータが利用できれば、さらに良い解析ができます。いつか解析できることを楽しみにしています。
治療のレベルはお察し。毎回更新したいです。
*1:打ち切りがある場合は5年後生存数を正確に得ることはできませんが、生存状況把握割合が非常に高いので、打ち切りに含まれる生存数を「打ち切り数*既知の生存数/対象数」で近似的に算出しました。
*2:本来は患者の年齢と生存か否かが結び付けれればさらに良い解析ができますが、現段階の情報ではこうやって多少強引に平均年齢を算出する他なさそうです。年代の人数は割と山型の分布なので、これで良さそうです。
*3:年代と同様です。しかし、癌種によっては「ステージIとIVは多いけどIIとIIIは少ない」ということがあり、進行度を要約するために平均+正規分布を使うのではなく、他の指標と分布を使った方がよいかもしれません。検討の余地があります。