- 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.
