.. _survivalrst:

=================
Analyse de survie
=================


.. only:: html

    **Links:** :download:`notebook <survival.ipynb>`, :downloadlink:`html <survival2html.html>`, :download:`PDF <survival.pdf>`, :download:`python <survival.py>`, :downloadlink:`slides <survival.slides.html>`, :githublink:`GitHub|_doc/notebooks/ml/survival.ipynb|*`


.. code:: ipython3

    from jyquickhelper import add_notebook_menu
    add_notebook_menu()






.. contents::
    :local:





Quelques données
----------------

On récupère les données disponibles sur *open.data.gouv.fr* `Données
hospitalières relatives à l’épidémie de
COVID-19 <https://www.data.gouv.fr/fr/datasets/donnees-hospitalieres-relatives-a-lepidemie-de-covid-19/>`__.
Ces données ne permettent pas de construire la courbe de
`Kaplan-Meier <https://fr.wikipedia.org/wiki/Estimateur_de_Kaplan-Meier>`__.
On sait combien de personnes rentrent et sortent chaque jour mais on ne
sait pas quand une personne qui sort un 1er avril est entrée.

.. code:: ipython3

    import numpy.random as rnd
    
    import pandas
    df = pandas.read_csv("https://www.data.gouv.fr/en/datasets/r/63352e38-d353-4b54-bfd1-f1b3ee1cabd7", sep=";")
    gr = df[["jour", "rad", "dc"]].groupby(["jour"]).sum()
    diff = gr.diff().reset_index(drop=False)
    diff.head()






.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>jour</th>
          <th>rad</th>
          <th>dc</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>2020-03-18</td>
          <td>NaN</td>
          <td>NaN</td>
        </tr>
        <tr>
          <th>1</th>
          <td>2020-03-19</td>
          <td>695.0</td>
          <td>207.0</td>
        </tr>
        <tr>
          <th>2</th>
          <td>2020-03-20</td>
          <td>806.0</td>
          <td>248.0</td>
        </tr>
        <tr>
          <th>3</th>
          <td>2020-03-21</td>
          <td>452.0</td>
          <td>151.0</td>
        </tr>
        <tr>
          <th>4</th>
          <td>2020-03-22</td>
          <td>608.0</td>
          <td>210.0</td>
        </tr>
      </tbody>
    </table>
    </div>



.. code:: ipython3

    
    def donnees_artificielles(hosp, mu=14, nu=21):
        dt = pandas.to_datetime(hosp['jour'])
        res = []
        for i in range(hosp.shape[0]):
            date = dt[i].dayofyear
            h = hosp.iloc[i, 1]
            delay = rnd.exponential(mu, int(h))
            for j in range(delay.shape[0]):
                res.append([date - int(delay[j]), date, 1])
            h = hosp.iloc[i, 2]
            delay = rnd.exponential(nu, int(h))
            for j in range(delay.shape[0]):
                res.append([date - int(delay[j]), date , 0])
        return pandas.DataFrame(res, columns=["entree", "sortie", "issue"])
    
    
    data = donnees_artificielles(diff[1:].reset_index(drop=True)).sort_values('entree')
    data.head()






.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>entree</th>
          <th>sortie</th>
          <th>issue</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>518488</th>
          <td>-200</td>
          <td>19</td>
          <td>0</td>
        </tr>
        <tr>
          <th>541408</th>
          <td>-192</td>
          <td>27</td>
          <td>0</td>
        </tr>
        <tr>
          <th>476735</th>
          <td>-187</td>
          <td>2</td>
          <td>0</td>
        </tr>
        <tr>
          <th>587013</th>
          <td>-185</td>
          <td>42</td>
          <td>0</td>
        </tr>
        <tr>
          <th>476057</th>
          <td>-180</td>
          <td>1</td>
          <td>0</td>
        </tr>
      </tbody>
    </table>
    </div>



Chaque ligne est une personne, ``entree`` est le jour d’entrée à
l’hôpital, ``sortie`` celui de la sortie, ``issue``, 0 pour décès, 1
pour en vie.

.. code:: ipython3

    data.describe()






.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>entree</th>
          <th>sortie</th>
          <th>issue</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>count</th>
          <td>624130.000000</td>
          <td>624130.000000</td>
          <td>624130.000000</td>
        </tr>
        <tr>
          <th>mean</th>
          <td>169.704510</td>
          <td>184.532815</td>
          <td>0.806729</td>
        </tr>
        <tr>
          <th>std</th>
          <td>125.420957</td>
          <td>124.343186</td>
          <td>0.394864</td>
        </tr>
        <tr>
          <th>min</th>
          <td>-200.000000</td>
          <td>1.000000</td>
          <td>0.000000</td>
        </tr>
        <tr>
          <th>25%</th>
          <td>53.000000</td>
          <td>84.000000</td>
          <td>1.000000</td>
        </tr>
        <tr>
          <th>50%</th>
          <td>133.000000</td>
          <td>144.000000</td>
          <td>1.000000</td>
        </tr>
        <tr>
          <th>75%</th>
          <td>301.000000</td>
          <td>315.000000</td>
          <td>1.000000</td>
        </tr>
        <tr>
          <th>max</th>
          <td>366.000000</td>
          <td>366.000000</td>
          <td>1.000000</td>
        </tr>
      </tbody>
    </table>
    </div>



