タグ別アーカイブ: .NET Framework

Cox回帰の結果を使ってMCMCしてみようとした

まず、Cox回帰は以下のように定義されます。すなわち、ある瞬間での瞬間死亡率を積分して、累積死亡率を出します。このとき、h_0(t)がベースラインハザード関数と呼ばれ、時間依存の死亡率を表します。そして、それに共変量の影響を加えますが影響は対数線形モデルで書けること、さらに時間によらず一定なこと(比例ハザード性)が仮定されています。

h(t)=\lim_{\Delta t \to 0} \frac{S(t)-S(t - \Delta t)}{\Delta t} \cdot \frac{1}{S(t)}
H(t)=\int_0^{t} h(t)du=-\log S(t)
h(t)=h_0(t) \cdot \mathrm{e}^{(\alpha + \beta_{1} X_{1} + \beta_{2} X_{2} + \ldots + \beta_{p} X_{p})} = h_0(t) \cdot RISK

たとえば、仮にCox回帰でゲームをやめてしまうモデルを作ったとすれば、お客さんがゲームをやめてしまう可能性はたとえば性別だとかそういったものを説明変数とした場合、その影響は時間によらず一定なことが求められます。というか、そうであることが仮定されているモデルだからですが。

裏を返せば、説明変数の影響は時間によらず一定なので、共変量の影響部分だけを切り離してスコア化できるモデルともとらえられると思います。つまり、先の式で言うと\mathrm{e}^{(\alpha + \beta_{1} X_{1} + \beta_{2} X_{2} + \ldots + \beta_{p} X_{p})}の部分を切り離して顧客のリスクとしてとらえられます。

もともと、Cox回帰は生存時間分析の一角で医学で用いられてきたため、用語として瞬間死亡確率、累積死亡確率というなかなかヘビーな単語が頻出します。今日では、離反予測などでも使われているのでイベントヒストリー分析なとども言います。

累積死亡率ですから、その区間の前までの累積死亡率からその前までの累積死亡率を引けばその区間での死亡率が求められるはずです。そうすると、その死亡率でさいころを振ってみてその区間で死亡するかどうかが求められます。その結果を集計すればその結果はほぼほぼ集計レベルでは元と近似するはずです。

ベースラインハザード関数の部分はノンパラメトリックなの推定ですから推定部分は定数としておくなり、データベースに格納するなりが考えられます。パラメトリックな共変量からの計算部分は関数としておけます。たとえば、Rのサンプルデータのkidneyのモデルはこのようにおけるでしょうか。

let KidneyModel sex diseaseGN diseaseAN diseasePKD = 
    let x00 = sex - 1.73684210526316  in
        let x01 = diseaseGN - 0.236842105263158 in
            let x02 = diseaseAN - 0.315789473684211  in
                let x03 = diseasePKD - 0.105263157894737  in
                    System.Math.Exp(-1.47739255480755 * x00 + 0.139220775300325  * x01 + 0.413157096502011  * x02 + (-1.36707122522468 ) * x03)

このコードはベースラインハザード関数の部分、先ほどの式ではh_0(t)の部分を含んでいないので式としては不完全です。最初で共変量から引いているのはRのcoxphは内部で共変量から平均を引き算しているため同様な処理を入れています。当然ながら、書き下した関数とRの結果がほぼほぼ誤差込みで一致するのは確認済みです。

Rのモデルの情報はこんな感じです。

Call:
coxph(formula = Surv(time, status) ~ sex + disease, data = kidney, 
    model = T)


             coef exp(coef) se(coef)     z       p
sex        -1.477     0.228    0.357 -4.14 3.5e-05
diseaseGN   0.139     1.149    0.364  0.38    0.70
diseaseAN   0.413     1.512    0.336  1.23    0.22
diseasePKD -1.367     0.255    0.589 -2.32    0.02

Likelihood ratio test=17.6  on 4 df, p=0.0015
n= 76, number of events= 58 

このとき、KidneyBaseline t なる時刻tを引数としてもち、時刻tのベースラインハザード曲線を返す関数があったならば、時刻tにおける累積死亡確率は

h(t)=h_0(t) \cdot \mathrm{e}^{(\alpha + \beta_{1} X_{1} + \beta_{2} X_{2} + \ldots + \beta_{p} X_{p})}

に基づき、(KidneyModel sex diseaseGN diseaseAN diseasePKD) * (KidneyBaseline t)になります。従って。累積死亡確率を求める関数 KidneyCumsumは以下のように定義されます。

let KidneyCumsum t sex diseaseGN diseaseAN diseasePKD = (KidneyModel sex diseaseGN diseaseAN diseasePKD) * (KidneyBaseline t)

その上で、区間t1からt2における死亡確率は(KidneyCumsum t2) – (KidneyCumsum t1)になります。従って、KidneyDeathProbabilityなる関数を定義するとすれば

