# Survival data

## Censoring

• Type I censoring: the censoring time is fixed
• Type II censoring
• Random censoring
• Right censoring
• Left censoring
• Interval censoring
• Truncation

The most common is called right censoring and occurs when a participant does not have the event of interest during the study and thus their last observed follow-up time is less than their time to event. This can occur when a participant drops out before the study ends or when a participant is event free at the end of the observation period.

• Event: Death, disease occurrence, disease recurrence, recovery, or other experience of interest
• Time: The time from the beginning of an observation period (such as surgery or beginning treatment) to (i) an event, or (ii) end of the study, or (iii) loss of contact or withdrawal from the study.
• Censoring / Censored observation: If a subject does not have an event during the observation time, they are described as censored. The subject is censored in the sense that nothing is observed or known about that subject after the time of censoring. A censored subject may or may not have an event after the end of observation time.

In R, "status" should be called event status. status = 1 means event occurred. status = 0 means no event (censored). Sometimes the status variable has more than 2 states. We can uses "status != 0" to replace "status" in Surv() function.

• status=0/1/2 for censored, transplant and dead in survival::pbc data.
• status=0/1/2 for censored, relapse and dead in randomForestSRC::follic data.

## How to explore survival data

• Create graph of length of time that each subject was in the study
library(survival)
# sort the aml data by time
aml <- aml[order(amltime),] with(aml, plot(time, type="h"))  • Create the life table survival object summary(aml.survfit) Call: survfit(formula = Surv(time, status == 1) ~ 1, data = aml) time n.risk n.event survival std.err lower 95% CI upper 95% CI 5 23 2 0.9130 0.0588 0.8049 1.000 8 21 2 0.8261 0.0790 0.6848 0.996 9 19 1 0.7826 0.0860 0.6310 0.971 12 18 1 0.7391 0.0916 0.5798 0.942 13 17 1 0.6957 0.0959 0.5309 0.912 18 14 1 0.6460 0.1011 0.4753 0.878 23 13 2 0.5466 0.1073 0.3721 0.803 27 11 1 0.4969 0.1084 0.3240 0.762 30 9 1 0.4417 0.1095 0.2717 0.718 31 8 1 0.3865 0.1089 0.2225 0.671 33 7 1 0.3313 0.1064 0.1765 0.622 34 6 1 0.2761 0.1020 0.1338 0.569 43 5 1 0.2208 0.0954 0.0947 0.515 45 4 1 0.1656 0.0860 0.0598 0.458 48 2 1 0.0828 0.0727 0.0148 0.462  • Kaplan-Meier curve for aml with the confidence bounds. plot(aml.survfit, xlab = "Time", ylab="Proportion surviving")  • Create aml life tables broken out by treatment (x, "Maintained" vs. "Not maintained") surv.by.aml.rx <- survfit(Surv(time, status == 1) ~ x, data = aml) summary(surv.by.aml.rx) Call: survfit(formula = Surv(time, status == 1) ~ x, data = aml) x=Maintained time n.risk n.event survival std.err lower 95% CI upper 95% CI 9 11 1 0.909 0.0867 0.7541 1.000 13 10 1 0.818 0.1163 0.6192 1.000 18 8 1 0.716 0.1397 0.4884 1.000 23 7 1 0.614 0.1526 0.3769 0.999 31 5 1 0.491 0.1642 0.2549 0.946 34 4 1 0.368 0.1627 0.1549 0.875 48 2 1 0.184 0.1535 0.0359 0.944 x=Nonmaintained time n.risk n.event survival std.err lower 95% CI upper 95% CI 5 12 2 0.8333 0.1076 0.6470 1.000 8 10 2 0.6667 0.1361 0.4468 0.995 12 8 1 0.5833 0.1423 0.3616 0.941 23 6 1 0.4861 0.1481 0.2675 0.883 27 5 1 0.3889 0.1470 0.1854 0.816 30 4 1 0.2917 0.1387 0.1148 0.741 33 3 1 0.1944 0.1219 0.0569 0.664 43 2 1 0.0972 0.0919 0.0153 0.620 45 1 1 0.0000 NaN NA NA  • Plot KM plot broken out by treatment plot(surv.by.aml.rx, xlab = "Time", ylab="Survival", col=c("black", "red"), lty = 1:2, main="Kaplan-Meier Survival vs. Maintenance in AML") legend(100, .6, c("Maintained", "Not maintained"), lty = 1:2, col=c("black", "red"))  • Perform the log rank test using the R function survdiff(). surv.diff.aml <- survdiff(Surv(time, status == 1) ~ x, data=aml) surv.diff.aml Call: survdiff(formula = Surv(time, status == 1) ~ x, data = aml) N Observed Expected (O-E)^2/E (O-E)^2/V x=Maintained 11 7 10.69 1.27 3.4 x=Nonmaintained 12 11 7.31 1.86 3.4 Chisq= 3.4 on 1 degrees of freedom, p= 0.07  ### Summary statistics • Kaplan-Meier Method and Log-Rank Test • Statistics • Table of status vs treatment (with proportion) • Table of treatment vs training/test • Life table • summary(survfit(Surv(time, status) ~ 1)) or summary(survfit(Surv(time, status) ~ treatment)) • KMsurv::lifetab() • Coxph function and visualize them using the ggforest package ### Some public data package data (sample size) survival pbc (418), ovarian (26), aml/leukemia (23), colon (1858), lung (228), veteran (137) pec GBSG2 (686), cost (518) randomForestSRC follic (541) KMsurv A LOT. tongue (80) survivalROC mayo (312) survAUC NA ## Kaplan & Meier and Nelson-Aalen: survfit.formula(), Surv() • Landmarks • Kaplan-Meier: 1958 • Nelson: 1969 • Cox and Brewlow: 1972 S(t) = exp(-Lambda(t)) • Aalen: 1978 Lambda(t) • D distinct times ${\displaystyle t_{1}. At time ${\displaystyle t_{i}}$ there are ${\displaystyle d_{i}}$ events. Let ${\displaystyle Y_{i}}$ be the number of individuals who are at risk at time ${\displaystyle t_{i}}$. The quantity ${\displaystyle d_{i}/Y_{i}}$ provides an estimate of the conditional probability that an individual who survives to just prior to time ${\displaystyle t_{i}}$ experiences the event at time ${\displaystyle t_{i}}$. The KM estimator of the survival function and the Nelson-Aalen estimator of the cumulative hazard (their relationship is given below) are define as follows (${\displaystyle t_{1}\leq t}$): {\displaystyle {\begin{aligned}{\hat {S}}(t)&=\prod _{t_{i}\leq t}[1-d_{i}/Y_{i}]\\{\hat {H}}(t)&=\sum _{t_{i}\leq t}d_{i}/Y_{i}\end{aligned}}} str(kidney) 'data.frame': 76 obs. of 7 variables: id     : num  1 1 2 2 3 3 4 4 5 5 ...
$time : num 8 16 23 13 22 28 447 318 30 12 ...$ status : num  1 1 1 0 1 1 1 1 1 1 ...
$age : num 28 28 48 48 32 32 31 32 10 10 ...$ sex    : num  1 1 2 2 1 1 2 2 1 1 ...
$disease: Factor w/ 4 levels "Other","GN","AN",..: 1 1 2 2 1 1 1 1 1 1 ...$ frail  : num  2.3 2.3 1.9 1.9 1.2 1.2 0.5 0.5 1.5 1.5 ...
kidney[order(kidney$time), c("time", "status")] kidney[kidney$time == 13, ] # one is dead and the other is alive
length(unique(kidney$time)) # 60 sfit <- survfit(Surv(time, status) ~ 1, data = kidney) sfit Call: survfit(formula = Surv(time, status) ~ 1, data = kidney) n events median 0.95LCL 0.95UCL 76 58 78 39 152 str(sfit) List of 13$ n        : int 76
$time : num [1:60] 2 4 5 6 7 8 9 12 13 15 ...$ n.risk   : num [1:60] 76 75 74 72 71 69 65 64 62 60 ...
$n.event : num [1:60] 1 0 0 0 2 2 1 2 1 2 ...$ n.censor : num [1:60] 0 1 2 1 0 2 0 0 1 0 ...
$surv : num [1:60] 0.987 0.987 0.987 0.987 0.959 ...$ type     : chr "right"
length(unique(kidney$time)) # [1] 60 all(sapply(sfit$time, function(tt) sum(kidney$time >= tt)) == sfit$n.risk) # TRUE
all(sapply(sfit$time, function(tt) sum(kidney$status[kidney$time == tt])) == sfit$n.event) # TRUE
all(sapply(sfit$time, function(tt) sum(1-kidney$status[kidney$time == tt])) == sfit$n.censor) #  TRUE
all(cumprod(1 - sfit$n.event/sfit$n.risk) == sfit$surv) # FALSE range(abs(cumprod(1 - sfit$n.event/sfit$n.risk) - sfit$surv))
# [1] 0.000000e+00 1.387779e-17