Il y a environ 80% de survie dans ces données.

.. code:: ipython3

    import numpy
    duree = data.sortie - data.entree
    deces = (data.issue == 0).astype(numpy.int32)

.. code:: ipython3

    import numpy
    import matplotlib.pyplot as plt
    from lifelines import KaplanMeierFitter
    fig, ax = plt.subplots(1, 1, figsize=(10, 4))
    kmf = KaplanMeierFitter()
    kmf.fit(duree, deces)
    kmf.plot(ax=ax)
    ax.legend();



.. image:: survival_9_0.png


Régression de Cox
-----------------

On reprend les données artificiellement générées et on ajoute une
variable identique à la durée plus un bruit mais quasi nul

.. code:: ipython3

    import pandas
    data_simple = pandas.DataFrame({'duree': duree, 'deces': deces,
                                    'X1': duree * 0.57 * deces + numpy.random.randn(duree.shape[0]),
                                    'X2': duree * (-0.57) * deces + numpy.random.randn(duree.shape[0])})
    data_simple.head()






.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>duree</th>
          <th>deces</th>
          <th>X1</th>
          <th>X2</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>518488</th>
          <td>219</td>
          <td>1</td>
          <td>125.653230</td>
          <td>-125.666662</td>
        </tr>
        <tr>
          <th>541408</th>
          <td>219</td>
          <td>1</td>
          <td>126.006024</td>
          <td>-125.327549</td>
        </tr>
        <tr>
          <th>476735</th>
          <td>189</td>
          <td>1</td>
          <td>107.920779</td>
          <td>-108.358230</td>
        </tr>
        <tr>
          <th>587013</th>
          <td>227</td>
          <td>1</td>
          <td>129.788930</td>
          <td>-130.045019</td>
        </tr>
        <tr>
          <th>476057</th>
          <td>181</td>
          <td>1</td>
          <td>103.642440</td>
          <td>-103.793008</td>
        </tr>
      </tbody>
    </table>
    </div>



.. code:: ipython3

    from sklearn.model_selection import train_test_split
    data_train, data_test = train_test_split(data_simple, test_size=0.8)

.. code:: ipython3

    from lifelines.fitters.coxph_fitter import CoxPHFitter
    cox = CoxPHFitter()
    cox.fit(data_train[['duree', 'deces', 'X1']], duration_col="duree", event_col="deces",
            show_progress=True)


.. parsed-literal::
    Iteration 1: norm_delta = 0.13943, step_size = 0.9000, log_lik = -250658.36250, newton_decrement = 889.93933, seconds_since_start = 0.0
    Iteration 2: norm_delta = 0.00660, step_size = 0.9000, log_lik = -249862.37270, newton_decrement = 2.81312, seconds_since_start = 0.0
    Iteration 3: norm_delta = 0.00073, step_size = 0.9000, log_lik = -249859.57376, newton_decrement = 0.03357, seconds_since_start = 0.1
    Iteration 4: norm_delta = 0.00000, step_size = 1.0000, log_lik = -249859.54017, newton_decrement = 0.00000, seconds_since_start = 0.1
    Convergence success after 4 iterations.




.. parsed-literal::

    <lifelines.CoxPHFitter: fitted with 124826 total observations, 100754 right-censored observations>



.. code:: ipython3

    cox.print_summary()