let KidneyDeathProbability t1 t2 sex diseaseGN diseaseAN diseasePKD = (KidneyCumsum t2) - (KidneyCumsum t1)

となります。

これをサンプルデータ一つ一つについて、順次計算し、同じ区間について集計すれば平均的なその区間における死者を推定することができます。Cox回帰は比例ハザード性を仮定しているので共変量の効き具合は時間にかかわらずという強力な仮定があるのでこの、仮定が崩れていると使えませんので注意がいります。

これらを踏まえてMCMCを実装するとこのようになります。

open Kidney;

type Row = {sex: float; diseaseGN: float; diseaseAN: float; diseasePKD: float; mutable dead: bool}

[]
let main argv = 
    let sr = new System.IO.StreamReader("kidneyp.csv", System.Text.Encoding.ASCII)
    let mutable buffer = sr.ReadLine()
    let getGN (buffer:string) =
        match buffer.Replace("\"", "") with
        | "GN" -> 1.0
        | _ -> 0.0
    let getAN (buffer:string) =
        match buffer.Replace("\"", "") with
        | "AN" -> 1.0
        | _ -> 0.0
    let getPKD (buffer:string) =
        match buffer.Replace("\"", "") with
        | "PKD" -> 1.0
        | _ -> 0.0

    let parse (buffer:string) =
        let items = buffer.Split [|','|]
        {sex = System.Double.Parse items.[5];diseaseGN = getGN items.[6];diseaseAN = getAN items.[6];diseasePKD = getPKD items.[6];dead = false}

    let patients = [|
        while not sr.EndOfStream do
            buffer <- sr.ReadLine()
            yield parse buffer
    |]

    let rg = new System.Random()

    let results = [|
        for i in 1..(KidneyModel.KidneyBaselineSize - 1) do
            let ps = [|
                for patient in patients do
                    if patient.dead = false then
                        let p = KidneyModel.KidneyDeathProbability (i - 1) i patient.sex patient.diseaseGN patient.diseaseAN patient.diseasePKD
                        let r = rg.NextDouble()
                        let c = match (r <= p) with
                        | true -> 1.0
                        | _ -> 0.0
                        if c > 0.0 then
                            patient.dead <- true
                        yield c
                    else
                        yield 0.0
                done
            |]
            yield Array.sum(ps)
        done
    |]

    let sw = new System.IO.StreamWriter("dead.csv", false, System.Text.Encoding.ASCII);
    for result in results do
        sw.WriteLine(result)
    done
    sw.Close()

    printfn "%A" argv
    0 // 整数の終了コードを返します

シミュレーションの結果は以下の通りです。

集計結果

集計結果

データサイエンティスト養成読本 R活用編

いささか、旧聞にはなるが上記、書籍を書いておきます。基本的には前に刊行されたムックの続きでRなどにフォーカスしたムックです。R活用編となってはいますが、商業書籍では珍しいJuliaの入門記事やRを.NET環境から利用するR.NETの記事など割と面白いものになっています。特にR.NETを含めた.NETのライブラリ群の配布は現在はNuGetであることは.NETの開発者の諸兄には自明なのですが、意外と気づかない向きもありますので書いておきます。

さて、私の目的はJuliaの入門記事であったわけですが、文法の解説に関しては割りとしっかり書かれているように思われます。とはいえ、ムックの一セクションなので食い足りないのも確かです。書籍ベースだと後は洋書の”Seven More Languages in Seven Weeks: Languages That Are Shaping the Future”と先に書評を書いた、同人誌の実例で学ぶJulia言語入門 v0.3.3/v0.4.0-dev 対応版あたり。

Jubatus Client for .NETを実装中

オンライン学習フレームワークとして知られる、Jubatusの.NET用のクライアントを実装中です。MessagePack RPC for CLIで楽ができるかと思ったのですが、どうも更新されていないようでMessagePackのシリアライザとの間でバージョンが厳しいことになっているようでRPCから再実装中です。ただ、現在作っているコードはJubatusとの通信しか想定していないのでクライアントライブラリにべったりと実装されています。現在は、プライベートなレポジトリでせこせこといっています。

とりあえず、get_configの実装はできているのでエビデンス代わりにぽとぺたと。もう少し、作ったらパブリックなレポジトリに展開します。

{
  "method": "PA",
  "converter": {
    "string_filter_types": {
      "detag": { "method": "regexp", "pattern": "<[^>]*>", "replace": "" }
    },
    "string_filter_rules": [
      { "key": "message", "type": "detag", "suffix": "-detagged" }
    ],
    "num_filter_types": {},
    "num_filter_rules": [],
    "string_types": {},
    "string_rules": [
      { "key": "message-detagged", "type": "space", "sample_weight": "bin", "global_weight": "bin"}
    ],
    "num_types": {},
    "num_rules": []
  },
  "parameter": {}
}