summary(sfit)
time n.risk n.event survival std.err lower 95% CI upper 95% CI
2     76       1    0.987  0.0131      0.96155        1.000
7     71       2    0.959  0.0232      0.91469        1.000
8     69       2    0.931  0.0297      0.87484        0.991
...
511      3       1    0.042  0.0288      0.01095        0.161
536      2       1    0.021  0.0207      0.00305        0.145
562      1       1    0.000     NaN           NA           NA

• Note that the KM estimate is left-continuous step function with the intervals closed at left and open at right. For ${\displaystyle t\in [t_{j},t_{j+1})}$ for a certain j, we have ${\displaystyle {\hat {S}}(t)=\prod _{i=1}^{j}(1-d_{i}/n_{i})}$ where ${\displaystyle d_{i}}$ is the number people who have an event during the interval ${\displaystyle [t_{i},t_{i+1})}$ and ${\displaystyle n_{i}}$ is the number of people at risk just before the beginning of the interval ${\displaystyle [t_{i},t_{i+1})}$.
• The product-limit estimator can be constructed by using a reduced-sample approach. We can estimate the ${\displaystyle P(T>t_{i}|T\geq t_{i})={\frac {Y_{i}-d_{i}}{Y_{i}}}}$ for ${\displaystyle i=1,2,\cdots ,D}$. ${\displaystyle S(t_{i})={\frac {S(t_{i})}{S(t_{i-1})}}{\frac {S(t_{i-1})}{S(t_{i-2})}}\cdots {\frac {S(t_{2})}{S(t_{1})}}{\frac {S(t_{1})}{S(0)}}S(0)=P(T>t_{i}|T\geq t_{i})P(T>t_{i-1}|T\geq t_{i-1})\cdots P(T>t_{2}|T\geq t_{2})P(T>t_{1}|T\geq t_{1})}$ because S(0)=1 and, for a discrete distribution, ${\displaystyle S(t_{i-1})=P(T>t_{i-1})=P(T\geq t_{i})}$.
• Self consistency. If we had no censored observations, the estimator of the survival function at a time t is the proportion of observations which are larger than t, that is, ${\displaystyle {\hat {S}}(t)={\frac {1}{n}}\sum I(X_{i}>t)}$.
• Curves are plotted in the same order as they are listed by print (which gives a 1 line summary of each). For example, -1 < 1 and 'Maintenance' < 'Nonmaintained'. That means, the labels list in the legend() command should have the same order as the curves.
• Kaplan and Meier is used to give an estimator of the survival function S(t)
• Nelson-Aalen estimator is for the cumulative hazard H(t). Note that ${\displaystyle 0\leq H(t)<\infty }$ and ${\displaystyle H(t)\rightarrow \infty }$ as t goes to infinity. So there is a constraint on the hazard function, see Wikipedia.

Note that S(t) is related to H(t) by ${\displaystyle H(t)=-ln[S(t)]}$ or ${\displaystyle S(t)=exp[-H(t)]}$. The two estimators are similar (see example 4.1A and 4.1B from Klein and Moeschberge).

The Nelson-Aalen estimator has two primary uses in analyzing data

1. Selecting between parametric models for the time to event
2. Crude estimates of the hazard rate h(t). This is related to the estimation of the survival function in Cox model. See 8.6 of Klein and Moeschberge.

The Kaplan–Meier estimator (the product limit estimator) is an estimator for estimating the survival function from lifetime data. In medical research, it is often used to measure the fraction of patients living for a certain amount of time after treatment.

Note that

• The "+" sign in the KM curves means censored observations (this convention matches with the output of Surv() function) and a long vertical line (not '+') means there is a dead observation at that time.
> aml[1:5,]
time status          x
1    9      1 Maintained
2   13      1 Maintained
3   13      0 Maintained
4   18      1 Maintained
5   23      1 Maintained
> Surv(aml$time, aml$status)[1:5,]
[1]  9  13  13+ 18  23

• If the last observation (longest survival time) is dead, the survival curve will goes down to zero. Otherwise, the survival curve will remain flat from the last event time.

Usually the KM curve of treatment group is higher than that of the control group.

The Y-axis (the probability that a member from a given population will have a lifetime exceeding time) is often called

• Cumulative probability
• Cumulative survival
• Percent survival
• Probability without event
• Proportion alive/surviving
• Survival
• Survival probability
> library(survival)
> str(aml$x) Factor w/ 2 levels "Maintained","Nonmaintained": 1 1 1 1 1 1 1 1 1 1 ... > plot(leukemia.surv <- survfit(Surv(time, status) ~ x, data = aml[7:17,] ) , lty=2:3, mark.time = TRUE) # a (small) subset, mark.time is used to show censored obs > aml[7:17,] time status x 7 31 1 Maintained 8 34 1 Maintained 9 45 0 Maintained 10 48 1 Maintained 11 161 0 Maintained 12 5 1 Nonmaintained 13 5 1 Nonmaintained 14 8 1 Nonmaintained 15 8 1 Nonmaintained 16 12 1 Nonmaintained 17 16 0 Nonmaintained > legend(100, .9, c("Maintenance", "No Maintenance"), lty = 2:3) # lty: 2=dashed, 3=dotted > title("Kaplan-Meier Curves\nfor AML Maintenance Study") # Cumulative hazard plot # Lambda(t) = -log(S(t)); # see https://en.wikipedia.org/wiki/Survival_analysis # http://statweb.stanford.edu/~olshen/hrp262spring01/spring01Handouts/Phil_doc.pdf plot(leukemia.surv <- survfit(Surv(time, status) ~ x, data = aml[7:17,] ) , lty=2:3, mark.time = T, fun="cumhaz", ylab="Cumulative Hazard")  # https://www.lexjansen.com/pharmasug/2011/CC/PharmaSUG-2011-CC16.pdf mydata <- data.frame(time=c(3,6,8,12,12,21),status=c(1,1,0,1,1,1)) km <- survfit(Surv(time, status)~1, data=mydata) plot(km, mark.time = T) survest <- stepfun(km$time, c(1, km$surv)) plot(survest) > str(km) List of 13$ n        : int 6
$time : num [1:5] 3 6 8 12 21$ n.risk   : num [1:5] 6 5 4 3 1
$n.event : num [1:5] 1 1 0 2 1$ n.censor : num [1:5] 0 0 1 0 0
$surv : num [1:5] 0.833 0.667 0.667 0.222 0$ type     : chr "right"
$std.err : num [1:5] 0.183 0.289 0.289 0.866 Inf$ upper    : num [1:5] 1 1 1 1 NA
$lower : num [1:5] 0.5827 0.3786 0.3786 0.0407 NA$ conf.type: chr "log"
$conf.int : num 0.95 > class(survest) [1] "stepfun" "function" > survest Step function Call: stepfun(km$time, c(1, km$surv)) x[1:5] = 3, 6, 8, 12, 21 6 plateau levels = 1, 0.83333, 0.66667, ..., 0.22222, 0 > str(survest) function (v) - attr(*, "class")= chr [1:2] "stepfun" "function" - attr(*, "call")= language stepfun(km$time, c(1, kmsurv))  ### Multiple curves Curves/groups are ordered. The first color in the palette is used to color the first level of the factor variable. This is same idea as ggsurvplot in the survminer package. This affects parameters like col and lty in plot() function. For example, • 1<2 • 'c' < 't' • 'control' < 'treatment' • 'Control' < 'Treatment' • 'female' < 'male'. For legend(), the first category in legend argument will appear at the top of the legend box. ### Inverse Probability of Censoring Weighted (IPCW) The plots below show by flipping the status variable, we can accurately recover the survival function of the censoring variable. See the R code here for superimposing the true exponential distribution on the KM plot of the censoring variable. require(survival) n = 10000 beta1 = 2; beta2 = -1 lambdaT = 1 # baseline hazard lambdaC = 2 # hazard of censoring set.seed(1234) x1 = rnorm(n,0) x2 = rnorm(n,0) # true event time # T = rweibull(n, shape=1, scale=lambdaT*exp(-beta1*x1-beta2*x2)) # Wrong T = Vectorize(rweibull)(n=1, shape=1, scale=lambdaT*exp(-beta1*x1-beta2*x2)) # method 1: exponential censoring variable C <- rweibull(n, shape=1, scale=lambdaC) time = pmin(T,C) status <- 1*(T <= C) mean(status) summary(T) summary(C) par(mfrow=c(2,1), mar = c(3,4,2,2)+.1) status2 <- 1-status plot(survfit(Surv(time, status2) ~ 1), ylab="Survival probability", main = 'Exponential censoring time') # method 2: uniform censoring variable C <- runif(n, 0, 21) time = pmin(T,C) status <- 1*(T <= C) status2 <- 1-status plot(survfit(Surv(time, status2) ~ 1), ylab="Survival probability", main = "Uniform censoring time")  ### stepfun() and plot.stepfun() ### Survival curves with number at risk at bottom: survminer package R function survminer::ggsurvplot() Paper examples ### Life table ## Alternatives to survival function plot https://www.rdocumentation.org/packages/survival/versions/2.43-1/topics/plot.survfit The fun argument, a transformation of the survival curve • fun = "event" or "F": f(y) = 1-y; it calculates P(T < t). This is like a t-year risk (Blanche 2018). • fun = "cumhaz": cumulative hazard function (f(y) = -log(y)); it calculates H(t). See Intuition for cumulative hazard function. ## Breslow estimate ## Logrank test ## Survival curve with confidence interval ## Parametric models and survival function for censored data Assume the CDF of survival time T is ${\displaystyle F(\cdot )}$ and the CDF of the censoring time C is ${\displaystyle G(\cdot )}$, {\displaystyle {\begin{aligned}P(T>t,\delta =1)&=\int _{t}^{\infty }(1-G(s))dF(s),\\P(T>t,\delta =0)&=\int _{t}^{\infty }(1-F(s))dG(s)\end{aligned}}} ### R ## Parametric models and likelihood function for uncensored data • Exponential. ${\displaystyle T\sim Exp(\lambda )}$. ${\displaystyle H(t)=\lambda t.}$ and ${\displaystyle ln(S(t))=-H(t)=-\lambda t.}$ • Weibull. ${\displaystyle T\sim W(\lambda ,p).}$ ${\displaystyle H(t)=\lambda ^{p}t^{p}.}$ and ${\displaystyle ln(-ln(S(t)))=ln(\lambda ^{p}t^{p})=const+pln(t)}$. See also accelerated life models where a set of covariates were used to model survival time. ## Survival modeling ### Accelerated life models - a direct extension of the classical linear model http://data.princeton.edu/wws509/notes/c7.pdf and also Kalbfleish and Prentice (1980). ${\displaystyle logT_{i}=x_{i}'\beta +\epsilon _{i}}$ Therefore • ${\displaystyle T_{i}=exp(x_{i}'\beta )T_{0i}}$. So if there are two groups (x=1 and x=0), and ${\displaystyle exp(\beta )=2}$, it means one group live twice as long as people in another group. • ${\displaystyle S_{1}(t)=S_{0}(t/exp(x'\beta ))}$. This explains the meaning of accelerated failure-time. Depending on the sign of ${\displaystyle \beta 'x}$, the time is either accelerated by a constant factor or degraded by a constant factor. If ${\displaystyle exp(\beta )=2}$, the probability that a member in group one (eg treatment) will be alive at age t is exactly the same as the probability that a member in group zero (eg control group) will be alive at age t/2. • The hazard function ${\displaystyle \lambda _{1}(t)=\lambda _{0}(t/exp(x'\beta ))/exp(x'\beta )}$. So if ${\displaystyle exp(\beta )=2}$, at any given age people in group one would be exposed to half the risk of people in group zero half their age. In applications, • If the errors are normally distributed, then we obtain a log-normal model for the T. Estimation of this model for censored data by maximum likelihood is known in the econometric literature as a Tobit model. • If the errors have an extreme value distribution, then T has an exponential distribution. The hazard ${\displaystyle \lambda }$ satisfies the log linear model ${\displaystyle \log \lambda _{i}=x_{i}'\beta }$. ### Proportional hazard models Note PH models is a type of multiplicative hazard rate models ${\displaystyle h(x|Z)=h_{0}(x)c(\beta 'Z)}$ where ${\displaystyle c(\beta 'Z)=\exp(\beta 'Z)}$. Assumption: Survival curves for two strata (determined by the particular choices of values for covariates) must have hazard functions that are proportional over time (i.e. constant relative hazard over time). Proportional hazards assumption meaning. The ratio of the hazard rates from two individuals with covariate value ${\displaystyle Z}$ and ${\displaystyle Z^{*}}$ is a constant function time. {\displaystyle {\begin{aligned}{\frac {h(t|Z)}{h(t|Z^{*})}}={\frac {h_{0}(t)\exp(\beta 'Z)}{h_{0}(t)\exp(\beta 'Z^{*})}}=\exp(\beta '(Z-Z^{*})){\mbox{ independent of time}}\end{aligned}}} Test the assumption ## Weibull and Exponential model to Cox model In summary: • Weibull distribution (Klein) ${\displaystyle h(t)=p\lambda (\lambda t)^{p-1}}$ and ${\displaystyle S(t)=exp(-\lambda t^{p})}$. If p >1, then the risk increases over time. If p<1, then the risk decreases over time. • Note that Weibull distribution has a different parametrization. See http://data.princeton.edu/pop509/ParametricSurvival.pdf#page=2. ${\displaystyle h(t)=\lambda ^{p}pt^{p-1}}$ and ${\displaystyle S(t)=exp(-(\lambda t)^{p})}$. R and wikipedia also follows this parametrization except that ${\displaystyle h(t)=pt^{p-1}/\lambda ^{p}}$ and ${\displaystyle S(t)=exp(-(t/\lambda )^{p})}$. • Exponential distribution ${\displaystyle h(t)}$ = constant (independent of t). This is a special case of Weibull distribution (p=1). • Weibull (and also exponential) distribution regression model is the only case which belongs to both the proportional hazards and the accelerated life families. {\displaystyle {\begin{aligned}{\frac {h(x|Z_{1})}{h(x|Z_{2})}}={\frac {h_{0}(x\exp(-\gamma 'Z_{1}))\exp(-\gamma 'Z_{1})}{h_{0}(x\exp(-\gamma 'Z_{2}))\exp(-\gamma 'Z_{2})}}={\frac {(a/b)\left({\frac {x\exp(-\gamma 'Z_{1})}{b}}\right)^{a-1}\exp(-\gamma 'Z_{1})}{(a/b)\left({\frac {x\exp(-\gamma 'Z_{2})}{b}}\right)^{a-1}\exp(-\gamma 'Z_{2})}}\quad {\mbox{which is independent of time x}}\end{aligned}}} f(t)=h(t)*S(t) h(t) S(t) Mean Exponential (Klein p37) ${\displaystyle \lambda \exp(-\lambda t)}$ ${\displaystyle \lambda }$ ${\displaystyle \exp(-\lambda t)}$ ${\displaystyle 1/\lambda }$ Weibull (Klein, Bender, wikipedia) ${\displaystyle p\lambda t^{p-1}\exp(-\lambda t^{p})}$ ${\displaystyle p\lambda t^{p-1}}$ ${\displaystyle exp(-\lambda t^{p})}$ ${\displaystyle {\frac {\Gamma (1+1/p)}{\lambda ^{1/p}}}}$ Exponential (R) ${\displaystyle \lambda \exp(-\lambda t)}$, ${\displaystyle \lambda }$ is rate ${\displaystyle \lambda }$ ${\displaystyle \exp(-\lambda t)}$ ${\displaystyle 1/\lambda }$ Weibull (R, wikipedia) ${\displaystyle {\frac {a}{b}}\left({\frac {t}{b}}\right)^{a-1}\exp(-({\frac {t}{b}})^{a})}$, ${\displaystyle a}$ is shape, and ${\displaystyle b}$ is scale ${\displaystyle {\frac {a}{b}}\left({\frac {t}{b}}\right)^{a-1}}$ ${\displaystyle \exp(-({\frac {t}{b}})^{a})}$ ${\displaystyle b\Gamma (1+1/a)}$ • Accelerated failure-time model. Let ${\displaystyle Y=\log(T)=\mu +\gamma 'Z+\sigma W}$. Then the survival function of ${\displaystyle T}$ at the covariate Z, {\displaystyle {\begin{aligned}S_{T}(t|Z)&=P(T>t|Z)\\&=P(Y>\ln t|Z)\\&=P(\mu +\sigma W>\ln t-\gamma 'Z|Z)\\&=P(e^{\mu +\sigma W}>t\exp(-\gamma 'Z)|Z)\\&=S_{0}(t\exp(-\gamma 'Z)).\end{aligned}}} where ${\displaystyle S_{0}(t)}$ denote the survival function T when Z=0. Since ${\displaystyle h(t)=-\partial \ln(S(t))}$, the hazard function of T with a covariate value Z is related to a baseline hazard rate ${\displaystyle h_{0}}$ by (p56 Klein) {\displaystyle {\begin{aligned}h(t|Z)=h_{0}(t\exp(-\gamma 'Z))\exp(-\gamma 'Z)\end{aligned}}} > mean(rexp(1000)^(1/2)) [1] 0.8902948 > mean(rweibull(1000, 2, 1)) [1] 0.8856265 > mean((rweibull(1000, 2, scale=4)/4)^2) [1] 1.008923  ### Graphical way to check Weibull, AFT, PH ### Weibull is related to Extreme value distribution ### Weibull distribution and bathtub ## CDF follows Unif(0,1) Take the Exponential distribution for example stem(pexp(rexp(1000))) stem(pexp(rexp(10000)))  Another example is from simulating survival time. Note that this is exactly Bender et al 2005 approach. See also the simsurv (newer) and survsim (older) packages. set.seed(100) #Define the following parameters outlined in the step: n = 1000 beta_0 = 0.5 beta_1 = -1 beta_2 = 1 b = 1.6 #This will be changed later as mentioned in Step 5 of documentation #Step 1 x_1<-rbinom(n, 1, 0.25) x_2<-rbinom(n, 1, 0.7) #Step 2 U<-runif(n, 0,1) T<-(-log(U)*exp(-(beta_0+beta_1*x_1+beta_2*x_2))) #Eqn (5) Fn <- ecdf(T) # https://stat.ethz.ch/R-manual/R-devel/library/stats/html/ecdf.html # verify F(T) or 1-F(T) ~ U(0, 1) hist(Fn(T)) # look at the plot of survival probability vs time plot(T, 1 - Fn(T))  ## Simulate survival data Note that status = 1 means an event (e.g. death) happened; Ti <= Ci. That is, the status variable used in R/Splus means the death indicator. y <- rexp(10) cen <- runif(10) status <- ifelse(cen < .7, 1, 0)  Note that for the exponential distribution, larger rate/${\displaystyle \lambda }$ corresponds to a smaller mean. This relation matches with the Cox regression where a large covariate corresponds to a smaller survival time. So the coefficient 3 in myrates in the below example has the same sign as the coefficient (2.457466 for censored data) in the output of the Cox model fitting. n <- 30 x <- scale(1:n, TRUE, TRUE) # create covariates (standardized) # the original example does not work on large 'n' myrates <- exp(3*x+1) set.seed(1234) y <- rexp(n, rate = myrates) # generates the r.v. cen <- rexp(n, rate = 0.5 ) # E(cen)=1/rate ycen <- pmin(y, cen) di <- as.numeric(y <= cen) survreg(Surv(ycen, di)~x, dist="weibull")coef[2]  # -3.080125
coxph(Surv(ycen, di)~x)$coef # 2.457466 # no censor survreg(Surv(y,rep(1, n))~x,dist="weibull")$coef[2]  # -3.137603
survreg(Surv(y,rep(1, n))~x,dist="exponential")$coef[2] # -3.143095 coxph(Surv(y,rep(1, n))~x)$coef  # 2.717794