.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <tbody>
        <tr>
          <th>model</th>
          <td>lifelines.CoxPHFitter</td>
        </tr>
        <tr>
          <th>duration col</th>
          <td>'duree'</td>
        </tr>
        <tr>
          <th>event col</th>
          <td>'deces'</td>
        </tr>
        <tr>
          <th>baseline estimation</th>
          <td>breslow</td>
        </tr>
        <tr>
          <th>number of observations</th>
          <td>124826</td>
        </tr>
        <tr>
          <th>number of events observed</th>
          <td>24072</td>
        </tr>
        <tr>
          <th>partial log-likelihood</th>
          <td>-249859.54</td>
        </tr>
        <tr>
          <th>time fit was run</th>
          <td>2021-02-24 23:48:57 UTC</td>
        </tr>
      </tbody>
    </table>
    </div><table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th style="min-width: 12px;"></th>
          <th style="min-width: 12px;">coef</th>
          <th style="min-width: 12px;">exp(coef)</th>
          <th style="min-width: 12px;">se(coef)</th>
          <th style="min-width: 12px;">coef lower 95%</th>
          <th style="min-width: 12px;">coef upper 95%</th>
          <th style="min-width: 12px;">exp(coef) lower 95%</th>
          <th style="min-width: 12px;">exp(coef) upper 95%</th>
          <th style="min-width: 12px;">z</th>
          <th style="min-width: 12px;">p</th>
          <th style="min-width: 12px;">-log2(p)</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>X1</th>
          <td>0.02</td>
          <td>1.02</td>
          <td>0.00</td>
          <td>0.02</td>
          <td>0.02</td>
          <td>1.02</td>
          <td>1.02</td>
          <td>42.23</td>
          <td>&lt;0.005</td>
          <td>inf</td>
        </tr>
      </tbody>
    </table><br><div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <tbody>
        <tr>
          <th>Concordance</th>
          <td>0.69</td>
        </tr>
        <tr>
          <th>Partial AIC</th>
          <td>499721.08</td>
        </tr>
        <tr>
          <th>log-likelihood ratio test</th>
          <td>1597.64 on 1 df</td>
        </tr>
        <tr>
          <th>-log2(p) of ll-ratio test</th>
          <td>inf</td>
        </tr>
      </tbody>
    </table>
    </div>


.. code:: ipython3

    cox2 = CoxPHFitter()
    cox2.fit(data_train[['duree', 'deces', 'X2']], duration_col="duree", event_col="deces",
            show_progress=True)
    cox2.print_summary()


.. parsed-literal::
    Iteration 1: norm_delta = 0.13946, step_size = 0.9000, log_lik = -250658.36250, newton_decrement = 888.92089, seconds_since_start = 0.0
    Iteration 2: norm_delta = 0.00667, step_size = 0.9000, log_lik = -249863.61089, newton_decrement = 2.86434, seconds_since_start = 0.0
    Iteration 3: norm_delta = 0.00074, step_size = 0.9000, log_lik = -249860.76079, newton_decrement = 0.03426, seconds_since_start = 0.1
    Iteration 4: norm_delta = 0.00000, step_size = 1.0000, log_lik = -249860.72650, newton_decrement = 0.00000, seconds_since_start = 0.1
    Convergence success after 4 iterations.





.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <tbody>
        <tr>
          <th>model</th>
          <td>lifelines.CoxPHFitter</td>
        </tr>
        <tr>
          <th>duration col</th>
          <td>'duree'</td>
        </tr>
        <tr>
          <th>event col</th>
          <td>'deces'</td>
        </tr>
        <tr>
          <th>baseline estimation</th>
          <td>breslow</td>
        </tr>
        <tr>
          <th>number of observations</th>
          <td>124826</td>
        </tr>
        <tr>
          <th>number of events observed</th>
          <td>24072</td>
        </tr>
        <tr>
          <th>partial log-likelihood</th>
          <td>-249860.73</td>
        </tr>
        <tr>
          <th>time fit was run</th>
          <td>2021-02-24 23:48:59 UTC</td>
        </tr>
      </tbody>
    </table>
    </div><table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th style="min-width: 12px;"></th>
          <th style="min-width: 12px;">coef</th>
          <th style="min-width: 12px;">exp(coef)</th>
          <th style="min-width: 12px;">se(coef)</th>
          <th style="min-width: 12px;">coef lower 95%</th>
          <th style="min-width: 12px;">coef upper 95%</th>
          <th style="min-width: 12px;">exp(coef) lower 95%</th>
          <th style="min-width: 12px;">exp(coef) upper 95%</th>
          <th style="min-width: 12px;">z</th>
          <th style="min-width: 12px;">p</th>
          <th style="min-width: 12px;">-log2(p)</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>X2</th>
          <td>-0.02</td>
          <td>0.98</td>
          <td>0.00</td>
          <td>-0.02</td>
          <td>-0.02</td>
          <td>0.98</td>
          <td>0.98</td>
          <td>-42.21</td>
          <td>&lt;0.005</td>
          <td>inf</td>
        </tr>
      </tbody>
    </table><br><div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <tbody>
        <tr>
          <th>Concordance</th>
          <td>0.69</td>
        </tr>
        <tr>
          <th>Partial AIC</th>
          <td>499723.45</td>
        </tr>
        <tr>
          <th>log-likelihood ratio test</th>
          <td>1595.27 on 1 df</td>
        </tr>
        <tr>
          <th>-log2(p) of ll-ratio test</th>
          <td>inf</td>
        </tr>
      </tbody>
    </table>
    </div>


.. code:: ipython3

    cox.predict_cumulative_hazard(data_test[:5])






