SurvivalAndDurationAnalysis_ja.ipynb

survivalanddurationanalysis.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "SurvivalAndDurationAnalysis_ja.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "toc_visible": true,
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/gist/AseiSugiyama/e0b875f4972481337277a661c0e10a0f/survivalanddurationanalysis.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "0EnOqUftToRD",
        "colab_type": "text"
      },
      "source": [
        "このノートブックは [Methods for Survival and Duration Analysis — statsmodels](https://www.statsmodels.org/dev/duration.html) の内容を日本語に翻訳し、社内向けの勉強会用に追記を行ったものです。"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ALs37FZEn2NA",
        "colab_type": "text"
      },
      "source": [
        "# 生存時間分析\n",
        "\n",
        "[statsmodels.duration](https://www.statsmodels.org/dev/duration.html#module-statsmodels.duration) は打ち切りのあるデータを扱うためのいくつかの標準的な手法を実装しています。これらの手法は主に、計測の開始から何らかの興味を持つイベントの発生までにかかる経過時間を含んでいます。典型的な例は医療研究における、被検者が何らかの症状の診断を受けてから、死亡 (あるいは、悪化や快癒) するまでの経過時間の分析です。\n",
        "\n",
        "現状では、右打ち切りデータのみ扱えます。右打ち切りデータとは、あるイベントがある時刻 t 以降に発生したことは知っているものの、具体的にどの時点でしたのかはわからないデータを言います。\n",
        "\n",
        "## 生存曲線の推定と推論\n",
        "\n",
        "`statsmodels.api.SurvfuncRight` クラスは、右打ち切りデータから生存関数の推定を行います。 `SurvfuncRight` はいくつかの推論処理を実装しており、生存分布の四分位点や、各点ごとの、または値域を通じた区間推定を行います。`duration.survdiff` 関数は検定を行い生存分布の比較を行います。\n",
        "\n",
        "ここでは、R データセットリポジトリで利用可能な [flchain](https://www.rdocumentation.org/packages/survival/versions/2.44-1.1/topics/flchain) [(Free Light Chain, 免疫グロブリン遊離L鎖)](http://www.eiken.co.jp/modern_media/backnumber/pdf/MM1209_02.pdf) に関する調査研究のデータを使って`SurvfuncRight` オブジェクトを生成ます。女性の被検者のデータだけを用いて生存曲線を推定してみましょう。"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "RBrZFweaMGbU",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Install required modules\n",
        "# !pip install statsmodels\n",
        "# !pip install pandas-profiling"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "8-RDFTTmN2g-",
        "colab_type": "text"
      },
      "source": [
        "まずはデータを読み込み、カラムと先頭のいくつかのデータを表示させましょう。\n",
        "データは 7874 名の被検者の情報を含んでいます。各カラムの意味は次のとおりです。\n",
        "\n",
        "- age : 年齢\n",
        "- sex : 性別\n",
        "- sample.yr : 血液を採取した年\n",
        "- kappa : 遊離グロブリンL鎖のうち、$\\kappa$タイプの割合\n",
        "- lambda : 遊離グロブリンL鎖のうち、$\\lambda$タイプの割合\n",
        "- flc.grp : 被検者の遊離グロブリンL鎖のグループ (もととなった研究で使われたもの)\n",
        "- creatinine : 血清中クレアチニン\n",
        "- mgus : 単クローン性免疫グロブリン血症 (monoclonal immunoglobulinemia : M蛋白血症) の罹患 (1 の場合発症 )\n",
        "- futime : 登録から死亡までの日数 (死亡日にサンプルが取得された被検者が3名いることに注意)\n",
        "- death : 0 は最後に連絡をとったときに生存していたことを、 1 は死亡していたことを示します。\n",
        "- chapter : 主な死因。分類は International Code of Diseases ICD-9 の分類に基づく。"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "_0LgD4UOMPTS",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import statsmodels.api as sm\n",
        "\n",
        "data = sm.datasets.get_rdataset(\"flchain\", \"survival\").data\n",
        "df = data.loc[data.sex == \"F\", :]\n",
        "df.head()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "nJDSShn8L56E",
        "colab_type": "text"
      },
      "source": [
        "もう少しデータの様子について理解を試みましょう。このステップを探索的データ解析（Exploratory data analysis) と呼びます。\n",
        "\n",
        "本来であれば、データについて丁寧に観察するステップですが、ここではいくつかの簡便な方法により行います。\n",
        "\n",
        "まずは特徴量の追加を行います。正常な人では $\\kappa$ と $\\lambda$ の比が一定である一方、M蛋白血症に罹患した場合には比が異常値を示すことが知られています。比を表すカラムを追加してみましょう。"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "JB7qvAOmV_GO",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "df = df.assign(kappa_lambda_ratio=df[\"kappa\"] / df['lambda'])\n",
        "df.head()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "0IDj20MsV3Rw",
        "colab_type": "text"
      },
      "source": [
        "次に、基礎統計量と呼ばれる平均値や標準偏差、最小値、最大値、中央値を確認します。"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "w949hQg5MnNs",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "df.describe()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "5xz8pXYlNrIY",
        "colab_type": "text"
      },
      "source": [
        "ここから、次のことがわかります。\n",
        "\n",
        "- データセットには 4350 名の女性が含まれている\n",
        "- 年齢は 50 歳以上と、比較的高齢な被検者が多い\n",
        "- `creatinine` には欠損がある\n",
        "\n",
        "各変数の分布を見てみましょう。"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "rlVduxwcTMh2",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "from pandas_profiling import ProfileReport\n",
        "\n",
        "ProfileReport(df)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "yaeMDCu_cZcf",
        "colab_type": "text"
      },
      "source": [
        "`kappa_lambda_ratio` と `mgus` が強く相関していることから、`kappa_lambda_ratio` がM蛋白血症の検査に使える可能性が考えられます。\n",
        "\n",
        "次に、生存曲線を作成していきましょう。"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "fzbyOni9MjZ8",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "sf = sm.SurvfuncRight(df[\"futime\"], df[\"death\"])"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Zsp_JGtAoxZ0",
        "colab_type": "text"
      },
      "source": [
        "得られた生存分布の主要な特徴は summary メソッドを呼ぶことで得られます。"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "5QKrF9sWMWpy",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "sf.summary()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "csnwLrN3o0rB",
        "colab_type": "text"
      },
      "source": [
        "生存分布で死亡率がある数値になる日数と、その信頼区間を求めることができます。この調査研究中で観測されたのは 30% の死亡のみのため、求められるのは死亡確率 0.3 までのみです。"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "YC5xXyhUMZun",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "sf.quantile(0.2)\n",
        "sf.quantile_ci(0.2)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "XbSaj3v4o3wn",
        "colab_type": "text"
      },
      "source": [
        "生存関数を描画するために、 `plot` メソッドを呼び出します。グラフの読み取り方は [Survival function - Wikipedia](https://en.wikipedia.org/wiki/Survival_function) が詳しいです。"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "urkZcWP5Mfk5",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "%config InlineBackend.figure_format = 'retina'\n",
        "import matplotlib as mpl\n",
        "import seaborn as sns\n",
        "\n",
        "sns.set_style(\"whitegrid\")\n",
        "mpl.rcParams['figure.figsize'] = [8.0, 8.0]\n",
        "\n",
        "sf.plot()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "95t-eDMMo9RB",
        "colab_type": "text"
      },
      "source": [
        "このデータセットには打ち切りが大量に含まれているため、打ち切りを除いて描画したい場合もあるでしょう。\n",
        "\n",
        "訳注 : Colab 上だともともと打ち切りが表示されていないのでなんのことかわからなくなってしまうのですが、純粋な Jupyter Notebook 上では次の画像が表示されます。ここでは純粋な Jupyter Notebook 上で打ち切りを表す縦ラインを消す方法について述べています。\n",
        "\n",
        "![](https://user-images.githubusercontent.com/30427632/71105723-ef8cca00-21b5-11ea-9499-730bc983f6b1.png)"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "3tTXppJNM_ns",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "fig = sf.plot()\n",
        "ax = fig.get_axes()[0]\n",
        "pt = ax.get_lines()[1]\n",
        "pt.set_visible(True)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "9cPFQihUpIPC",
        "colab_type": "text"
      },
      "source": [
        "95% 信頼区間を追加することもできます。通常、信頼区間の幅は分布の中央の部分にのみプロットされます。\n",
        "\n",
        "訳注 : コメントしている箇所を実行することで上記の文章の意味がわかります。"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "JoGeD7RaNLOU",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import matplotlib.pyplot as plt\n",
        "cycle = plt.rcParams['axes.prop_cycle'].by_key()['color']\n",
        "\n",
        "fig = sf.plot()\n",
        "lcb, ucb = sf.simultaneous_cb()\n",
        "ax = fig.get_axes()[0]\n",
        "ax.fill_between(sf.surv_times, lcb, ucb, color=cycle[0], alpha=0.2)\n",
        "ax.set_xlim(365, 365*10)\n",
        "ax.set_ylim(0.7, 1)\n",
        "# ax.set_xlim(0, 5500)\n",
        "# ax.set_ylim(0, 1)\n",
        "ax.set_ylabel(\"Proportion alive\")\n",
        "ax.set_xlabel(\"Days since enrollment\")\n",
        "fig"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Mwjx8NQwpRwD",
        "colab_type": "text"
      },
      "source": [
        "一つのグラフ上に2つのグループの生存曲線 (男性と女性) を書くこともできます。"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "JWqY3B_RNPTK",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "gb = data.groupby(\"sex\")\n",
        "# ax = plt.axes()\n",
        "sexes = []\n",
        "cycle = plt.rcParams['axes.prop_cycle'].by_key()['color']\n",
        "\n",
        "for index, g in enumerate(gb):\n",
        "    sexes.append(g[0])\n",
        "    sf = sm.SurvfuncRight(g[1][\"futime\"], g[1][\"death\"])\n",
        "    fig = sf.plot() if index == 0 else sf.plot(ax)\n",
        "    lcb, ucb = sf.simultaneous_cb()\n",
        "    ax = fig.get_axes()[0]\n",
        "    ax.fill_between(sf.surv_times, lcb, ucb, color=cycle[index], alpha=0.2)\n",
        "\n",
        "li = ax.get_lines()\n",
        "li[1].set_visible(False)\n",
        "li[3].set_visible(False)\n",
        "plt.figlegend((li[0], li[2]), sexes, \"center right\")\n",
        "ax.set_ylabel(\"Proportion alive\")\n",
        "ax.set_ylim(0.6, 1)\n",
        "ax.set_xlabel(\"Days since enrollment\")\n",
        "ax.set_xlim(365, 365*10)\n",
        "fig"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "2j56adwBpZUY",
        "colab_type": "text"
      },
      "source": [
        "`survdiff`を用いることで2つの生存分布について比較することもできます。`survdiff` はいくつかの標準的なノンパラメトリックな手法が実装されています。デフォルトはログランク検定です。"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "6CBRrJRTNUKM",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "stat, pv = sm.duration.survdiff(data.futime, data.death, data.sex)\n",
        "\n",
        "# stat : χ二乗検定の統計量\n",
        "# pv : p value (p値) 信頼水準が 95 % の場合 0.05 よりも小さければ有意\n",
        "stat, pv"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Td832x5Tpco_",
        "colab_type": "text"
      },
      "source": [
        "survdiff では次のようにして他の検定手法を用いることもできます。"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "FTvEbDLeR9j-",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Fleming-Harrington with p=1, i.e. weight by pooled survival time\n",
        "stat, pv = sm.duration.survdiff(data.futime, data.death, data.sex, weight_type='fh', fh_p=1)\n",
        "\n",
        "# Gehan-Breslow, weight by number at risk\n",
        "stat, pv = sm.duration.survdiff(data.futime, data.death, data.sex, weight_type='gb')\n",
        "\n",
        "# Tarone-Ware, weight by the square root of the number at risk\n",
        "stat, pv = sm.duration.survdiff(data.futime, data.death, data.sex, weight_type='tw')"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ifbX5Zs2poHz",
        "colab_type": "text"
      },
      "source": [
        "# 回帰\n",
        "\n",
        "ハザード回帰モデル (Cox 回帰モデル) は打ち切りのあるデータに対して適用可能な手法です。線形モデルや一般化線形モデルにおける共変量のように、イベントが発生するまでの時間にバリエーションをもたせることができます。これらのモデルは共変量の効果をハザード比として表します。ハザード比 (ある瞬間のイベント比) は共編料の値に応じて乗算される値のことです。\n",
        "\n",
        "訳注 : ハザード比はある瞬間における死亡率の比のことです。"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "HYwLErqZSIU7",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import statsmodels.api as sm\n",
        "import statsmodels.formula.api as smf\n",
        "\n",
        "data = sm.datasets.get_rdataset(\"flchain\", \"survival\").data\n",
        "del data[\"chapter\"]\n",
        "data = data.dropna()\n",
        "data = data.assign(lam=data[\"lambda\"])\n",
        "data = data.assign(female=(data[\"sex\"] == \"F\").astype(int))\n",
        "data = data.assign(year=data[\"sample.yr\"] - min(data[\"sample.yr\"]))\n",
        "status = data[\"death\"].values\n",
        "\n",
        "mod = smf.phreg(\"futime ~ age + female + creatinine + \"\n",
        "                \"np.sqrt(kappa) + np.sqrt(lam) + year + mgus\",\n",
        "                data, status=status, ties=\"efron\")\n",
        "rslt = mod.fit()\n",
        "rslt.summary()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Z1jh567aqA-B",
        "colab_type": "text"
      },
      "source": [
        "より詳細な例については [Examples](https://www.statsmodels.org/dev/examples/index.html#statsmodels-examples) を参照してください。\n",
        "\n",
        "Wiki にもノートブックのサンプルコードがいくつかあります。 [Wiki notebooks for PHReg and Survival Analysis](https://github.com/statsmodels/statsmodels/wiki/Examples#survival-analysis)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "B6aTLkJNqhVX",
        "colab_type": "text"
      },
      "source": [
        "# License\n",
        "\n",
        "Copyright (C) 2006, Jonathan E. Taylor\n",
        "All rights reserved.\n",
        "\n",
        "Copyright (c) 2006-2008 Scipy Developers.\n",
        "All rights reserved.\n",
        "\n",
        "Copyright (c) 2009-2018 statsmodels Developers.\n",
        "All rights reserved.\n",
        "\n",
        "\n",
        "Redistribution and use in source and binary forms, with or without\n",
        "modification, are permitted provided that the following conditions are met:\n",
        "\n",
        "  a. Redistributions of source code must retain the above copyright notice,\n",
        "     this list of conditions and the following disclaimer.\n",
        "  b. Redistributions in binary form must reproduce the above copyright\n",
        "     notice, this list of conditions and the following disclaimer in the\n",
        "     documentation and/or other materials provided with the distribution.\n",
        "  c. Neither the name of statsmodels nor the names of its contributors\n",
        "     may be used to endorse or promote products derived from this software\n",
        "     without specific prior written permission.\n",
        "\n",
        "\n",
        "THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS \"AS IS\"\n",
        "AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE\n",
        "IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE\n",
        "ARE DISCLAIMED. IN NO EVENT SHALL STATSMODELS OR CONTRIBUTORS BE LIABLE FOR\n",
        "ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL\n",
        "DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR\n",
        "SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER\n",
        "CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT\n",
        "LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY\n",
        "OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH\n",
        "DAMAGE."
      ]
    }
  ]
}