Cox regression

Cox regression är synonymt med överlevnadsanalys. Detta beror på metodens många fördelar gentemot andra regressionsmodeller. Alternativen till Cox regression, exempelvis Weibull regression, ger under vissa omständigheter mer korrekta modeller, på bekostnad av ökad komplexitet. Cox regression är en enkel och kraftfull modell för att studera överlevnad. Skillnaden mellan Cox regression och Kaplan-Meier estimator är att förstnämnda är en regressionsmodell, vilket innebär att vi kan studera hur flera prediktorer påverkar överlevnaden. Cox regressioner kan, likt vanlig multipel regression, inkludera en eller flera prediktorer.

Kapitlen Överlevnadsanalys och Kaplan-Meier Estimator bör läsas innan detta kapitel. Förståelse för alla nedanstående formler är inte nödvändigt för att använda Cox regression. Efter den matematiska diskussionen nedan följer förklaringar och praktiska exempel på utförande och tolkning av Cox regressioner.

Med Cox regression studeras sambandet mellan prediktorer (kovariater) och survival time (observationstid). För förstå Cox regression tar vi avstamp i hazard function, \(h(t)\), och survival function, \(S(t)\).

• \(S(t)\) – Överlevnadsfunktionen (survival function) är sannolikheten för att en person skall vara vid liv efter tidpunkt \(t\).
  • Tolkning: Vad är sannolikheten för att jag lever om 1 år?
• \(\lambda(t)\) – Hazardfunktionen (hazard function, hazard rate) är den momentana risken (sannolikheten) för att händelsen skall inträffa vid tidpunkt \(t\), förutsatt att personen överlevt till \(t\).
  • Tolkning: Vad är sannolikheten för att jag dör om 1 år, givet att jag överlevt fram till nu?

\(S(t)\) och \(\lambda(t)\) har egenskaper som gör de lämpliga för att studera överlevnad. Inte minst framgår detta av de relevanta frågor som funktionerna besvarar (se tolkning ovan).

Survival time (överlevnadstid, T) är en kontinuerlig variabel (ofta specificerad som antal dagar) och kan därför beskrivas med en kumulativ fördelningsfunktion (cumulative distribution function, CDF):

\(F(x) = P(X≤x)\)

\(F(x)\) beskriver sannolikheten för att \(X\) är mindre eller lika med \(x\). Överlevnadstid betecknas med \(T\) och kan också beskrivas med en kumulativ fördelningsfunktion:

\(F(t) = P(T≤t)\)

\(F(t)\) beskriver sannolikheten för att händelsen inträffar senast på tidpunkt \(t\). Sannolikheten för att ett nyfött barn skall dö vid 90 års ålder är cirka 1%, men om barnet tar sig till 89 års ålder är sannolikheten för död på 90 års dagen 10%. Med \(F(t)\) kan vi skatta \(S(t)\), som är sannolikheten för att individen skall överleva tidpunkt \(t\):

\(S(t) = 1 - F(t) = P(T > t)\)

\(P(T > t)\) är sannolikheten för att överleva tidpunkt \(t\).

Hazardfunktionen (hazard function) har följande formel:

\(\lambda(t) = \lim_{\Delta t \to 0} \frac{P \left( t \leq T < t + \Delta t \vert T \geq t \right)}{\Delta t}\)

\(\lim_{\Delta t \to 0}\) innebär att kvoten gäller då risk (hazard) studeras under en mycket kort tidsperiod, dvs då \(\Delta t\) (skillnaden i \(t\)) går mot \(0\). Täljaren i \(\lambda(t)\) innehåller | vilket betyder att det är en betingad sannolikhet. Täljaren beskriver sannolikheten för att händelsen skall inträffa under ett litet tidsintervall mellan \(t\) och \(t+\Delta t\), betingat på att individen överlevde fram till \(t\). Nämnaren indikerar tidsintervallets längd (\(\Delta t\)). Hazardfunktionen (\(\lambda(t)\)) är således ekvivalent med event rate (hastigheten med vilken händelser inträffar) för individer som ännu inte drabbats av händelsen.

Täljaren i \(\lambda(t)\) utgör kvoten mellan sannolikheten för att T ska vara i intervallet \(t\) till \(t+\Delta t\) och att \(P(T \geq t)\), varav sistnämnda är ekvivalent med \(S(t)\) och förstnämnda kan omformuleras som \(f(t)\Delta t\) för små \(\Delta t\). Genom att dividera den kvoten med \(\Delta t\) erhåller vi följande formel:

