giovedì 3 settembre 2009

Analisi della sopravvivenza: curve di Kaplan-Meier e modello di Cox

Le funzioni di sopravvivenza possono essere studiate con due metodi: il metodo attuariale e il metodo di Kaplan-Meier. Quest'ultimo viene utilizzato quando:
- si conosce il tempo esatto in cui avviene l'evento considerato (e non l'intervallo, come nel metodo atturiale)
- i persi al follow-up (truncated o censored, in inglese) partecipano alla determinazione della funzione di sopravvivenza fino al loro ritiro, perchè considerati a rischio fino a quel momento
Così la curva di Kaplan-Meier (e quindi la probabilità di non-evento) cambia ogni volta che accade l'evento nel gruppo considerato.

Supponiamo di voler analizzare la curva di sopravvivenza di 10 soggetti, con il metodo di Kaplan-Meier. Supponiamo che questo gruppo sia affetto da una certa malattia, e siamo quindi interessati a stimare la probabilità che un soggetto ha di sopravvivere fino ad un certo tempo.
I dati vengono raccolti annualmente, e sono i seguenti (per assurdo ammettiamo che gli eventi trocato o deceduto avvengano esattamente alla fine di ogni anno; in realtà si dovrebbe considerare il momento esatto, come si fa nell'esempio successivo):


Tempo Numero.Morti Numero.Troncati
1 0 1
3 1 0
4 0 1
5 2 0
6 0 1
7 2 1
8 1 0


Una tabella di questo genere viene letta nel modo seguente: alla fine del primo anno non è morto nessun paziente, ed è stato perso al follow-up 1 paziente; al termine del terzo anno è morto 1 paziente; e così via.

In R il pacchetto per l'analisi delle funzione di sopravvivenza è la library survival. Il formato con cui sono stati presentati i dati, non è però quello ottimale per l'analisi in R. La funzione che tra breve useremo richiede in input questi stessi dati, presentati in questo modo:


time status
1 0
3 1
4 0
5 1
5 1
6 0
7 1
7 1
7 0
8 1


Per ogni paziente in osservazione, viene definito uno status:
0 per i troncati
1 per i deceduti
Questa tabella viene letta nel modo seguente: al termine del primo anno un paziente è stato perso al follow-up; se in un certo tempo, più di un paziente viene perso o è deceduto, si scrive una riga per ognuno di essi. Quindi osservando il settimo anno, si scrivono 3 righe, due relative ai 2 pazienti deceduti, e una relativa al paziente troncato
In pratica ogni volta che accade un certo evento (paziente troncato o deceduto), si scrive una riga per paziente, specificando il tempo in cui è avvenuto e cosa è avvenuto.

Creiamo quindi un dataframe in R con i seguenti dati:


time <- c(1,3,4,5,5,6,7,7,7,8)
status <- c(0,1,0,1,1,0,1,1,0,1)
popol <- data.frame(time, status)


A questo punto possiamo generare la curva di Kaplan-Meier:


library(survival)
fit <- survfit(Surv(time, status) ~ 1, data=popol)
plot(fit)
text( fit$time, 0, format(fit$n.risk), cex = 0.7 )




Osserviamo il grafico. Ciò che interessa a noi è la linea continua, che indica la probabilità di sopravvivenza in funzione del tempo; le linee discontinue fanno riferimento all'intervallo di confidenza della probabilità di sopravvivenza; se non vogliamo l'intervallo di confidenza, scriviamo plot(fit, conf.int=F). I trattini che possiamo vedere, in corrispondenza dei tempi 1,4,6 indicano i soggetti che sono stati censurati (persi). I numeri posizionati col comando text sull'asse X indicano il numero di soggetti in osservazione ai vari tempi.
Con la funzione summary(fit) otteniamo le probabilità in corrispondenza dei punti dove la curva cambia:


fit <- survfit(Surv(time, status) ~ 1, data=popol, se.fit=TRUE, conf.int=.95, type="kaplan-meier")
summary(fit)

Call: survfit(formula = Surv(time, status) ~ 1, data = popol, se.fit = TRUE,
conf.int = 0.95)

time n.risk n.event survival std.err lower 95% CI upper 95% CI
3 9 1 0.889 0.105 0.706 1.000
5 7 2 0.635 0.169 0.377 1.000
7 4 2 0.317 0.180 0.105 0.964
8 1 1 0.000 NaN NA NA


La colonna survival della tabella indica la probabilità di sopravvivenza (l'asse Y del grafico di prima). Ad esempio, deduciamo da questa che la probabilità (cumulativa, e non puntuale) di sopravvivere al termine del 5° anno è del 63.5%. Essendo il gruppo A costituito da 10 soggetti, che al termine dell'ottavo anno sono tutti scomparsi (persi oppure deceduti), risulta ovvio che la probabilità di sopravvivere al termine dell'ottavo anno è pari a zero.




Come ulteriore esercizio delle funzioni viste fin qui, utilizziamo i dati estratti dalla tabella pubblicata a pagina 2 di questo pdf.
Riassiumiamo rapidamente i comandi per generare la Kaplan-Meier survival curve, e per ottenere l'ultima colonna (Cumulative proportion surviving).


time <- c(0.909, 1.112, 1.322, 1.328, 1.536, 2.713, 2.741, 2.743, 3.524, 4.079)
status <- c(1,1,0,1,1,1,0,1,0,0)
table1 <- data.frame(time, status)

fit <- survfit(Surv(time, status) ~ 1, data=table1)

summary(fit)

Call: survfit(formula = Surv(time, status) ~ 1, data = table1)

time n.risk n.event survival std.err lower 95% CI upper 95% CI
0.909 10 1 0.900 0.0949 0.732 1.000
1.112 9 1 0.800 0.1265 0.587 1.000
1.328 7 1 0.686 0.1515 0.445 1.000
1.536 6 1 0.571 0.1638 0.326 1.000
2.713 5 1 0.457 0.1662 0.224 0.932
2.743 3 1 0.305 0.1666 0.104 0.890


plot(fit, col="red", lwd=2, conf.int=F)







Vediamo adesso come confrontare due curve di sopravvivenza di Kaplan-Meier utilizzando il Log-rank test (o test di Mantel-Cox).
Si dispone di due gruppi (A e B), e dei seguenti dati:


Gruppo A
Tempo Numero.Morti Numero.Troncati
1 0 1
3 1 0
4 0 1
5 2 0
6 0 1
7 2 1
8 1 0

Gruppo B
Tempo Numero.Morti Numero.Troncati
2 2 0
3 0 1
4 1 0
6 0 2
7 1 0
10 1 0


Come abbiamo visto prima, questa forma non è quella idonea per le funzioni di R. Leggiamo le due tabelle in questo modo:


Gruppo A
1 soggetto troncato al tempo 1
1 soggetto morto al tempo 3
1 soggetto troncato al tempo 4
2 soggetti morti al tempo 5
1 soggetto troncato al tempo 6
2 soggetti morti e 1 soggetto troncato al tempo 7
1 soggetto morto al tempo 8

Gruppo B
2 soggetti morti al tempo 2
1 soggetto troncato al tempo 3
1 soggetto morto al tempo 4
2 soggetti troncati al tempo 6
1 soggetto morto al tempo 7
1 soggetto morto al tempo 10


Traduciamo quindi le tabelle nel seguente formato:


Gruppo A
time status
1 0
3 1
4 0
5 1
5 1
6 0
7 1
7 1
7 0
8 1

Gruppo B
time status
2 1
2 1
3 0
4 1
6 0
6 0
7 1
10 1


Adesso creiamo un dataframe contenente questi valori, più la variabile che specifica il gruppo di appartenenza (che trasformiamo in un fattore):


time <- c(1,3,4,5,5,6,7,7,7,8,2,2,3,4,6,6,7,10)
status <- c(0,1,0,1,1,0,1,1,0,1,1,1,0,1,0,0,1,1)
groups <- c(rep("GroupA",10), rep("GroupB",8))
pop <- data.frame(time, status, groups)
pop$groups <- factor(pop$groups)


Adesso disegnano su un unico grafico le due curve di Kaplan-Meier. Dapprima fittiamo il modello: anzichè specificare ~1, il modello deve essere in funzione del gruppo di appartenenza, e quindi abbiamo ~groups:


fit <- survfit(Surv(time, status) ~ groups, data=pop)
summary(fit)

Call: survfit(formula = Surv(time, status) ~ groups, data = pop)

groups=GroupA
time n.risk n.event survival std.err lower 95% CI upper 95% CI
3 9 1 0.889 0.105 0.706 1.000
5 7 2 0.635 0.169 0.377 1.000
7 4 2 0.317 0.180 0.105 0.964
8 1 1 0.000 NaN NA NA

groups=GroupB
time n.risk n.event survival std.err lower 95% CI upper 95% CI
2 8 2 0.75 0.153 0.5027 1
4 5 1 0.60 0.182 0.3315 1
7 2 1 0.30 0.231 0.0664 1
10 1 1 0.00 NaN NA NA

plot(fit, col=2:3, lwd=2)




Ora confrontiamo le due curve con il Log-rank test:


survdiff(Surv(time, status) ~ groups, data=pop)

Call:
survdiff(formula = Surv(time, status) ~ groups, data = pop)

N Observed Expected (O-E)^2/E (O-E)^2/V
groups=GroupA 10 6 6.05 0.000364 0.00110
groups=GroupB 8 5 4.95 0.000445 0.00110

Chisq= 0 on 1 degrees of freedom, p= 0.974


Commentiamo l'output: la colonna N indica il numero di soggetti analizzati per ciascun gruppo; la colonna observed è il numero di morti osservate (O); la colonna expected è il numero di morti attese per ciascun gruppo (E). Infine vengono riportati i valori per il calcolo del chi-quadro della statistica test.
Il chi-quadro calcolato è pari alla somma dei valori della colonna (O-E)^2/E; esso ci viene dato anche alla fine, e ha valore 0 (approssimato, in realtà è 0.000809). Il relativo p-value (0.974) è abbondantemente superiore a 0.05, quindi accettiamo l'ipotesi nulla: l'andamento della sopravvivenza nei due gruppi è significativamente uguale.

Per calcolare il relative risk del gruppoA rispetto al gruppoB, si applica questa formula:



Calcoliamo in R utilizzando i valori della tabella ottenuta:


RR <- (6/6.05) / (5/4.95)
RR
[1] 0.9818182


Avendo ottenuto un Log-rank test non significativo, è abbastanza prevedibile ottenere anche un relative-risk molto vicino ad 1 (il rischio di morte del gruppo A è pressochè uguale a quello del gruppo B).




Vediamo ora un altro esempio di confronto di due curve KM. Utilizzo qui i dati presenti in questo link. Vediamo subito che uno dei modi di presentare i dati è il seguente:


Group 1: 143, 165, 188, 188, 190, 192, 206, 208, 212, 216, 220, 227, 230, 235, 246, 265, 303, 216*, 244*
Group 2: 142, 157, 163, 198, 205, 232, 232, 232, 233, 233, 233, 233, 239, 240, 261, 280, 280, 295, 295, 323, 204*, 344*


Per ogni individuo in analisi, viene fornito il tempo di sopravvivenza prima dell'evento, mentre vengono asteriscati i soggetti persi al follow-up (censored). Traduciamo questi dati in R:


Time1 <- c(143, 165, 188, 188, 190, 192, 206, 208, 212, 216, 220, 227, 230, 235, 246, 265, 303, 216, 244)
Time2 <- c(142, 157, 163, 198, 205, 232, 232, 232, 233, 233, 233, 233, 239, 240, 261, 280, 280, 295, 295, 323, 204, 344)

Status1 <- c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0)
Status2 <- c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0)

Gru1 <- rep("Group1", length(Status1))
Gru2 <- rep("Group2", length(Status2))

Time <- c(Time1,Time2)
Status <- c(Status1,Status2)
Group <- factor(c(Gru1,Gru2))


Osserviamo il formato della funzione Surv:


> Surv(Time1, Status1)
[1] 143 165 188 188 190 192 206 208 212 216 220 227 230 235 246 265
[17] 303 216+ 244+

> Surv(Time2, Status2)
[1] 142 157 163 198 205 232 232 232 233 233 233 233 239 240 261 280
[17] 280 295 295 323 204+ 344+


Come possiamo vedere, la funzione Surv trasforma i dati nel formato descritto prima. Adesso fittiamo il modello:


KM.fit <- survfit(Surv(Time,Status) ~ Group)
summary(KM.fit)

Call: survfit(formula = Surv(Time, Status) ~ Group)

Group=Group1
time n.risk n.event survival std.err lower 95% CI upper 95% CI
143 19 1 0.947 0.0512 0.8521 1.000
165 18 1 0.895 0.0704 0.7669 1.000
188 17 2 0.789 0.0935 0.6259 0.996
190 15 1 0.737 0.1010 0.5632 0.964
192 14 1 0.684 0.1066 0.5041 0.929
206 13 1 0.632 0.1107 0.4480 0.890
208 12 1 0.579 0.1133 0.3946 0.850
212 11 1 0.526 0.1145 0.3435 0.806
216 10 1 0.474 0.1145 0.2949 0.761
220 8 1 0.414 0.1145 0.2412 0.712
227 7 1 0.355 0.1124 0.1911 0.661
230 6 1 0.296 0.1082 0.1447 0.606
235 5 1 0.237 0.1015 0.1023 0.548
246 3 1 0.158 0.0934 0.0495 0.504
265 2 1 0.079 0.0728 0.0130 0.481
303 1 1 0.000 NaN NA NA

Group=Group2
time n.risk n.event survival std.err lower 95% CI upper 95% CI
142 22 1 0.9545 0.0444 0.87136 1.000
157 21 1 0.9091 0.0613 0.79656 1.000
163 20 1 0.8636 0.0732 0.73151 1.000
198 19 1 0.8182 0.0822 0.67189 0.996
205 17 1 0.7701 0.0904 0.61180 0.969
232 16 3 0.6257 0.1051 0.45020 0.870
233 13 4 0.4332 0.1082 0.26548 0.707
239 9 1 0.3850 0.1063 0.22408 0.662
240 8 1 0.3369 0.1034 0.18465 0.615
261 7 1 0.2888 0.0992 0.14731 0.566
280 6 2 0.1925 0.0864 0.07991 0.464
295 4 2 0.0963 0.0647 0.02580 0.359
323 2 1 0.0481 0.0469 0.00712 0.326


Se siamo interessati a conoscere la mediana del tempo di sopravvivenza per ciascun gruppo, e il suo intervallo di confidenza, scriviamo:


KM.fit

Call: survfit(formula = Surv(Time, Status) ~ Group)

records n.max n.start events median 0.95LCL 0.95UCL
Group=Group1 19 19 19 17 216 206 265
Group=Group2 22 22 22 20 233 232 280