.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>621725</th>
          <th>110352</th>
          <th>139986</th>
          <th>72623</th>
          <th>248121</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0.0</th>
          <td>0.008806</td>
          <td>0.008673</td>
          <td>0.008543</td>
          <td>0.008717</td>
          <td>0.008661</td>
        </tr>
        <tr>
          <th>1.0</th>
          <td>0.018218</td>
          <td>0.017942</td>
          <td>0.017673</td>
          <td>0.018035</td>
          <td>0.017918</td>
        </tr>
        <tr>
          <th>2.0</th>
          <td>0.026735</td>
          <td>0.026330</td>
          <td>0.025936</td>
          <td>0.026466</td>
          <td>0.026295</td>
        </tr>
        <tr>
          <th>3.0</th>
          <td>0.036002</td>
          <td>0.035457</td>
          <td>0.034926</td>
          <td>0.035640</td>
          <td>0.035409</td>
        </tr>
        <tr>
          <th>4.0</th>
          <td>0.045543</td>
          <td>0.044854</td>
          <td>0.044182</td>
          <td>0.045085</td>
          <td>0.044793</td>
        </tr>
        <tr>
          <th>...</th>
          <td>...</td>
          <td>...</td>
          <td>...</td>
          <td>...</td>
          <td>...</td>
        </tr>
        <tr>
          <th>184.0</th>
          <td>2.055736</td>
          <td>2.024623</td>
          <td>1.994277</td>
          <td>2.035070</td>
          <td>2.021872</td>
        </tr>
        <tr>
          <th>189.0</th>
          <td>2.092002</td>
          <td>2.060340</td>
          <td>2.029459</td>
          <td>2.070971</td>
          <td>2.057541</td>
        </tr>
        <tr>
          <th>197.0</th>
          <td>2.138687</td>
          <td>2.106318</td>
          <td>2.074747</td>
          <td>2.117186</td>
          <td>2.103457</td>
        </tr>
        <tr>
          <th>201.0</th>
          <td>2.205513</td>
          <td>2.172133</td>
          <td>2.139576</td>
          <td>2.183341</td>
          <td>2.169182</td>
        </tr>
        <tr>
          <th>217.0</th>
          <td>2.330629</td>
          <td>2.295356</td>
          <td>2.260952</td>
          <td>2.307199</td>
          <td>2.292237</td>
        </tr>
      </tbody>
    </table>
    <p>165 rows × 5 columns</p>
    </div>



.. code:: ipython3

    cox.predict_survival_function(data_test[:5])






.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>621725</th>
          <th>110352</th>
          <th>139986</th>
          <th>72623</th>
          <th>248121</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0.0</th>
          <td>0.991233</td>
          <td>0.991365</td>
          <td>0.991494</td>
          <td>0.991321</td>
          <td>0.991377</td>
        </tr>
        <tr>
          <th>1.0</th>
          <td>0.981947</td>
          <td>0.982218</td>
          <td>0.982482</td>
          <td>0.982127</td>
          <td>0.982242</td>
        </tr>
        <tr>
          <th>2.0</th>
          <td>0.973619</td>
          <td>0.974013</td>
          <td>0.974398</td>
          <td>0.973881</td>
          <td>0.974048</td>
        </tr>
        <tr>
          <th>3.0</th>
          <td>0.964638</td>
          <td>0.965164</td>
          <td>0.965677</td>
          <td>0.964988</td>
          <td>0.965211</td>
        </tr>
        <tr>
          <th>4.0</th>
          <td>0.955478</td>
          <td>0.956137</td>
          <td>0.956780</td>
          <td>0.955916</td>
          <td>0.956196</td>
        </tr>
        <tr>
          <th>...</th>
          <td>...</td>
          <td>...</td>
          <td>...</td>
          <td>...</td>
          <td>...</td>
        </tr>
        <tr>
          <th>184.0</th>
          <td>0.127999</td>
          <td>0.132044</td>
          <td>0.136112</td>
          <td>0.130671</td>
          <td>0.132407</td>
        </tr>
        <tr>
          <th>189.0</th>
          <td>0.123440</td>
          <td>0.127411</td>
          <td>0.131407</td>
          <td>0.126063</td>
          <td>0.127768</td>
        </tr>
        <tr>
          <th>197.0</th>
          <td>0.117809</td>
          <td>0.121685</td>
          <td>0.125588</td>
          <td>0.120370</td>
          <td>0.122034</td>
        </tr>
        <tr>
          <th>201.0</th>
          <td>0.110194</td>
          <td>0.113934</td>
          <td>0.117705</td>
          <td>0.112664</td>
          <td>0.114271</td>
        </tr>
        <tr>
          <th>217.0</th>
          <td>0.097235</td>
          <td>0.100726</td>
          <td>0.104251</td>
          <td>0.099540</td>
          <td>0.101040</td>
        </tr>
      </tbody>
    </table>
    <p>165 rows × 5 columns</p>
    </div>