\(\lambda(t) = \frac{f(t)}{S(t)}\)

Detta är den simpla hazardfunktionen. Den förtäljer att hastigheten med vilken händelser inträffar vid tidpunkt \(t\) är ekvivalent med tätheten (eng. density) för events vid tidpunkt \(t\), dividerat med sannolikheten för att överleva fram till tidpunkt \(t\) utan att drabbas av händelsen. Täljaren \(f(t)\) är en täthetsfunktion (probability density function, PDF). Eftersom \(–f(t)\) är derivatan till \(S(t)\), kan formeln skrivas som följer:

\(\lambda(t) = - \frac{d}{dt}\log S(t)\)

Om denna integreras från \(0\) till \(t\), och \(S(0)=1\), kan följande formel användas för att beskriva sannolikheten för att överleva till tidpunkt \(t\) som en funktion av hazard vid alla tidpunkter fram till \(t\).

\(S(t) = e^{ - \int_0^t \lambda(x) \Delta x }\)

Integralfunktionen ovan är cumulative hazard (cumulative risk) och formuleras som följer:

\(\Lambda(t) = \int_0^t \lambda(x) \Delta x = -log[1-F(t)] = -log[S(t)]\)

Kumulativ hazard är den ackumulerade risken en individ förvärvar när den passerar från \(0\) till \(t\).

Med \(S(t)\) kan vi alltid beräkna \(f(t)\) och därmed erhålla \(\lambda(t)\). Med \(\lambda(t)\) kan vi alltid beräkna kumulativ hazard och exponentiera den för att erhålla \(S(t)\).

En individ som överlevt fram till tidpunkt \(T\) är som regel mer intresserad av sin betingade sannolikhet för överlevnad, snarare än den obetingade. I praxis betyder detta att en individ som överlevt fram till tidpunkt \(T\) efter en cancerdiagnos är mer intresserad av sin förväntade överlevnad från tidpunkt \(T\). Dessutom är hazardfunktionen en enkel och pragmatisk beräkning när observationerna är högercensurerade.

Cox regression

Med Cox regression beskrivs hazard för en individ med prediktorerna (kovariater) \(X_1, X_2, ..., X_K\) som följer:

\(\lambda(t) = \lambda_0(t) e^{\beta_1 X_1 + \beta_2 X_2 + \ldots + \beta_K X_K}\)

Där \(\lambda_0(t)\) är baseline hazard, som är individens risk när alla variabler har värdet 0 (för kontinuerliga variabler) eller är på referensnivån (för kategoriska variabler), och \(e^{\beta_1 X_1 + \beta_2 X_2 + \ldots + \beta_K X_K}\) är exponentierade prediktorer. Den enda kvantiteten ovan som inbegriper \(t\) är baseline hazard. Som framgår nedan kommer baseline hazard automatiskt att elimineras och därmed blir hazard oberoende av \(t\).

Om alla prediktorer har värdet 0 elimineras hela exponenten och formeln blir som följer:

\(\lambda(t) = \lambda_0(t)\).

Ovanstående formel är definitionen av baseline hazard, vilket alltså är hazard när alla prediktorer har värdet 0. Baseline hazard är ointressant i Cox regression. Det görs inga antaganden om baseline hazard och den får variera godtyckligt. Förklaringen till detta är att syftet med Cox regression är att jämföra exponeringar (exempelvis behandlingar, interventioner, etc) och när två individer jämförs så eliminerar de varandras baseline hazard, vilket demonstreras nedan.

Om båda sidor av ekvationen logaritmeras erhåller vi en modell där prediktorernas effekt är additiva (dvs de har en oberoende inverkan på utfallet som studeras) och linjära, vilket påminner om klassisk linjär regression.

\(\log \lambda(t) = \log [ \lambda_0(t) e^{\beta_1 X_1 + \beta_2 X_2 + \ldots + \beta_K X_K}]\)

\(=\)

\(\log \lambda(t) = \log [ \lambda_0(t)] + {\beta_1 X_1 + \beta_2 X_2 + \ldots + \beta_K X_K}\)

Ponera att vi studerar associationen mellan prediktorn kön och överlevnad. Kön är definierad som följer:

\(
X = \begin{cases} 1 = {Man, referensgrupp} \\ 2 = {Kvinna}
\end{cases}
\)

För denna variabel är man referensgrupp. Cox regressionen för denna modell är som följer:

\(\lambda(t) = \lambda_0(t) e^{\beta_1 X_1}\)

