Un titolo dice tante cose ma non tutto.

Il nome Survival Functions ci ricorda che si tratta di strumenti inizialmente sviluppati ed utilizzati nel campo della Medicina. Ma, come spero di chiarire con alcuni buoni esempi, in realtà possono essere applicati in tanti differenti campi, ed io ritengo che debbano rientrare tra gli strumenti padroneggiati dai Data Scientist di oggi.

Un primo esempio.

In medicina purtroppo accade che un intervento effettuato su un paziente possa non essere risolutivo. O che dopo l'intervento vi sia una probabilità significativa che il paziente sopravviva per un tempo non molto lungo. In questo caso il medico, e a volte il paziente o la sua famiglia, vogliono sapere quali sono le probabilità di sopravvivenza nel tempo (a 6 mesi, un anno, cinque anni, etc).

Possiamo schematizzare la situazione in questo modo:

1. Al momento t = 0 è effettuato l'intervento (es: un trapianto di un organo importante)
2. Vogliamo conoscere qual è la probabilità che il paziente sopravviva per più di un tempo t (per un valore arbitrario di t)

Il metodo che è stato utilizzato tante volte nel passato, per pervenire a modelli che consentano una tale previsione, è di effettuare uno studio controllato su un insieme il più possibile ampio di pazienti nella stessa situazione, registrando per ciascun paziente la data dell'intervento e, laddove possibile, la data di morte.

Per rendere più leggera ed ottimistica la lettura, anticipo che tali strumenti possono essere utilizzati per casi differenti dalla "sopravvivenza vera e propria". Ad esempio per prevedere quanto tempo passa tra la nascita del primo e del secondo figlio, quanto tempo passa tra la perdita del lavoro ed un nuovo impiego, etc. Per rimanere fedeli al titolo, continuerò ad eusare l'esempio originario

Non sempre è possibile registrare accuratamente i dati. Vi sono più ragioni possibili per cui tali dati non sono tutti disponibili e ciò conduce ad un aspetto, denominato "right censoring", che va considerato nel calcolo delle probabilità.

Nel seguito cercherò di chiarire i termini ed i concetti.

Survival Function.

In uno scenario quale quello brevemente delineato, si definisce Survival Function S(t) la funzione che, per un qualsivoglia valore di t, fornisce la probabilità che il paziente sopravviva per il tempo t. Ovvero, la probabilità che l'evento di morte avvenga ad un tempo T > t

S(t) = Pr(T > t)

La Survival Function per come la definiamo e stimiamo in questo articolo è relativa all'intera popolazione. Su questo punto vi tornerò alla fine dell'articolo.

Del tutto in generale possiamo subito elencare alcune delle proprietà matematiche che una generica Survival Function deve avere.

Deve essere:

0 =< S(t) <= 1

Inoltre:

S(0) = 1

E' sempre una funzione non crescente, ovvero

u < v implica S(u) >= S(v)

Questa proprietà deriva dal fatto che la Survival Function da una probabilità cumulativa (di sopravvivenza fino al tempo t) e non la probabilità puntuale che l'evento si verifichi al tempo t (quest'ultimo concetto è legato all'hazard function).

Infine deve tendere a zero

lim S(t) = 0 per t -> ∞

 

Ecco il grafico della Survival Function per un dataset reale (spiegherò poi da quali dati è stata ricavata):

Fig. 2

 

Come stimare la Survival Function.

La prima cosa che dobbiamo chiederci è: come deve essere fatto un dataset per poter stimare la Survival Function S(t), per ogni tempo t?

Dobbiamo ovviamente avere un numero sufficientemente grande di rilevazioni provenienti da una popolazione omogenea. Tali rilevazioni mi devono dire, per ciascun paziente, quanto tempo è passato tra l'evento iniziale (es: l'intervento, oppure la diagnosi) e la morte del paziente stesso.

Ovviamente, per i vari pazienti l'intervento (o l'evento di inizio, più in generale) non cade nello stesso giorno e, quindi, quello che conta realmente è il tempo trascorso tra l'evento iniziale (detto a volte anche birth, nascita) e l'evento finale (death).

Il dataset ideale, frutto e risultato dello studio ideale, contiene dunque un numero abbastanza grande di record del tipo:

id	time (mesi)
1	12
2	14
3	13

tab. 1

ove id è l'identificativo del paziente, che non serve a determinare la probabilità di sopravvivenza ma, soltanto, a controllare che non vi siano duplicati, e time che è il tempo trascorso tra inizio e fine (ad esempio, espresso in mesi).

A questo punto per stimare la S(t) dovremmo soltanto estrarre dalla tabella la lista ordinata dei tempi ti e contare, per ciascun valore di ti, la frazione di pazienti ancora vivi.

Ma vi è un problema: gli studi effettuati nel mondo reale mostrano di frequenza eccezioni. Ad esempio, e si spera che accada, molti pazienti possono essere sopravvissuti dopo il termine dello studio. Per questi pazienti l'intervallo tra inizio e fine quanto vale?

Questa osservazione porta ad introdurre il concetto di "censored data" e, in particolare, di "right censored", che dobbiamo definire bene e lo faremo nel capitolo seguente.  

Altra osservazione, più banale direi, è che, per forza di cose, la funzione la possiamo stimare soltanto in corrispondenza dei tempi contenuti nel dataset (nella tabella 1). Quindi per un insieme discreto di valori. Ciò complica se dobbiamo, per fare un esempio che poi si comprenderà meglio quando parleremo di Hazard Function, se dobbiamo calcolare la derivata della funzione.

Right Censoring: di cosa si tratta e perchè non va ignorato.

In uno studio statistico reale vi possono esser tante ragioni per cui i dati contenuti nel dataset di cui abbiamo parlato possono essere imprecisi o addirittura mancanti.

Proviamo ad elencarne alcuni:

1. Un paziente potrebbe abbandonare lo studio, ad esempio perchè cambia nazione e non è più possibile contattarlo.
2. Un paziente potrebbe morire per ragioni del tutto estranee allo studio in oggetto (es: un incidente stradale);
3. Lo studio ha un termine temporale e, per fortuna, numerosi pazienti sopravvivono oltre la fine dello studio.

Se teniamo conto di queste evenienze, il nostro dataset della tabella 1 deve diventare qualcosa di simile alla tabella 2:

id	time (mesi)	esito
1	12	1
2	14	1
3	13	0
4	16	1
5	11	0

tab. 2

Come interpretiamo i dati? 

Dunque: time ora diviene il tempo relativo all'ultimo evento noto. Esito può valere 0, 1, dove 1 indica la morte del paziente.

Se teniamo conto del fatto che lo studio ha una durata temporale finita ci rendiamo conto che le righe che contengono esito = 0 contengono dei dati mancanti (missing data). Per tali pazienti sappiamo con certezza che sono in vita al tempo time, ma non sappiamo cosa è accaduto dopo. Tali dati si chiamano "censored data". Inoltre, poichè l'informazione mancante è sulla destra della scala temporale si parla di "right censored data".

Come gestiamo i "right censored data" dal punto di vista della stima della Survival Function?

Vediamo: prima del time indicato nella riga abbiamo delle informazioni certe: il paziente è vivo e quindi, senza nessun dubbio, utilizziamo tale informazione per il computo della S(T) per t <= time. Ma per t > time non sappiamo nulla. IL paziente potrebbe essere morto subito dopo (dead immediately) oppure potrebbe, all'altro estremo, essere vivo fino ad oggi.

Se ci fermiamo a ragionarci su un attimo ci rendiamo conto che:

1. Se decidessimo di cancellare dal dataset tutte le righe "right censored", avremmo per molti t una sottostima di S(t)
2. In generale il valore di S(T) è compreso tra quanto otteniamo assumendo che tutti i pazienti "right censored" muoiono immediatamente dopo il time e quanto otteniamo assumendo che vivano indefinitivamente.

Con la formula di Kaplan-Meier possiamo ottenere una stima che tiene conto meglio dei right censored data e quindi fornisce una stima più accurata. 

La formula di Kaplan-Meier, uno stimatore non-parametrico.

Una tecnica per stimare la Survival Function è di utilizzare la formula di Kaplan-Meier, che ci fornisce uno stimatore non-parametrico.

 

S(t) = ∏ (1 - di/ni) 

ove:

• il prodotto è esteso a tutti i tempi ti <= t
• di è il numero di "pazienti morti al tempo ti"
• ni è "il numero di pazienti noti come vivi al tempo ti"