# See the pdf note for the rest of code

{\displaystyle {\begin{aligned}\lambda =exp(-intercept)\end{aligned}}}
> futime <- rexp(1000, 5)
> survreg(Surv(futime,rep(1,1000))~1,dist="exponential")coef (Intercept) -1.618263 > exp(1.618263) [1] 5.044321  {\displaystyle {\begin{aligned}\gamma &=1/scale\\\alpha &=exp(-(Intercept)*\gamma )\end{aligned}}} > survreg(Surv(futime,rep(1,1000))~1,dist="weibull") Call: survreg(formula = Surv(futime, rep(1, 1000)) ~ 1, dist = "weibull") Coefficients: (Intercept) -1.639469 Scale= 1.048049 Loglik(model)= 620.1 Loglik(intercept only)= 620.1 n= 1000  • rsurv() function from the ipred package • Use Weibull distribution to model survival data. We assume the shape is constant across subjects. We then allow the scale to vary across subjects. For subject ${\displaystyle i}$ with covariate ${\displaystyle X_{i}}$, ${\displaystyle \log(scale_{i})}$ = ${\displaystyle \beta 'X_{i}}$. Note that if we want to make the ${\displaystyle \beta }$ sign to be consistent with the Cox model, we want to use ${\displaystyle \log(scale_{i})}$ = ${\displaystyle -\beta 'X_{i}}$ instead. • http://sas-and-r.blogspot.com/2010/03/example-730-simulate-censored-survival.html. Assuming shape=1 in the Weibull distribution, then the hazard function can be expressed as a proportional hazard model ${\displaystyle h(t|x)=1/scale={\frac {1}{\lambda /e^{\beta 'x}}}={\frac {e^{\beta 'x}}{\lambda }}=h_{0}(t)\exp(\beta 'x)}$ n = 10000 beta1 = 2; beta2 = -1 lambdaT = .002 # baseline hazard lambdaC = .004 # hazard of censoring set.seed(1234) x1 = rnorm(n,0) x2 = rnorm(n,0) # true event time T = Vectorize(rweibull)(n=1, shape=1, scale=lambdaT*exp(-beta1*x1-beta2*x2)) # No censoring event2 <- rep(1, length(T)) coxph(Surv(T, event2)~ x1 + x2) # coef exp(coef) se(coef) z p # x1 1.99825 7.37613 0.01884 106.07 <2e-16 # x2 -1.00200 0.36715 0.01267 -79.08 <2e-16 # # Likelihood ratio test=15556 on 2 df, p=< 2.2e-16 # n= 10000, number of events= 10000 # Censoring C = rweibull(n, shape=1, scale=lambdaC) #censoring time time = pmin(T,C) #observed time is min of censored and true event = time==T # set to 1 if event is observed coxph(Surv(time, event)~ x1 + x2) # coef exp(coef) se(coef) z p # x1 2.01039 7.46622 0.02250 89.33 <2e-16 # x2 -0.99210 0.37080 0.01552 -63.95 <2e-16 # # Likelihood ratio test=11321 on 2 df, p=< 2.2e-16 # n= 10000, number of events= 6002 mean(event) # [1] 0.6002  • https://stats.stackexchange.com/a/135129 (Bender's inverse probability method). Let ${\displaystyle h_{0}(t)=\lambda \rho t^{\rho -1}}$ where shape 𝜌>0 and scale 𝜆>0. Following the inverse probability method, a realisation of 𝑇∼𝑆(⋅|𝐱) is obtained by computing ${\displaystyle t=\left(-{\frac {\log(v)}{\lambda \exp(x'\beta )}}\right)^{1/\rho }}$ with 𝑣 a uniform variate on (0,1). Using results on transformations of random variables, one may notice that 𝑇 has a conditional Weibull distribution (given 𝐱) with shape 𝜌 and scale 𝜆exp(𝐱′𝛽). # N = sample size # lambda = scale parameter in h0() # rho = shape parameter in h0() # beta = fixed effect parameter # rateC = rate parameter of the exponential distribution of censoring variable C simulWeib <- function(N, lambda, rho, beta, rateC) { # covariate --> N Bernoulli trials x <- sample(x=c(0, 1), size=N, replace=TRUE, prob=c(0.5, 0.5)) # Weibull latent event times v <- runif(n=N) Tlat <- (- log(v) / (lambda * exp(x * beta)))^(1 / rho) # censoring times C <- rexp(n=N, rate=rateC) # follow-up times and event indicators time <- pmin(Tlat, C) status <- as.numeric(Tlat <= C) # data set data.frame(id=1:N, time=time, status=status, x=x) } # Test set.seed(1234) betaHat <- rate <- rep(NA, 1e3) for(k in 1:1e3) { dat <- simulWeib(N=100, lambda=0.01, rho=1, beta=-0.6, rateC=0.001) fit <- coxph(Surv(time, status) ~ x, data=dat) rate[k] <- mean(datstatus == 0)
betaHat[k] <- fit$coef } mean(rate) # [1] 0.12287 mean(betaHat) # [1] -0.6085473  ## Standardize covariates coxph() does not have an option to standardize covariates but glmnet() does. library(glmnet) library(survival) N=1000;p=30 nzc=p/3 beta <- c(rep(1, 5), rep(-1, 5)) set.seed(1234) x=matrix(rnorm(N*p),N,p) x[, 1:5] <- x[, 1:5]*2 x[, 6:10] <- x[, 6:10] + 2 fx=x[,seq(nzc)] %*% beta hx=exp(fx) ty=rexp(N,hx) tcens <- rep(0,N) y=cbind(time=ty,status=1-tcens) # y=Surv(ty,1-tcens) with library(survival) coxph(Surv(ty, 1-tcens) ~ x) %>% coef %>% head(10) # x1 x2 x3 x4 x5 x6 x7 # 0.6076146 0.6359927 0.6346022 0.6469274 0.6152082 -0.6614930 -0.5946101 # x8 x9 x10 # -0.6726081 -0.6275205 -0.7073704 xscale <- scale(x, TRUE, TRUE) # halve the covariate values coxph(Surv(ty, 1-tcens) ~ xscale) %>% coef %>% head(10) # double the coef # xscale1 xscale2 xscale3 xscale4 xscale5 xscale6 xscale7 # 1.2119940 1.2480628 1.2848646 1.2857796 1.1959619 -0.6431946 -0.5941309 # xscale8 xscale9 xscale10 # -0.6723137 -0.6188384 -0.6793313 set.seed(1) fit=cv.glmnet(x,y,family="cox", nfolds=10, standardize = TRUE) as.vector(coef(fit, s = "lambda.min"))[seq(nzc)] # [1] 0.9351341 0.9394696 0.9187242 0.9418540 0.9111623 -0.9303783 # [7] -0.9271438 -0.9597583 -0.9493759 -0.9386065 set.seed(1) fit=cv.glmnet(x,y,family="cox", nfolds=10, standardize = FALSE) as.vector(coef(fit, s = "lambda.min"))[seq(nzc)] # [1] 0.9357171 0.9396877 0.9200247 0.9420215 0.9118803 -0.9257406 # [7] -0.9232813 -0.9554017 -0.9448827 -0.9356009 set.seed(1) fit=cv.glmnet(xscale,y,family="cox", nfolds=10, standardize = TRUE) as.vector(coef(fit, s = "lambda.min"))[seq(nzc)] # [1] 1.8652889 1.8436015 1.8601198 1.8719515 1.7712951 -0.9046420 # [7] -0.9263966 -0.9593383 -0.9362407 -0.9014015 set.seed(1) fit=cv.glmnet(xscale,y,family="cox", nfolds=10, standardize = FALSE) as.vector(coef(fit, s = "lambda.min"))[seq(nzc)] # [1] 1.8652889 1.8436015 1.8601198 1.8719515 1.7712951 -0.9046420 # [7] -0.9263966 -0.9593383 -0.9362407 -0.9014015  ## Predefined censoring rates ## Cross validation • Cross validation in survival analysis by Verweij & van Houwelingen, Stat in medicine 1993. • Using cross-validation to evaluate predictive accuracy of survival risk classifiers based on high-dimensional data. Simon et al, Brief Bioinform. 2011 ## Competing risk ## Survival rate terminology • Disease-free survival (DFS): the period after curative treatment [disease eliminated] when no disease can be detected • Progression-free survival (PFS), overall survival (OS). PFS is the length of time during and after the treatment of a disease, such as cancer, that a patient lives with the disease but it does not get worse. See an use at NCI-MATCH trial. • Time to progression: The length of time from the date of diagnosis or the start of treatment for a disease until the disease starts to get worse or spread to other parts of the body. In a clinical trial, measuring the time to progression is one way to see how well a new treatment works. Also called TTP. • Metastasis-free survival (MFS) time: the period until metastasis is detected • Understanding Statistics Used to Guide Prognosis and Evaluate Treatment (DFS & PFS rate) ## Books ## HER2-positive breast cancer # Cox proportional hazards model and the partial log-likelihood function Let Yi denote the observed time (either censoring time or event time) for subject i, and let Ci be the indicator that the time corresponds to an event (i.e. if Ci = 1 the event occurred and if Ci = 0 the time is a censoring time). The hazard function for the Cox proportional hazard model has the form ${\displaystyle \lambda (t|X)=\lambda _{0}(t)\exp(\beta _{1}X_{1}+\cdots +\beta _{p}X_{p})=\lambda _{0}(t)\exp(X\beta ^{\prime }).}$ This expression gives the hazard at time t for an individual with covariate vector (explanatory variables) X. Based on this hazard function, a partial likelihood (defined on hazard function) can be constructed from the datasets as ${\displaystyle L(\beta )=\prod \limits _{i:C_{i}=1}{\frac {\theta _{i}}{\sum _{j:Y_{j}\geq Y_{i}}\theta _{j}}},}$ where θj = exp(Xj β) and X1, ..., Xn are the covariate vectors for the n independently sampled individuals in the dataset (treated here as column vectors). This pdf or this note give a toy example The corresponding log partial likelihood is ${\displaystyle \ell (\beta )=\sum _{i:C_{i}=1}\left(X_{i}\beta ^{\prime }-\log \sum _{j:Y_{j}\geq Y_{i}}\theta _{j}\right).}$ This function can be maximized over β to produce maximum partial likelihood estimates of the model parameters. The partial score function is ${\displaystyle \ell ^{\prime }(\beta )=\sum _{i:C_{i}=1}\left(X_{i}-{\frac {\sum _{j:Y_{j}\geq Y_{i}}\theta _{j}X_{j}}{\sum _{j:Y_{j}\geq Y_{i}}\theta _{j}}}\right),}$ and the Hessian matrix of the partial log likelihood is ${\displaystyle \ell ^{\prime \prime }(\beta )=-\sum _{i:C_{i}=1}\left({\frac {\sum _{j:Y_{j}\geq Y_{i}}\theta _{j}X_{j}X_{j}^{\prime }}{\sum _{j:Y_{j}\geq Y_{i}}\theta _{j}}}-{\frac {\sum _{j:Y_{j}\geq Y_{i}}\theta _{j}X_{j}\times \sum _{j:Y_{j}\geq Y_{i}}\theta _{j}X_{j}^{\prime }}{[\sum _{j:Y_{j}\geq Y_{i}}\theta _{j}]^{2}}}\right).}$ Using this score function and Hessian matrix, the partial likelihood can be maximized using the Newton-Raphson algorithm. The inverse of the Hessian matrix, evaluated at the estimate of β, can be used as an approximate variance-covariance matrix for the estimate, and used to produce approximate standard errors for the regression coefficients. If X is age, then the coefficient is likely >0. If X is some treatment, then the coefficient is likely <0. ## Compare the partial likelihood to the full likelihood ## z-column (Wald statistic) from R's coxph() ## How exactly can the Cox-model ignore exact times? library(survival) survfit(Surv(time, status) ~ x, data = aml) fit <- coxph(Surv(time, status) ~ x, data = aml) coef(fit) # 0.9155326 min(diff(sort(unique(aml$time)))) # 1

# Shift survival time for some obs but keeps the same order
# make sure we choose obs (n=20 not works but n=21 works) with twins
rbind(order(aml$time), sort(aml$time), aml$time[order(aml$time)])
# [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16]
# [1,]   12   13   14   15    1   16    2    3   17     4     5    18    19     6    20     7
# [2,]    5    5    8    8    9   12   13   13   16    18    23    23    27    28    30    31
# [3,]    5    5    8    8    9   12   13   13   16    18    23    23    27    28    30    31
# [,17] [,18] [,19] [,20] [,21] [,22] [,23]
# [1,]    21     8    22     9    23    10    11
# [2,]    33    34    43    45    45    48   161
# [3,]    33    34    43    45    45    48   161

aml$time2 <- aml$time
aml$time2[order(aml$time)[1:21]] <- aml$time[order(aml$time)[1:21]] - .9
fit2 <- coxph(Surv(time2, status) ~ x, data = aml); fit2
coef(fit2) #      0.9155326
coef(fit) == coef(fit2) # TRUE

aml$time3 <- aml$time
aml$time3[order(aml$time)[1:20]] <- aml$time[order(aml$time)[1:20]] - .9
fit3 <- coxph(Surv(time3, status) ~ x, data = aml); fit3
coef(fit3) #      0.8891567
coef(fit) == coef(fit3) # FALSE


## Partial likelihood when there are ties; hypothesis testing: Likelihood Ratio Test, Wald Test & Score Test

In R's coxph(): Nearly all Cox regression programs use the Breslow method by default, but not this one. The Efron approximation is used as the default here, it is more accurate when dealing with tied death times, and is as efficient computationally.

http://sfb649.wiwi.hu-berlin.de/fedc_homepage/xplore/tutorials/xaghtmlnode28.html (include the case when there is a partition of parameters). The formulas for 3 tests are also available on Appendix B of Klein book.

The following code does not test for models. But since there is only one coefficient, the results are the same. If there is more than one variable, we can use anova(model1, model2) to run LRT.

library(KMsurv)
# No ties. Section 8.2
data(btrial)
str(btrial)
# 'data.frame':	45 obs. of  3 variables:
# $time : int 19 25 30 34 37 46 47 51 56 57 ... #$ death: int  1 1 1 1 1 1 1 1 1 1 ...
# $im : int 1 1 1 1 1 1 1 1 1 1 ... table(subset(btrial, death == 1)$time)
# death time is unique
coxph(Surv(time, death) ~ im, data = btrial)
#     coef exp(coef) se(coef)    z     p
# im 0.980     2.665    0.435 2.25 0.024
# Likelihood ratio test=4.45  on 1 df, p=0.03
# n= 45, number of events= 24

# Ties, Section 8.3
data(kidney)
str(kidney)
# 'data.frame':	119 obs. of  3 variables:
# $time : num 1.5 3.5 4.5 4.5 5.5 8.5 8.5 9.5 10.5 11.5 ... #$ delta: int  1 1 1 1 1 1 1 1 1 1 ...
# $type : int 1 1 1 1 1 1 1 1 1 1 ... table(subset(kidney, delta == 1)$time)
# 0.5  1.5  2.5  3.5  4.5  5.5  6.5  8.5  9.5 10.5 11.5 15.5 16.5 18.5 23.5 26.5
# 6    1    2    2    2    1    1    2    1    1    1    2    1    1    1    1

# Default: Efron method
coxph(Surv(time, delta) ~ type, data = kidney)
# coef exp(coef) se(coef)     z    p
# type -0.613     0.542    0.398 -1.54 0.12
# Likelihood ratio test=2.41  on 1 df, p=0.1
# n= 119, number of events= 26
summary(coxph(Surv(time, delta) ~ type, data = kidney))
# n= 119, number of events= 26
# coef exp(coef) se(coef)      z Pr(>|z|)
# type -0.6126    0.5420   0.3979 -1.539    0.124
#
# exp(coef) exp(-coef) lower .95 upper .95
# type     0.542      1.845    0.2485     1.182
#
# Concordance= 0.497  (se = 0.056 )
# Rsquare= 0.02   (max possible= 0.827 )
# Likelihood ratio test= 2.41  on 1 df,   p=0.1
# Wald test            = 2.37  on 1 df,   p=0.1
# Score (logrank) test = 2.44  on 1 df,   p=0.1

# Breslow method
summary(coxph(Surv(time, delta) ~ type, data = kidney, ties = "breslow"))
# n= 119, number of events= 26
#         coef exp(coef) se(coef)      z Pr(>|z|)
# type -0.6182    0.5389   0.3981 -1.553     0.12
#
#       exp(coef) exp(-coef) lower .95 upper .95
# type    0.5389      1.856     0.247     1.176
#
# Concordance= 0.497  (se = 0.056 )
# Rsquare= 0.02   (max possible= 0.827 )
# Likelihood ratio test= 2.45  on 1 df,   p=0.1
# Wald test            = 2.41  on 1 df,   p=0.1
# Score (logrank) test = 2.49  on 1 df,   p=0.1

# Discrete/exact method
summary(coxph(Surv(time, delta) ~ type, data = kidney, ties = "exact"))
#         coef exp(coef) se(coef)      z Pr(>|z|)
# type -0.6294    0.5329   0.4019 -1.566    0.117
#
#      exp(coef) exp(-coef) lower .95 upper .95
# type    0.5329      1.877    0.2424     1.171
#
# Rsquare= 0.021   (max possible= 0.795 )
# Likelihood ratio test= 2.49  on 1 df,   p=0.1
# Wald test            = 2.45  on 1 df,   p=0.1
# Score (logrank) test = 2.53  on 1 df,   p=0.1


## Hazard (function) and survival function

A hazard is the rate at which events happen, so that the probability of an event happening in a short time interval is the length of time multiplied by the hazard.