Hazard ratio (HR, hazardkvot)

Hazard ratio (HR) beräknas som kvoten mellan gruppernas hazard, där referensgruppen (män) utgör nämnaren:

\(HR = \frac{\lambda_{Kvinna}(t)}{\lambda_{Man}(t)} = \frac{\lambda_0(t)e^{\beta_1 X_1}}{\lambda_0(t)}=e^{\beta_1 X_1}\)

Effekten av \(\beta\) är samma vid alla \(t\). Formeln innebär och fungerar endast om hazard ratio inte varierar över tid. Cox modell kallas även Cox proportional hazards model av denna anledningen.

\( \lambda(t) = \lambda_0(t) e^{\beta_1 X_1 + \ldots + \beta_K X_K}\)

Cox regression innebär alltså att hazard (risk) vid tidpunkt \(t\) är produkten av baseline hazard och exponenten av prediktorernas sammanlagda effekt.

Cox regression i praxis

Baseline hazard i Cox regression

\(\lambda_0(t)\) är baseline hazard och denna kvantitet kan vara svår att förstå. Baseline hazard är den risk som individen har när värdet på alla prediktorer är 0. Två exempel följer:

• Ponera att vi studerar hur ålder påverkar dödlighet. I en Cox regression där ålder är den enda prediktorn skulle baseline hazard vara den risk som människan har vid ålder = 0 år. Efter födseln tillkommer risk när åldern ökar.
• Ponera att vi studerar hur blodtryck påverkar dödlighet. I en sådan modell skulle baseline hazard vara individens risk vid blodtryck = 0, vilket inte är biologiskt meningsfullt (ingen levande person har blodtryck = 0). Detta innebär att baseline hazard är individens risk när alla prediktorer ställs till 0, vilket är detsamma som att modellen saknar prediktorer.

Prediktorer (kovariater) i Cox regression

\(\beta_1 X_1 + \beta_2 X_2 + \ldots + \beta_K X_K\) är de prediktorer (kovariater) som används för att modellera en individs hazard function. Samlingen av prediktorer liknar således en vanlig multipel (linjär) regression. Prediktorerna kan vara kontinuerliga eller kategoriska variabler. Prediktorerna är utövar en additiv effekt, vilket innebär att vare enskild prediktor har en oberoende (av de andra prediktorerna) effekt på utfallet (\(Y\)). Vidare är sambandet mellan prediktorerna och utfallet linjärt. Antagandena om additivitet och linearitet kan modifieras (se nedan).

Proportional hazards

Som nämnt ovan är hazard (risk) vid tidpunkt \(t\) produkten av baseline hazard och exponenten av prediktorernas sammanlagda effekt. Detta stämmer under antagandet att hazard är proportionerlig över tid. Proportional hazard innebär att ratio (kvoten) mellan två individers hazard function skall vara konstant över tid. Detta gäller för alla karaktäristika som definierar individerna. Ett exempel följer.

BLOCKERAD KOD, BILD ELLER TEXT


Du måste logga in för att komma åt allt material.

Kvoten mellan två hazard functions är hazard ratio (HR). Det är alltså den relativa skillnaden i risk för två individer. Eftersom \(\lambda_0(t)\) ingår i både täljaren och nämnaren så elimineras den. I ekvationen finns inte heller tid. Således är det endast värdena på prediktorerna som särskiljer individ 1 och individ 2; varken baseline hazard eller tid kan påverka hazard ratio, som därför är konstant mellan individ 1 och individ 2. Det innebär också att om någon av prediktorerna (\(\beta_1 X_1 + \beta_2 X_2 + \ldots + \beta_K X_K\)) skulle ha en varierande effekt över tid, så kommer kravet om proportional hazards inte vara tillfredsställt och modellen blir ogiltig.

Hazard ratio (hazardkvot)

Hazard ratio beräknas genom att dividera hazard för en individ med hazard för en annan individ (se formel ovan). Hazard ratio beräknas när en händelse inträffar. Om en händelse inträffar vid tidpunkt t så beräknas hazard ratio genom att dividera hazard för individen som erfarde händelsen med hazard för övriga som klarade sig fram till tidpunkt t utan att erfara händelsen.  Baseline hazard utgår (se ovan) från kvoten likaså är tidpunkten irrelevant eftersom alla individer jämförs vid samma tidpunkt. Det som karaktäriserar individerna blir således deras prediktorer (x₁, x₂ … x_k). Således kan effekten av prediktorerna analyseras.

I en Cox regression beräknas således koefficienterna för prediktorerna (\(\beta_1 X_1 + \beta_2 X_2 + \ldots + \beta_K X_K\)), vilket är dessa prediktorers inverkan på överlevnaden. Den hazard ratio som erhålls när två individers hazard divideras är logaritmerad, vilket också gör att den är svårtolkad. För att åtgärda detta så exponentieras koefficienten. Om regressionskoefficienten är 0.5 så blir den exponerade koefficienten e^0.5 = 1.65. När vi använder termen hazard ratio så referar vi egentligen till den exponentierade koefficienten.

Tolkning av hazard ratio

• För en kontinuerlig prediktor så är hazard ratio den ökning/minskning i hazard som följer när prediktorn \(X_1\) ökar med en enhet.
• För en kategorisk variabel så är hazard ratio skillnaden i hazard mellan två (eller flera grupper) där en av grupperna är referensgrupp.

  För att bättre förstå detta behöver vi ett exempel och vi skapar själva detta exempel i R.

Exempel på Cox regression i R

Vi kommer att använda data från en klinisk studie om coloncancer (tjocktarmscancer). Dessa data finns i paketet survival, som också innehåller många viktiga funktioner för överlevnadsanalys.

# Aktivera tidyverse
library(tidyverse)

# Aktivera survival-paketet
library(survival)

# Importera exempeldata
data(cancer)
# Inspektera första raderna
head(colon)

  id study      rx sex age obstruct perfor adhere nodes status differ extent surg node4 time etype
1  1     1 Lev+5FU   1  43        0      0      0     5      1      2      3    0     1 1521     2
2  1     1 Lev+5FU   1  43        0      0      0     5      1      2      3    0     1  968     1
3  2     1 Lev+5FU   1  63        0      0      0     1      0      2      3    0     0 3087     2
4  2     1 Lev+5FU   1  63        0      0      0     1      0      2      3    0     0 3087     1
5  3     1     Obs   0  71        0      0      1     7      1      2      2    0     1  963     2
6  3     1     Obs   0  71        0      0      1     7      1      2      2    0     1  542     1

Vi kodar om vissa variabler så att de korrekt klassificeras som kategoriska:

# Kodar om vissa kategoriska varaibler
colon <- colon |> 
  mutate(sex = recode_factor(sex, '1'='Man', '0'='Kvinna'),
         perfor = recode_factor(perfor, '1'='Yes', '0'='No'),
         obstruct = recode_factor(obstruct, '1'='Yes', '0'='No'))

Nu skapar vi en Cox regression med hjälp av funktionen coxph(). Resultatet sparas i objektet min_modell. Notera att vi måste specificera ett överlevnadsobjekt med hjälp av funktionen Surv(). I Surv() specificerar vi vilka variabler som utgör observationstid (time) och utfallet (status). I detta fallet specificerar vi att status==1 indikerar vilket värde på status som innebär att händelsen inträffat.

BLOCKERAD KOD, BILD ELLER TEXT


Du måste logga in för att komma åt allt material.

Call:
coxph(formula = Surv(time, status == 1) ~ sex + age + nodes + 
    perfor + obstruct, data = colon)

  n= 1822, number of events= 897 
   (36 observations deleted due to missingness)

                coef exp(coef)  se(coef)      z Pr(>|z|)    
sexKvinna   0.057961  1.059674  0.067077  0.864  0.38753    
age         0.001363  1.001364  0.002822  0.483  0.62911    
nodes       0.088342  1.092362  0.006342 13.929  < 2e-16 ***
perforNo   -0.223009  0.800107  0.181919 -1.226  0.22025    
obstructNo -0.249065  0.779529  0.082675 -3.013  0.00259 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

           exp(coef) exp(-coef) lower .95 upper .95
sexKvinna     1.0597     0.9437    0.9291    1.2086
age           1.0014     0.9986    0.9958    1.0069
nodes         1.0924     0.9154    1.0789    1.1060
perforNo      0.8001     1.2498    0.5601    1.1429
obstructNo    0.7795     1.2828    0.6629    0.9167

Concordance= 0.635  (se = 0.009 )
Likelihood ratio test= 145.1  on 5 df,   p=<2e-16
Wald test            = 202.8  on 5 df,   p=<2e-16
Score (logrank) test = 206.2  on 5 df,   p=<2e-16