Osservando il prodotto che compare nella formula di Kaplan-Meier ricaviamo una proprietà interessante, che può essere molto utile per ridurre la complessità computazionale necessaria per stimare la S(t) su tutti i tempi della tabella 1.

Se indichiamo con tm e tm+1 due tempi consecutivi, risulta

S(tm+1) = S(tm) * (1 - dm/nm)

quindi, se dobbiamo calcolare tutta la sequenza dei tm, possiamo iterare calcolando e memorizzando gli S(tm) in un array. L'elemento successivo richiederà un solo prodotto, a partire dal precedente.

Qui riporto un frammento di codice Python che implementa il calcolo 

<inserire gist qui !!!>

Come ho detto, lo stimatore è non parametrico. La formula fornisce un algoritmo che non prevede parametri da ricavare da un fit sui dati. Per inciso vi sono altri stimatori che prevedono un tale fit e quindi sono detti, appunto, parametrici.

Applicazione al dataset AML.

AML è acronimo (in inglese) di Leucemia Mieloide Acuta. Una neoplasia grave che produce una rapida crescita del numero di globuli bianchi anomali.

E' una malattia che se non curata adeguatamente può produrre un esito fatale anche in poche settimane e, purtroppo, in pazienti anziani, che non possono sopportare la chemioterapia, comunque non permette una prognosi favorevole, con sopravvivenze di pochi mesi. Ovvio che su tale malattia si siano prodotti tanti studi volti a stimare la probabilità di sopravvivenza nel tempo.

AML è un dataset noto, analizzato in molti testi di medicina e fornito con il linguaggio R. Una descrizione accurata del dataset può essere trovata nella pagina Wikipedia: https://en.wikipedia.org/wiki/Survival_analysis#Example:_Acute_myelogenous_leukemia_survival_data

Ne riporto un estratto:

 

La prima colonna è un id.

Time è espresso in mesi. 

Cens vale zero se il dato è da considerare censored (il paziente  al tempo riportato è in remissione dei sintomi). UNo indica che si è avuta recrudescenza della malattia. 

Group sta ad indicare se per il paziente è stata effettuata o meno chemioterapia.

Il dataset contiene 24 record, possiamo utilizzarlo per stimare la survival function.

Nel mio repository github med4ainotes, nella directory corso2, potete trovare sia il dataset che un Notebook che calcola la Survival Function su tale dataset.

Ecco il grafico, che confronta il risultato ottenuto dal mio codice custom, confrontato con il risultato prodotto dal package lifeline, che calcola anche l'intervallo di confidenza.

 

(in colore arancione la curva ottenuta dal codice custom, in blu quella prodotta da lifelines).

Dalla curva si vede che la probabilità di non avere sintomi scende rapidamente a circa lo 0.1 in circa 45 mesi. La retta finale è dovuta alla scala ed al fatto che per t > 45 i dati sono pochissimi.

Ovviamente, il grafico è relativo al dataset che è abbastanza datato (1997), spero che oggi la prognosi sia più favorevole, per il progresso delle cure. 

Un'ulteriore considerazione: la colonna "group" permetterebbe di analizzare separatamente la probabilità nei due gruppi (pazienti con chemio mantenuta e non).

 

Lifelines, un package Python specializzato per la Survival Analysis.

Un package Python completo e ben fatto, che fornisce un ampio insieme di strumenti per la Survival Analysis, è lifelines.

La documentazione è alla URL: https://lifelines.readthedocs.io/en/latest/index.html

Il link al repository github, ove si possono trovare esempi: https://github.com/CamDavidsonPilon/lifelines/

Esempio concreto di applicazione ad un caso estraneo alla Medicina.

Le tecniche di Survival Analysis possono essere applicate con successo anche a campi differenti dalla Medicina. Ad esempio per prevedere che un generico cliente di un'azienda, che ha sottoscritto a t = 0 dei servizi, "sopravviva come cliente" e continui ad usare e pagare tali servizi.

Sempre nel mio repository (corso2) ho collocato un Notebook in cui analizzo il dataset di Kaggle "Telco Customer Churn" e ne derivo la Survival Function. Il grafico lo riporto per comodità del lettore: