Bartosz KostkaSurvival Analysis
# -*- coding: utf-8 -*-
%pylab inline
import pandas as pd
import seaborn as sns
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.
dane = pd.read_csv('pharynx.csv', names=["INST", "SEX", "TX", "GRADE", "AGE", "COND", "SITE", "T_STAGE", "N_STAGE", "ENTRY_DT", "STATUS", "TIME"])
dane.head()
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.
dane["AGE"].hist()
Sprawdzamy dane statystyczne dot. czasu przeżycia, bądź wycofania się z badania.
dst = dane.groupby('STATUS')
dst["TIME"].describe()
dane["TIME"].hist(by=dane["STATUS"])
Przy pomocy estymatora Kaplana Meiera estymujemy funkcję przeżycia.
from lifelines import KaplanMeierFitter
kmf = KaplanMeierFitter()
kmf.fit(dane["TIME"], event_observed=dane["STATUS"])
kmf.plot()
W badaniu chcieliśmy porównać sposoby leczenia, co też czynimy.
dtx = dane.groupby(['STATUS', 'TX'])
dtx["TIME"].describe()
danett = dane[dane["TX"]==1]
danest = dane[dane["TX"]==2]
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")
Tworzymy model używając regresji liniowej (nie wiem jeszcze, która wersja jest ładniejsza).
from sklearn.linear_model import LinearRegression
from sklearn.feature_selection import chi2
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)
from lifelines import AalenAdditiveFitter
aaf = AalenAdditiveFitter(coef_penalizer=1.0, fit_intercept=True)
aaf.fit(dane, 'TIME', event_col='STATUS')
aaf.cumulative_hazards_.head()
from lifelines import CoxPHFitter
cf = CoxPHFitter()
cf.fit(dane, 'TIME', event_col='STATUS')
cf.print_summary()