Ovanstående utskrift är förvisso informativ men den saknar en del information och är inte tillräckligt estetiskt tilltalande för att medbringas i en rapport. Med hjälp av paketet gtsummary kan vi göra en estetiskt tilltalande tabell som är lätt att modifiera. Funktionerna i gtsummary fungerar som resterande funktioner i tidyverse, vilket innebär att man stegvis adderar lager med komponenter, enligt nedanstående. Du kan förmodligen gissa vad de olika funktionerna adderar till tabellen.

BLOCKERAD KOD, BILD ELLER TEXT


Du måste logga in för att komma åt allt material.

Tolkning av Cox regression

Ovanstående modell har fem prediktorer: sex, age, nodes, perfor och obstruct. Sex är en kategorisk prediktor med nivåerna Man och Kvinna. Age och nodes är kontinuerliga prediktorer. Hazard ratio (HR) presenteras tillsammans med ett 95% konfidensintervall (95% CI). Hazard ratios ovan tolkas som följer:

• För män finns inte hazard ratio angivet. Detta beror på att män är referensnivån för variabeln sex. Alla kategoriska variabler har en referensnivå med vilken de övriga nivåerna jämförs.
• Hazard ratio för kvinnor är 1.06, vilket innebär att den relativa risken är 1.06 gånger högre (dvs 6% högre) för kvinnor jämfört med män. Konfidensintervallet är 0.93 till 1.21. Det betyder att med 95% sannolikhet är den sanna risken för kvinnor mellan 0.93 (vilket innebär att risken är 7% lägre för kvinnor jämfört med män) och 1.21 (vilket innebär att risken är 21% högre för kvinnor jämfört med män). Detta konfidensintervall inbegriper alltså siffran 1.0, vilket innebär att det finns en sannolikhet för både lägre och högre risk för kvinnor jämfört med män. Därför är variabeln sex inte statistiskt signifikant. Det framgår även av P-värdet som är 0.4.
• Av ovanstående följer att om ett konfidensintervall inbegriper siffran 1.0 så är prediktorn inte statistiskt signifikant.
• Av ovanstående följer också att en kategorisk prediktor måste ha en referensnivå med vilken övriga nivåer jämförs.
• För variabeln age är hazard ratio 1.00 och konfidensintervallet inbegriper också 1.00, vilket gör att variabeln inte är en statistiskt signifikant prediktor.
• För variabeln nodes är hazard ratio 1.09, vilket innebär att varje gång nodes ökar med 1 enhet, så ökar risken för händelsen (död) med 9%. Konfidensintervallet är 1.08 till 1.11, vilket innebär att det sanna värdet med 95% sannolikhet är mellan 8% och 11% ökad risk för händelsen varje gång nodes stiger med 1 enhet. P-värdet är <0.001, vilket är tydligt signifikant.
• För varje prediktor i modellen är hazard ratio justerat för övriga variabler i modellen. Effekten av att öka nodes med 1 enhet är 1.09 (9% ökad risk), justerat för övriga prediktorer i modellen.
• Denna modellen inkluderar alltså även age och meal.cal, utöver sex. Då är effekten av sex justerad för effekten av de övriga prediktorerna i modellen (de övriga prediktorerna “hålls konstanta”).

Konfidensintervallet i Cox regression

Oftast används Walds konfidensintervall, vilket är en enkel beräkning med standard error (SE):

Nedre intervall: HR – 1.96×SE

Övre intervall: HR + 1.96×SE

Där HR är hazard ratio och SE är standard error.

Mer om proportional hazards (PH)

Proportional hazard kan också förklaras genom interaktioner mellan en prediktor och tid (uppföljningstid). Om en prediktor uppfyller kriteriet för proportional hazard så skall det inte finnas någon statistiskt signifikant interaktion mellan den prediktorn och uppföljningstid (survival time). Man kan således skapa en interaktionsterm (prediktor × uppföljningstid) och utvärdera om den är statistiskt signifikant. Om interaktionen inte är signifikant så innebär det att prediktorn utövar samma effekt på hazard function vid alla värden på t (uppföljningstid) och då är kriteriet för proportional hazard uppfyllt. Om interaktionen däremot är signifikant så innebär det att effekten av prediktorn varierar över tid och då strider den variabeln mot antagandet om proportional hazard. Det finns ytterligare sätt att utvärdera proportional hazards, och dessa är som följer:

• Utvärderig av Schoenfeld residualer
• Visualisera log(-log S(t)) mot t

Egentligen bör man föredra de två sistnämnda metoderna framför interaktionstermer. I R finns enkla funktioner som presenterar både Schoenfeld residualer och grafer över log(-log S(t)) mot t.