${\displaystyle h(t)=\lim _{\Delta t\to 0}{\frac {P(t\leq T

Therefore

${\displaystyle H(x)=\int _{0}^{x}h(u)d(u)=-ln[S(x)].}$

or

${\displaystyle S(x)=e^{-H(x)}}$

Hazards (or probability of hazards) may vary with time, while the assumption in proportional hazard models for survival is that the hazard is a constant proportion.

Examples:

• If h(t)=c, S(t) is exponential. f(t) = c exp(-ct). The mean is 1/c.
• If ${\displaystyle \log h(t)=c+\rho t}$, S(t) is Gompertz distribution.
• If ${\displaystyle \log h(t)=c+\rho \log(t)}$, S(t) is Weibull distribution.
• For Cox regression, the survival function can be shown to be ${\displaystyle S(t|X)=S_{0}(t)^{\exp(X\beta )}}$.
{\displaystyle {\begin{aligned}S(t|X)&=e^{-H(t)}=e^{-\int _{0}^{t}h(u|X)du}\\&=e^{-\int _{0}^{t}h_{0}(u)exp(X\beta )du}\\&=e^{-\int _{0}^{t}h_{0}(u)du\cdot exp(X\beta )}\\&=S_{0}(t)^{exp(X\beta )}\end{aligned}}}

Alternatively,

{\displaystyle {\begin{aligned}S(t|X)&=e^{-H(t)}=e^{-\int _{0}^{t}h(u|X)du}\\&=e^{-\int _{0}^{t}h_{0}(u)exp(X\beta )du}\\&=e^{-H_{0}(t)\cdot exp(X\beta )}\end{aligned}}}

where the cumulative baseline hazard at time t, ${\displaystyle H_{0}(t)}$, is commonly estimated through the non-parametric Breslow estimator.

## Sample size calculators

### How many events are required to fit the Cox regression reliably?

If we have only 1 covariate and the covariate is continuous, we need at least 2 events (one for the baseline hazard and one for beta).

If the covariate is discrete, we need at least one event from (each of) two groups in order to fit the Cox regression reliably. For example, if status=(0,0,0,1,0,1) and x=(0,0,1,1,2,2) works fine.

library(survival)
#   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

ova <- ovarian # n=26
ova$time <- ova$futime
ova$status <- 0 ova$status[1:4] <- 1
coxph(Surv(time, status) ~ rx, data = ova) # OK
summary(survfit(Surv(time, status) ~ rx, data =ova))
#                 rx=1
#  time n.risk n.event survival std.err lower 95% CI upper 95% CI
#    59     13       1    0.923  0.0739        0.789            1
#   115     12       1    0.846  0.1001        0.671            1
#   156     11       1    0.769  0.1169        0.571            1
#                 rx=2
#     time  n.risk  n.event  survival  std.err lower 95% CI upper 95% CI
# 421.0000 10.0000   1.0000    0.9000   0.0949       0.7320       1.0000

# Suspicious Cox regression result due to 0 sample size in one group
ova$status <- 0 ova$status[1:3] <- 1
coxph(Surv(time, status) ~ rx, data = ova)
#         coef exp(coef)  se(coef) z p
# rx -2.13e+01  5.67e-10  2.32e+04 0 1
#
# Likelihood ratio test=4.41  on 1 df, p=0.04
# n= 26, number of events= 3
# Warning message:
# In fitter(X, Y, strats, offset, init, control, weights = weights,  :
#   Loglik converged before variable  1 ; beta may be infinite.

summary(survfit(Surv(time, status) ~ rx, data = ova))
#                rx=1
# time n.risk n.event survival std.err lower 95% CI upper 95% CI
#   59     13       1    0.923  0.0739        0.789            1
#  115     12       1    0.846  0.1001        0.671            1
#  156     11       1    0.769  0.1169        0.571            1
#                rx=2
# time n.risk n.event survival std.err lower 95% CI upper 95% CI


## Extract p-values

fit <- coxph(Surv(futime, fustat) ~ age, data = ovarian)

# method 1:
beta <- coef(fit)
se <- sqrt(diag(vcov(fit)))
1 - pchisq((beta/se)^2, 1)

# method 2: https://www.biostars.org/p/65315/
coef(summary(fit))[, "Pr(>|z|)"]


## Expectation of life & expected future lifetime

• The average lifetime is the same as the area under the survival curve.
{\displaystyle {\begin{aligned}\mu &=\int _{0}^{\infty }tf(t)dt\\&=\int _{0}^{\infty }S(t)dt\end{aligned}}}

by integrating by parts making use of the fact that -f(t) is the derivative of S(t), which has limits S(0)=1 and ${\displaystyle S(\infty )=0}$. The average lifetime may not be bounded if you have censored data, there's censored observations that last beyond your last recorded death.

${\displaystyle {\frac {1}{S(t_{0})}}\int _{0}^{\infty }t\,f(t_{0}+t)\,dt={\frac {1}{S(t_{0})}}\int _{t_{0}}^{\infty }S(t)\,dt,}$

## Hazard Ratio (exp^beta) vs Relative Risk

1. https://en.wikipedia.org/wiki/Hazard_ratio
2. Hazard represents the instantaneous event rate, which means the probability that an individual would experience an event (e.g. death/relapse) at a particular given point in time after the intervention, assuming that this individual has survived to that particular point of time without experiencing any event.
3. Hazard ratio is a measure of an effect of an intervention of an outcome of interest over time.
4. Hazard ratio = hazard in the intervention group / Hazard in the control group
5. A hazard ratio is often reported as a “reduction in risk of death or progression” – This risk reduction is calculated as 1 minus the Hazard Ratio (exp^beta), e.g., HR of 0.84 is equal to a 16% reduction in risk. See www.time4epi.com and stackexchange.com.
6. Hazard ratio and its confidence can be obtained in R by using the summary() method; e.g. fit <- coxph(Surv(time, status) ~ x); summary(fit)$conf.int; confint(fit) 7. The coefficient beta represents the expected change in log hazard if X changes by one unit and all other variables are held constant in Cox models. See Variable selection – A review and recommendations for the practicing statistician by Heinze et al 2018. Another example (John Fox, Cox Proportional-Hazards Regression for Survival Data) is assuming Y ~ age + prio + others. • If exp(beta_age) = 0.944. It means an additional year of age reduces the hazard by a factor of .944 on average, or (1-.944)*100 = 5.6 percent. • If exp(beta_prio) = 1.096, it means each prior conviction increases the hazard by a factor of 1.096, or 9.6 percent. Odds Ratio, Hazard Ratio and Relative Risk by Janez Stare For two groups that differ only in treatment condition, the ratio of the hazard functions is given by ${\displaystyle e^{\beta }}$, where ${\displaystyle \beta }$ is the estimate of treatment effect derived from the regression model. See here. Compute ratio ratios from coxph() in R (Hint: exp(beta)). Prognostic index is defined on http://www.math.ucsd.edu/~rxu/math284/slect6.pdf#page=2. Basics of the Cox proportional hazards model. Good prognostic factor (b<0 or HR<1) and bad prognostic factor (b>0 or HR>1). Variable selection: variables were retained in the prediction models if they had a hazard ratio of <0.85 or >1.15 (for binary variables) and were statistically significant at the 0.01 level. see Development and validation of risk prediction equations to estimate survival in patients with colorectal cancer: cohort study. library(KMsurv) # No ties. Section 8.2 data(btrial) coxph(Surv(time, death) ~ im, data = btrial) summary(coxph(Surv(time, death) ~ im, data = btrial))$conf.int
#     exp(coef) exp(-coef) lower .95 upper .95
# im  2.664988  0.3752362  1.136362  6.249912


So the hazard ratio and its 95% ci can be obtained from the 1st, 3rd and 4th element in conf.int.

### Hazard Ratio and death probability

Suppose S0(t)=.2 (20% survived at time t) and the hazard ratio (hr) is 2 (a group has twice the chance of dying than a comparison group), then (Cox model is assumed)

1. S1(t)=S0(t)hr = .22 = .04 (4% survived at t)
2. The corresponding death probabilities are 0.8 and 0.96.
3. If a subject is exposed to twice the risk of a reference subject at every age, then the probability that the subject will be alive at any given age is the square of the probability that the reference subject (covariates = 0) would be alive at the same age. See p10 of this lecture notes.
4. exp(x*beta) is the relative risk associated with covariate value x.

### Hazard Ratio Forest Plot

The forest plot quickly summarizes the hazard ratio data across multiple variables –If the line crosses the 1.0 value, the hazard ratio is not significant and there is no clear advantage for either arm.

## Estimate baseline hazard ${\displaystyle h_{0}(t)}$, Breslow cumulative baseline hazard ${\displaystyle H_{0}(t)}$, baseline survival function ${\displaystyle S_{0}(t)}$ and the survival function ${\displaystyle S(t)}$

Google: how to estimate baseline hazard rate

basehaz(coxph(Surv(time,status)~1,data=aml))


because the (Breslow) hazard estimator for a Cox model reduces to the Nelson-Aalen estimator when there are no covariates. You can also compute it from information returned by survfit().

fit <- survfit(Surv(time, status) ~ 1, data = aml)
cumsum(fit$n.event/fit$n.risk) # the Nelson-Aalen estimator for the times given by fit$times -log(fit$surv)   # cumulative hazard


### Manually compute

Breslow estimator of the baseline cumulative hazard rate reduces to the Nelson-Aalen estimator ${\displaystyle \sum _{t_{i}\leq t}{\frac {d_{i}}{Y_{i}}}}$ (${\displaystyle Y_{i}}$ is the number at risk at time ${\displaystyle t_{i}}$) when there are no covariates present; see p283 of Klein 2003.

{\displaystyle {\begin{aligned}{\hat {H}}_{0}(t)&=\sum _{t_{i}\leq t}{\frac {d_{i}}{W(t_{i};b)}},\\W(t_{i};b)&=\sum _{j\in R(t_{i})}\exp(b'z_{j})\end{aligned}}}

where ${\displaystyle t_{1} denotes the distinct death times and ${\displaystyle d_{i}}$ be the number of deaths at time ${\displaystyle t_{i}}$. The estimator of the baseline survival function ${\displaystyle S_{0}(t)=\exp[-H_{0}(t)]}$ is given by ${\displaystyle {\hat {S}}_{0}(t)=\exp[-{\hat {H}}_{0}(t)]}$. Below we use the formula to compute the cumulative hazard (and survival function) and compare them with the result obtained using R's built-in functions. The following code is a modification of the snippet from the post Breslow cumulative hazard and basehaz().

bhaz <- function(beta, time, status, x) {
# time can be duplicated
# x (covariate) must be continuous
data <- data.frame(time,status,x)
data <- data[order(data$time), ] dt <- unique(data$time)
k    <- length(dt)
risk <- exp(data.matrix(data[,-c(1:2)]) %*% beta)
h    <- rep(0,k)

for(i in 1:k) {
h[i] <- sum(data$status[data$time==dt[i]]) / sum(risk[data$time>=dt[i]]) } return(data.frame(h, dt)) } # Example 1 'ovarian' which has unique survival time all(table(ovarian$futime) == 1) # TRUE

fit <- coxph(Surv(futime, fustat) ~ age, data = ovarian)
# 1. compute the cumulative baseline hazard
# 1.1 manually using Breslow estimator at the observed time
h0 <- bhaz(fit$coef, ovarian$futime, ovarian$fustat, ovarian$age)
H0 <- cumsum(h0$h) # 1.2 use basehaz (always compute at the observed time) # since we consider the baseline, we need to use centered=FALSE H0.bh <- basehaz(fit, centered = FALSE) cbind(H0, h0$dt, H0.bh)
range(abs(H0 - H0.bh$hazard)) # [1] 6.352747e-22 5.421011e-20 # 2. compute the estimation of the survival function # 2.1 manually using Breslow estimator at t = observed time (one dim, sorted) # and observed age (another dim): # S(t) = S0(t) ^ exp(bx) = exp(-H0(t)) ^ exp(bx) S1 <- outer(exp(-H0), exp(coef(fit) * ovarian$age), "^")
dim(S1) # row = times, col = age
# 2.2 How about considering times not at observed (e.g. h0$dt - 10)? # Let's look at the hazard rate newtime <- h0$dt - 10
H0 <- sapply(newtime, function(tt) sum(h0$h[h0$dt <= tt]))
S2 <- outer(exp(-H0),  exp(coef(fit) * ovarian$age), "^") dim(S2) # row = newtime, col = age # 2.3 use summary() and survfit() to obtain the estimation of the survival S3 <- summary(survfit(fit, data.frame(age = ovarian$age)), times = h0$dt)$surv
dim(S3)  # row = times, col = age
range(abs(S1 - S3)) # [1] 2.117244e-17 6.544321e-12
# 2.4 How about considering times not at observed (e.g. h0$dt - 10)? # Note that we cannot put times larger than the observed S4 <- summary(survfit(fit, data.frame(age = ovarian$age)), times = newtime)$surv range(abs(S2 - S4)) # [1] 0.000000e+00 6.544321e-12  # Example 2 'kidney' which has duplicated time fit <- coxph(Surv(time, status) ~ age, data = kidney) # manually compute the breslow cumulative baseline hazard # at the observed time h0 <- with(kidney, bhaz(fit$coef, time, status, age))
H0 <- cumsum(h0$h) # use basehaz (always compute at the observed time) # since we consider the baseline, we need to use centered=FALSE H0.bh <- basehaz(fit, centered = FALSE) head(cbind(H0, h0$dt, H0.bh))
range(abs(H0 - H0.bh$hazard)) # [1] 0.000000000 0.005775414 # manually compute the estimation of the survival function # at t = observed time (one dim, sorted) and observed age (another dim): # S(t) = S0(t) ^ exp(bx) = exp(-H0(t)) ^ exp(bx) S1 <- outer(exp(-H0), exp(coef(fit) * kidney$age), "^")
dim(S1) # row = times, col = age
# How about considering times not at observed (h0$dt - 1)? # Let's look at the hazard rate newtime <- h0$dt - 1
H0 <- sapply(newtime, function(tt) sum(h0$h[h0$dt <= tt]))
S2 <- outer(exp(-H0),  exp(coef(fit) * kidney$age), "^") dim(S2) # row = newtime, col = age # use summary() and survfit() to obtain the estimation of the survival S3 <- summary(survfit(fit, data.frame(age = kidney$age)), times = h0$dt)$surv
dim(S3)  # row = times, col = age
range(abs(S1 - S3)) # [1] 0.000000000 0.002783715
# How about considering times not at observed (h0$dt - 1)? # We cannot put times larger than the observed S4 <- summary(survfit(fit, data.frame(age = kidney$age)), times = newtime)$surv range(abs(S2 - S4)) # [1] 0.000000000 0.002783715  • basehaz() (an alias for survfit) from stackexchange.com and here. basehaz() has a parameter centered which by default is TRUE. Actually basehaz() gives cumulative hazard H(t). See here. That is, exp(-basehaz(fit)$hazard) is the same as summary(survfit(fit))$surv. basehaz() function is not useful. fit <- coxph(Surv(futime, fustat) ~ age, data = ovarian) > fit Call: coxph(formula = Surv(futime, fustat) ~ age, data = ovarian) coef exp(coef) se(coef) z p age 0.1616 1.1754 0.0497 3.25 0.0012 Likelihood ratio test=14.3 on 1 df, p=0.000156 n= 26, number of events= 12 # Note the default 'centered = TRUE' for basehaz() > exp(-basehaz(fit)$hazard) # exp(-cumulative hazard)
[1] 0.9880206 0.9738738 0.9545899 0.9334790 0.8973620 0.8624781 0.8243117
[8] 0.8243117 0.8243117 0.7750981 0.7750981 0.7244924 0.6734146 0.6734146
[15] 0.5962187 0.5204807 0.5204807 0.5204807 0.5204807 0.5204807 0.5204807
[22] 0.5204807 0.5204807 0.5204807 0.5204807 0.5204807
> dim(ovarian)
[1] 26  6
> exp(-basehaz(fit)$hazard)[ovarian$fustat == 1]
[1] 0.9880206 0.9738738 0.9545899 0.8973620 0.8243117 0.8243117 0.7750981
[8] 0.7750981 0.5204807 0.5204807 0.5204807 0.5204807
> summary(survfit(fit))$surv [1] 0.9880206 0.9738738 0.9545899 0.9334790 0.8973620 0.8624781 0.8243117 [8] 0.7750981 0.7244924 0.6734146 0.5962187 0.5204807 > summary(survfit(fit, data.frame(age=mean(ovarian$age))),
time=ovarian$futime[ovarian$fustat == 1])$surv # Same result as above > summary(survfit(fit, data.frame(age=mean(ovarian$age))),
time=ovarian$futime)$surv
[1] 0.9880206 0.9738738 0.9545899 0.9334790 0.8973620 0.8624781 0.8243117
[8] 0.8243117 0.8243117 0.7750981 0.7750981 0.7244924 0.6734146 0.6734146
[15] 0.5962187 0.5204807 0.5204807 0.5204807 0.5204807 0.5204807 0.5204807
[22] 0.5204807 0.5204807 0.5204807 0.5204807 0.5204807


## Predicted survival probability in Cox model: survfit.coxph(), plot.survfit() & summary.survfit( , times)

For theory, see section 8.6 Estimation of the survival function in Klein & Moeschberger.

plot.survfit(). fun="log" to plot log survival curve, fun="event" plot cumulative events, fun="cumhaz" plots cumulative hazard (f(y) = -log(y)).

The plot function below will draw 4 curves: ${\displaystyle S_{0}(t)^{\exp({\hat {\beta }}_{age}*60)}}$, ${\displaystyle S_{0}(t)^{\exp({\hat {\beta }}_{age}*60+{\hat {\beta }}_{stageII})}}$, ${\displaystyle S_{0}(t)^{\exp({\hat {\beta }}_{age}*60+{\hat {\beta }}_{stageIII})}}$,