In [13]:
# -*- coding: utf-8 -*-
%pylab inline
import pandas as pd
import seaborn as sns
Populating the interactive namespace from numpy and matplotlib

Poniższe dane przedstawiają dane dużej próby klinicznej, która zorganizowana była przez grupę leczenia onkologii przed radioterapię w Stanach Zjednoczonych.

Badania dotyczyły pacjentów z rakiem kolczystokomórkowym skóry w 15 miejscach w ustach i gardle człowika z 16 ośrodków (poniższe dane reprezentują tylko miejsca w części ustnej gardła, z 6 największych instytucji).

Pacjenci poddani badaniu zostali losowo podzieleni na jedną z dwóch grup leczenia: radioterapię oraz radioterapię razem z chemioterapią. Jednym z celów badań było porównanie obu tych sposobów leczenia wobec śmiertelności pacjentów.

In [3]:
dane = pd.read_csv('pharynx.csv', names=["INST", "SEX", "TX", "GRADE", "AGE", "COND", "SITE", "T_STAGE", "N_STAGE", "ENTRY_DT", "STATUS", "TIME"])
dane.head()
Out[3]:
INST SEX TX GRADE AGE COND SITE T_STAGE N_STAGE ENTRY_DT STATUS TIME
1 2 2 1 1 51 1 2 3 1 2468 1 631
2 2 1 2 1 65 1 4 2 3 2968 1 270
3 2 1 1 2 64 2 1 3 3 3368 1 327
4 2 1 1 1 73 1 1 4 0 5768 1 243
5 5 1 2 2 64 1 1 4 3 9568 1 916

Opis zmiennych:

  • INST - numer instytucji biorącej udział w badniu
  • SEX - płeć 1=mężczyzna, 2=kobieta
  • TX - leczenie 1=standardowe, 2=testowe
  • GRADE - typ raka 1=złośliwy, 2=średni, 3=niezłośliwy, 9=brak danych
  • AGE - wiek w czasie diagnozy
  • COND - kondycja pacjenta: 1=brak niepełnosprawności, 2=ograniczona sprawność, 3=potrzebuje pomocy w codziennych czynnościach, 4=nie rusza się z łóżka, 9=brak danych
  • SITE - miejsce wykrycia raka: 1=faucial arch, 2=tonsillar fossa, 3=posterior pillar, 4=pharyngeal tongue, 5=posterior wall
  • T_STAGE - 1=pierwotny guz mniejszy niż 2cm w najdłuższej średnicy, 2=pierwotny guz między 2 a 4cm w najdłuższej średnicy z minimalnym przenikaniem w głębokości, 3=pierwotny guz mierzy więcej niż 4cm, 4=wielki, inwazyjny guz
  • N_STAGE - 0=brak medycznej ewidencji o przerzutach, 1=pojedynczy zajęty węzeł chłonny o średnicy 3cm lub mniejszy, nie operowany, 2=pojedynczy przypadek zajętego wezła chłonnego większego niż 3cm w średnicy, nie operowany, 3=wiele zajętych węzłów chłonnych lub operatowane
  • ENTRY_DT - data wstąpeinia do badania format: dddyy
  • STATUS - 0=ocenzurowane, 1=martwy
  • TIME - czas przeżycia od diagnozy

Poniższy histogram pokazuje histogram wieku ludzi, którzy poddali się badaniu.

In [25]:
dane["AGE"].hist()
Out[25]:
<matplotlib.axes._subplots.AxesSubplot at 0x7f137e6defd0>

Sprawdzamy dane statystyczne dot. czasu przeżycia, bądź wycofania się z badania.

In [109]:
dst = dane.groupby('STATUS')
dst["TIME"].describe()
Out[109]:
STATUS       
0       count      53.000000
        mean     1013.716981
        std       408.369006
        min        90.000000
        25%       733.000000
        50%       943.000000
        75%      1317.000000
        max      1823.000000
1       count     142.000000
        mean      388.908451
        std       269.939497
        min        11.000000
        25%       205.750000
        50%       325.500000
        75%       525.750000
        max      1574.000000
dtype: float64
In [110]:
dane["TIME"].hist(by=dane["STATUS"])
Out[110]:
array([<matplotlib.axes._subplots.AxesSubplot object at 0x7fbb2c440310>,
       <matplotlib.axes._subplots.AxesSubplot object at 0x7fbb2c4de0d0>], dtype=object)

Przy pomocy estymatora Kaplana Meiera estymujemy funkcję przeżycia.

In [113]:
from lifelines import KaplanMeierFitter
kmf = KaplanMeierFitter()
kmf.fit(dane["TIME"], event_observed=dane["STATUS"])
kmf.plot()
Out[113]:
<matplotlib.axes._subplots.AxesSubplot at 0x7fbb2c399b10>

W badaniu chcieliśmy porównać sposoby leczenia, co też czynimy.

In [116]:
dtx = dane.groupby(['STATUS', 'TX'])
dtx["TIME"].describe()
Out[116]:
STATUS  TX       
0       1   count      27.000000
            mean     1020.814815
            std       359.267207
            min        90.000000
            25%       808.000000
            50%       943.000000
            75%      1278.500000
            max      1609.000000
        2   count      26.000000
            mean     1006.346154
            std       461.004290
            min       182.000000
            25%       618.000000
            50%       994.000000
            75%      1317.000000
            max      1823.000000