För att utvärdera antagandet om PH kan följande kod användas:

BLOCKERAD KOD, BILD ELLER TEXT


Du måste logga in för att komma åt allt material.

Ovanstående grafer visar hur \(beta\) varierar över tid för varje enskild prediktor. Det skall inte finnas någon lutning eller systematisk trend bort från den horisontella linjen. Om så är fallet är effekten av \(beta\) konstant över tid och antagandet om proportional hazard är uppfyllt.

Om ett variabel strider mot antagandet om PH kan detta åtgärdas genom något av följande:

• Inkludera en interaktion mellan variabeln och tid i modellen.
• Stratifiera på den variabeln som strider mot PH. Nackdelen med detta är att man inte får en hazard ratio för en variabel som används för att stratifiera. Detta är dock bara ett problem om man är intresserad av hazard ratio för just den variabeln.

Stratifierad Cox regression

Cox regressionen erbjuder möjligheten att justera för variabler som inte ingår i modellen och detta görs genom att stratifiera. Man får inte förväxla stratifiering vid Cox regression med stratifiering vid vanlig linjär (multipel) regression.

För vanlig linjär (multipel) regression så innebär stratifiering att man skapar separata modeller utifrån en kategorisk variabel. Ponera att vi studerar hur ålder påverkar blodtryck och vi har även information om vikt. Om vi vill veta hur ålder påverkar blodtryck i olika viktkategorier så kan vi stratifiera analysen genom att göra en regression för normalviktiga, en regression för överviktiga och en regression för personer med fetma. Detta är alltså tre separata modeller som ger separata koefficienter för ålder.

För Cox regression kan stratifieringen hanteras annorlunda. Cox regression tillåter oss stratifiera data (observationerna) utifrån en kategorisk variabel. Baseline hazard tillåts variera mellan dessa strata. Modellen kan sedan estimera en sammanlagd (kumulativ) effekt av ålder i alla strata. Då erhåller vi endast en koefficient för ålder, istället för separata koefficienter i varje strata. Man kan säga att effekten av ålder “poolas” (sammanvägs) i alla strata.

Det finns flera situationer då det är klokt att statifiera modellen. Om en variabel strider mot antagandet om proportional hazard så kan man kategorisera variabeln till två eller flera lämpliga nivåer och stratifiera modellen på den variabeln. Då har man justerat för variabeln och man behöver inte oroa sig för att den strider mot proportional hazards. Nackdelen är att man inte erhåller några regressionskoefficienter för variabeln (men om den ändå inte är av intresse för analysen så är detta ingen större nackdel.

Med föjande kod stratifieras Cox regressionen på kön (sex):

min_modell <- coxph(Surv(time, status==1) ~ strata(sex) + age + nodes + perfor + obstruct, data=colon)

Vi skriver ut resultatet med funktionen tbl_regression(), enligt föjande:

BLOCKERAD KOD, BILD ELLER TEXT


Du måste logga in för att komma åt allt material.

Som framgår av ovanstående har stratifieringen inte påverkat hazard kvoterna för övriga prediktorer i modellen.

Referenser

Germán Rodríguez, Survival Models, Princeton university.

Cox DR and Oakes D (1984). Analysis of Survival Data. Chapman and Hall, New York.

Kalbfleisch JD and Prentice RL (2002). The Statistical Analysis of Failure Time Data. Second Edition. Wiley, New York.

Klein and Moeschberger (1997). Survival Analysis: Techniques for Censored and Truncated Data. Springer-Verlag, New York.

Lagakos SW (1979). General right censoring and its impact on the analysis of survival data, Biometrics 35:139-156.

Lawless JF (2003). Statistical Models and Methods for Lifetime Data. Wiley, New York.

Clark, T., Bradburn, M., Love, S., & Altman, D. (2003). Survival analysis part I: Basic concepts and first analyses. 232-238. ISSN 0007-0920.

M J Bradburn, T G Clark, S B Love, & D G Altman. (2003). Survival Analysis Part II: Multivariate data analysis – an introduction to concepts and methods. British Journal of Cancer, 89(3), 431-436.

Bradburn, M., Clark, T., Love, S., & Altman, D. (2003). Survival analysis Part III: Multivariate data analysis – choosing a model and assessing its adequacy and fit. 89(4), 605-11.

Clark, T., Bradburn, M., Love, S., & Altman, D. (2003). Survival analysis part IV: Further concepts and methods in survival analysis. 781-786. ISSN 0007-0920.