Per il gruppo1 la mediana è 216, con intervallo di confidenza (206,265), mentre per il gruppo2 è 233, con CI = (232,280). L'algoritmo per calcolare l'intervallo di confidenza della mediana può essere scelto con il comando conf.type="" (nell'help della funzione survfit sono descritti i vari algoritmi che possiamo scegliere).
In R non ho trovato funzioni per calcolare la media; ma del resto la media nell'analisi delle curve di sopravvivenza di Kaplan-Meier non viene accettata da buona parte degli statistici.

Adesso poniamo in grafico le curve di Kaplan-Meier:


plot(KM.fit, conf.int=F, col=2:3, lwd=2, xlab="time", ylab="Survival probability")
mtext("K-M survival curve for Group1 and Group2", 3, line=2, cex=1.5)
abline(h=0)

text(300, 1, expression(bold("Group 1 - median = 216")), col=2, cex=.8, lwd=2)
text(300, 0.95, expression(bold("Group 2 - median = 233")), col=3, cex=.8, lwd=2)




Adesso eseguiamo il Log-rank test sui due gruppi:


survdiff(Surv(Time, Status) ~ Group)

Call:
survdiff(formula = Surv(Time, Status) ~ Group)

N Observed Expected (O-E)^2/E (O-E)^2/V
Group=Group1 19 17 12.2 1.885 3.12
Group=Group2 22 20 24.8 0.928 3.12

Chisq= 3.1 on 1 degrees of freedom, p= 0.0771


Abbiamo ottenuto p-value > 0.05, quindi l'andamento della sopravvivenza nei due gruppi è significativamente (95%) uguale. Il rischio relativo del gruppo1 rispetto al gruppo2 è:


RR <- (17/12.2) / (20/24.8)
RR
[1] 1.727869


Il gruppo1 ha un rischio di morte 1.73 volte più elevato del gruppo2.
Se vogliamo esaminare la differenza delle probabilità di sopravvivenza calcolata per ogni tempo esaminato, e il suo intervallo in confidenza, utilizziamo la funzione survdiffplot disponibile nel pacchetto Design:


library(Design)
survdiffplot(KM.fit)




Per ottenere il grafico dell'hazard in funzione del tempo, dobbiamo utilizzare le funzioni coxph() e quindi basehaz() e plottare il risultato:


temp <- coxph(Surv(Time, Status)~strata(Group))
plot(basehaz(temp))





Proportional hazard model: analisi della sopravvivenza col modello di Cox

Il Log-rank test che abbiamo visto precedentemente permette di confrontare l'andamento della sopravvivenza in più gruppi, che potrebbero anche essere suddivisi in base ad un certo predittore (ad esempio si potrebbero confrontare le curve di sopravvivenza nei maschi e nelle femmine). Esso però non ci fornisce alcuna idea sul peso che assume la regola con la quale dividiamo i gruppi (ossia non ci dice, ad esempio, quanto influisce la variabile "sesso" sulla sopravvivenza).
Il modello di Cox è il modello di regressione sviluppato proprio per poter analizzare il peso di certe variabili sulla sopravvivenza. Esso altro non è che una regressione multipla, di cui segue gli stessi principi.
In presenza di un certo numero di fattori di rischio (probabili), chiamate covariate, è possibile stimare il rischio di morte (hazard) ad un certo tempo t, in presenza delle covariate z. La funzione di rischio (che rappresenta la probabilità di morire al tempo t per un individuo con vettore di covariate z, sopravvissuto fino al tempo t) viene generalmente fattorizzata in due parti:


Qui è la funzione di rischio di base (baseline hazard) per un individuo con vettore di covariate z=0 che dipende solo da t ed è lasciata non specificata; è la funzione hazard, con beta i regressori (il peso di ciascuna covariate z).
Nel modello più semplice, la funzione hazard H è quella esponenziale:


Il modello di Cox è definito semi-parametrico. Ma cominciamo subito con un esempio.

Carichiamo il pacchetto survival, e utilizziamo il dataset ovarian, ivi disponibile:


library(survival)
Carico il pacchetto richiesto: splines

ovarian
futime fustat age resid.ds rx ecog.ps
1 59 1 72.3315 2 1 1
2 115 1 74.4932 2 1 1
3 156 1 66.4658 2 1 2
4 421 0 53.3644 2 2 1
5 431 1 50.3397 2 1 1
6 448 0 56.4301 1 1 2
7 464 1 56.9370 2 2 2
8 475 1 59.8548 2 2 2
9 477 0 64.1753 2 1 1
10 563 1 55.1781 1 2 2
11 638 1 56.7562 1 1 2
12 744 0 50.1096 1 2 1
13 769 0 59.6301 2 2 2
14 770 0 57.0521 2 2 1
15 803 0 39.2712 1 1 1
16 855 0 43.1233 1 1 2
17 1040 0 38.8932 2 1 2
18 1106 0 44.6000 1 1 1
19 1129 0 53.9068 1 2 1
20 1206 0 44.2055 2 2 1
21 1227 0 59.5890 1 2 2
22 268 1 74.5041 2 1 2
23 329 1 43.1370 2 1 1
24 353 1 63.2192 1 2 2
25 365 1 64.4247 2 2 1
26 377 0 58.3096 1 2 1


Le colonne futime e fustat sono rispettivamente il tempo di analisi dei 26 soggetti, e lo status (morto e censored); rx rappresenta i due gruppi. Vogliamo stimare il peso dei regressori age (età) nei due gruppi rx. Utilizziamo la funzione coxph():


cph1 <- coxph(Surv(futime,fustat) ~ rx+age, data=ovarian)
summary(cph1)

Call:
coxph(formula = Surv(futime, fustat) ~ rx + age, data = ovarian)

n= 26

coef exp(coef) se(coef) z Pr(>|z|)
rx -0.80397 0.44755 0.63205 -1.272 0.20337
age 0.14733 1.15873 0.04615 3.193 0.00141 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

exp(coef) exp(-coef) lower .95 upper .95
rx 0.4475 2.234 0.1297 1.545
age 1.1587 0.863 1.0585 1.268

Rsquare= 0.457 (max possible= 0.932 )
Likelihood ratio test= 15.89 on 2 df, p=0.0003551
Wald test = 13.47 on 2 df, p=0.00119
Score (logrank) test = 18.56 on 2 df, p=9.34e-05


Scriviamo subito l'equazione del nostro modello:


Qual è il rischio relativo di un soggetto appartenente al gruppo rx=1 di 50 anni, rispetto ad un soggetto appartenente al gruppo rx=2 di 60 anni?


Possiamo mettere in grafico il modello, con la funzione plot(survfit(cph1)). Vengono consideato tutte le covariate prese in considerazione dal modello.

Nell'output vengono anche forniti i valori z e p-value per verificare se i coefficienti beta delle covariate sono significativi. La colonna zeta è il test di Wald per ogni regressore, ottenuta dividendo il valore beta per lo standard error del coefficiente stesso.
In questo caso è significativo solamente il regressore relativo alla covariata age.
Il likelihood ratio test, il Wald test e la statistica score chi-square che osserviamo alla fine rappresentano i test statistici per verificare se tutti i regressori considerati insieme sono differenti da zero (serve quindi per validare il modello fittato). In questo caso p-value < 0.05 in tutti i test, e quindi l'ipotesi nulla viene rifiutata.




Analizziamo un altro esempio, utilizzando il dataset qui disponibile. Poichè nel txt non è specificato il nome delle colonne, lo assegnamo noi con la funzione names():


data <- read.table("http://www.sph.emory.edu/~cdckms/CoxPH/Anderson_i.txt", F, sep=",")
names(data) <- c("pred1", "pred2", "pred3", "time", "status")
data[1:10, 1:5]

pred1 pred2 pred3 time status
1 0 1.45 0 35 0
2 0 1.47 0 34 0
3 0 2.20 0 32 0
4 0 2.53 0 32 0
5 0 1.78 0 25 0
6 0 2.57 0 23 1
7 0 2.32 0 22 1
8 0 2.01 0 20 0
9 0 2.05 0 19 0
10 0 2.16 0 17 0


Fittiamo il modello di regressione di Cox, conoscendo che esiste un'interazione tra pred2 e pred3:


coxph1 <- coxph(Surv(time,status)~pred1+pred2*pred3, data=data)
summary(coxph1)

Call:
coxph(formula = Surv(time, status) ~ pred1 + pred2 * pred3, data = data)

n= 42

coef exp(coef) se(coef) z Pr(>|z|)
pred1 8.01642 3030.32167 3.06634 2.614 0.00894 **
pred2 1.84990 6.35920 0.44718 4.137 3.52e-05 ***
pred3 -4.12837 0.01611 1.88621 -2.189 0.02862 *
pred2:pred3 0.60680 1.83455 0.29746 2.040 0.04136 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

exp(coef) exp(-coef) lower .95 upper .95
pred1 3.030e+03 0.00033 7.4372579 1.235e+06
pred2 6.359e+00 0.15725 2.6470659 1.528e+01
pred3 1.611e-02 62.07646 0.0003995 6.496e-01
pred2:pred3 1.835e+00 0.54509 1.0240708 3.286e+00

Rsquare= 0.704 (max possible= 0.988 )
Likelihood ratio test= 51.15 on 4 df, p=2.072e-10
Wald test = 35.56 on 4 df, p=3.57e-07
Score (logrank) test = 71.79 on 4 df, p=9.548e-15


Il p-value di tutti i regressori è minore di 0.05, così come quelli relativi ai test sul modello completo.

Una delle assunzioni del modello è il "proportional hazards" (i coefficienti del modello non dipendono dal tempo), il quale può essere esaminato con la funzione cox.zph(), per effettuare il test di rischio proporzionale:


cox.zph(coxph1)

rho chisq p
pred1 0.1303 0.587 0.444
pred2 0.2117 1.089 0.297
pred3 -0.1165 0.544 0.461
pred2:pred3 0.0936 0.394 0.530
GLOBAL NA 1.716 0.788


Essendo il p-value > 0.05 per tutte le covariate e per il modello globale, l'assunzione è stata verificata.
Possiamo anche verificare graficamente i residui (di Schoenfeld) plottando la funzione:


par(mfrow=c(2,2))
plot(cox.zph(coxph1))




I punti rispettano una certa linearità (come richiesto dall'assunzione); se i punti fossero stati maggiormente dispersi, non avremmo potuto verificare l'assuzione di proporzionalità, e quindi si dovrebbe procedere con trasformazioni della variabile.

Un'altra verifica che andrebbe affettuata è quella della non-linearità delle covariate (salvo quelle dicotomiche). In questo pdf nella sezione 5.3 sono presenti i test da effettuare.

13 commenti:

  1. Spiegazione molto esaustiva, però forse ho perso un passaggio: tramite quale comando posso inserire il numero iniziale di pazienti per realizzare una curva di Kaplan Meier? Perché, omettendo questo valore, ho ottenuto una curva diversa da quella che mi ero calcolata usando Excell.
    Grazie e complimenti per il lavoro!!

    RispondiElimina
  2. Rispondo subito alla domanda: in R in realtà non c'è bisogno di specificare il numero di soggetti presi in esame, perchè questo valore viene automaticamente computato dalla funzioni viste, contando il numero di righe del dataframe che creiamo. Se ad esempio nel nostro lavoro abbiamo un gruppo iniziale di 10 soggetti, dovremo creare dieci righe (colonne time e status) che descrivono ciascun soggetto. Questo perchè, tecnicamente, lavori sulla sopravvivenza dovrebbero concludersi solo nel momento in cui tutti i soggetti analizzati vengono persi al follow up, oppure vanno incontro all'evento da noi ricercato (se ad esempio decidiamo di considerare 100 soggetti, non potremo interrompere il lavoro, ad esempio, con ancora 20 soggetti non-persi e senza-evento: in questo caso la curva andrebbe calcolata su 80 soggetti).
    Spero così di aver risolto il tuo dubbio. A disposizione :)

    RispondiElimina
  3. Grazie mille per la spiegazione! :-) Avevo omesso di censurare i pazienti che sono arrivati fino alla fine del follow up, quindi -ovviamente- i conti non tornavano.
    Vorrei approfittare della tua disponibilità per chiederti un'altra cosa, anche se esula dalla spiegazione. Se volessi applicare il modello di Cox per studiare una coorte di circa 90 pazienti, quante variabili potrei considerare al massimo? Grazie ancora!

    RispondiElimina
  4. Spiegazione molto chiara, ma poco esaustiva per quanto riguarda il calcolo del log rank test. La differenza tra le morti attese e le morti osservate e' la stessa per i due gruppi. Questo porta ad utilizzare una formula che temo non sia utilizzabile nel caso, molto piu' comune, in cui la suddetta differenza non dia lo stesso valore numerico per i due gruppi. Qualcuno sa dirmi che formula bisogna utilizzare in questo caso?
    Grazie e complimenti per il sito!
    Fabio

    RispondiElimina
  5. io devo la vita a questa spiegazione! un'unica domanda:io ho utilizzato la curva di K-M per valutare l'efficia di un farmaco quindi mi aspetto una curva che non sia descrescente ma CRESCENTE....come faccio?? esistono dei comandi?...
    se mi rispondesse sappia che le invierò un regalo per la mia laurea (che è l'11 ottobre)!

    RispondiElimina
  6. @brasella
    sono molto contento di esserti stato d'aiuto!
    Sinceramente non so dirti come capovolgere il grafico, in quanto non m'è mai capitato un problema del genere. Potresti provare ad usare la funzione reverse() e vedere cosa succede, se l'input dei dati è effettivamente invertito. Oppure la stessa funzione reverse, la puoi applicare nel grafico soltanto (ma in questo caso dovresti modificare il codice della funzione).
    Ma se il tuo problema è solamente "grafico", mentre i valori che ottieni sono tutti giusti, beh, ti consiglio di usare mspaint, o inkscape, e invertire i tuoi dati. Soluzione spartana, ma è inutile "sbattere la testa" a trovare una soluzione troppo complessa =)
    Diversamente non sarei dirti, in quanto non ho mai affrontato il problema personalmente.

    In bocca al lupo per la tua laurea!

    RispondiElimina
  7. Ciao! ottimo post!
    solo una cosa non mi e' chiara, una volta stimato il modello di Cox. Io vorrei poter dire qual'e la probabilita' per un soggetto x con caratteristiche y,z di avere l'evento.

    Grazie!

    RispondiElimina
  8. bellissime spiegazioni per l'uso di R... cose che neanche credevo possibili.
    Io ho una domanda di statistica, giacchè mi sembri esperto.
    Se analizzo l'effetto di un farmaco sulla comparsa di malattia (recidiva) con le curve di sopravvivenza trovo un log rank significativo, se lo analizzo con un chi quadro (malati si malati no - tp si tp no) non mi da significatività (p>0,05). Perchè c'è discordanza???

    RispondiElimina
  9. Ciao ottimo articolo mi è stato molto utile! volevo chiedere però se c'è un metodo per plottare le curve in modo che nell'asse x sia utilizzato un intervallo di tempo di 6 mesi (di default segna ogni 5 mesi). grazie!!

    RispondiElimina
  10. Ciao, scusate non riesco a fare il bootstrap del modello multivariato di Cox con R, chi mi può aiutare?

    RispondiElimina
  11. Grazie per il chiarissimo e coinciso articolo.

    Vorrei suggerire la seguente funzione
    summary(fit,time=c(24,60))

    per il calcolo della sopravvivenza a un determinato valore temporale. Nell'esempio a 24 e 60 mesi.

    Confermi la validità di questa funzione?

    Grazie ancora

    RispondiElimina
  12. ciao! ma se all'analisi di kaplan tra due gruppi elativamente ad un evento x non c è differenza, ovvero log rank non significativo
    e facendo la cox invece una covariata risulta statisticamente significativa con hr intervallo ecc. COME SI INTERPRETA? diventa un fattore di rischio predittivo dell evento in x indipendentemente dall 'appartenza ad un gruppo ????

    RispondiElimina