1       1   count      73.000000
            mean      443.931507
            std       316.968511
            min        81.000000
            25%       213.000000
            50%       363.000000
            75%       599.000000
            max      1574.000000
        2   count      69.000000
            mean      330.695652
            std       194.999646
            min        11.000000
            25%       192.000000
            50%       291.000000
            75%       446.000000
            max       916.000000
dtype: float64
In [117]:
danett = dane[dane["TX"]==1]
danest = dane[dane["TX"]==2]
In [118]:
kmf.fit(danest["TIME"], event_observed=danest["STATUS"], label="normal treatment")
ax = kmf.plot()
kmf.fit(danett["TIME"], event_observed=danett["STATUS"], label="test treatment")
kmf.plot(ax=ax, c="r")
Out[118]:
<matplotlib.axes._subplots.AxesSubplot at 0x7fbb2c1ba5d0>

Tworzymy model używając regresji liniowej (nie wiem jeszcze, która wersja jest ładniejsza).

In [58]:
from sklearn.linear_model import LinearRegression
from sklearn.feature_selection import chi2
In [93]:
model = LinearRegression()
ozn = ["INST", "SEX", "TX", "GRADE", "AGE", "COND", "SITE", "T_STAGE", "N_STAGE"]
X = vstack([dane[x] for x in ozn])
model.fit(X=X.T, y=dane["TIME"])
tests = chi2(X=X.T, y=dane["TIME"])
war = pd.DataFrame(vstack([model.coef_, tests[0], tests[1]]).T, columns=["beta", "chi2", "p-value"], index=ozn)
display(war)
beta chi2 p-value
INST 6.846629 136.631934 9.875441e-01
SEX 87.372344 25.473029 1.000000e+00
TX -68.345070 30.965517 1.000000e+00
GRADE 20.920558 63.768930 1.000000e+00
AGE 0.677060 369.817989 7.394772e-16
COND -105.408137 100.581784 9.999991e-01
SITE 16.383039 120.035011 9.995786e-01
T_STAGE -100.103428 38.525493 1.000000e+00
N_STAGE -61.271206 137.720000 9.851548e-01
In [6]:
from lifelines import AalenAdditiveFitter
In [11]:
aaf = AalenAdditiveFitter(coef_penalizer=1.0, fit_intercept=True)
aaf.fit(dane, 'TIME', event_col='STATUS')
aaf.cumulative_hazards_.head()
 [-----------------100%-----------------] 142 of 142 complete in 0.2 sec
Out[11]:
INST SEX TX GRADE AGE COND SITE T_STAGE N_STAGE ENTRY_DT baseline
11 0.000455 0.005074 0.018043 -0.012589 -0.000061 0.009409 0.010918 0.013836 0.006880 -2.405703e-07 -0.086005
15 0.003675 -0.005402 0.027258 -0.015880 -0.000093 0.012477 0.017400 0.021719 0.009615 -5.090018e-07 -0.127921
38 0.006789 0.010413 0.036698 -0.016270 -0.000048 0.014433 0.023921 0.020561 0.000907 -6.198836e-07 -0.164098
74 0.009123 0.003714 0.044493 -0.016419 -0.000525 0.025700 0.019189 0.025721 0.005818 -1.139921e-06 -0.161061
81 0.011971 -0.002002 0.035889 -0.022151 -0.000600 0.025274 0.016909 0.030665 0.004741 -6.442594e-07 -0.146757
In [4]:
from lifelines import CoxPHFitter
In [5]:
cf = CoxPHFitter()
cf.fit(dane, 'TIME', event_col='STATUS')
cf.print_summary() 
n=195, number of events=142

               coef  exp(coef)  se(coef)          z         p  lower 0.95  upper 0.95     
INST     -2.210e-02  9.781e-01 9.121e-02 -2.423e-01 8.086e-01  -2.009e-01   1.567e-01     
SEX      -8.488e-02  9.186e-01 8.783e-02 -9.664e-01 3.338e-01  -2.571e-01   8.730e-02     
TX        8.647e-02  1.090e+00 8.690e-02  9.951e-01 3.197e-01  -8.389e-02   2.568e-01     
GRADE    -9.875e-02  9.060e-01 1.000e-01 -9.873e-01 3.235e-01  -2.948e-01   9.734e-02     
AGE       1.927e-02  1.019e+00 1.013e-01  1.902e-01 8.492e-01  -1.793e-01   2.179e-01     
COND      2.176e-01  1.243e+00 6.417e-02  3.391e+00 6.955e-04   9.183e-02   3.434e-01  ***
SITE     -6.959e-02  9.328e-01 9.249e-02 -7.524e-01 4.518e-01  -2.509e-01   1.117e-01     
T_STAGE   2.259e-01  1.253e+00 9.449e-02  2.391e+00 1.681e-02   4.067e-02   4.112e-01    *
N_STAGE   2.185e-01  1.244e+00 9.254e-02  2.361e+00 1.822e-02   3.708e-02   3.999e-01    *
ENTRY_DT  1.258e-02  1.013e+00 9.363e-02  1.344e-01 8.931e-01  -1.710e-01   1.961e-01     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Concordance = 0.648
/home/kostka/anaconda2/lib/python2.7/site-packages/lifelines/fitters/coxph_fitter.py:285: FutureWarning: sort(columns=....) is deprecated, use sort_values(by=.....)
  df.sort(duration_col, inplace=True)

© Bartosz Kostka 2013-2017