---

title: "Penalized logistic regression on time-to-event data using casebase sampling"

author: "Jesse Islam"

date: "2/27/2020"

mathspec: true

theme: "metropolis"

output:

  beamer_presentation:

    incremental: FALSE

    fig_caption: false

header-includes:

  \usepackage{xcolor}

---

```{r setup, include=FALSE}

source("packages.R")

set.seed(1)

knitr::opts_chunk$set( warning=FALSE,message=FALSE,echo=FALSE)

knitr::opts_chunk$set(dev = 'pdf')

```

# Popular methods in time-to-event analysis 

* In disease etiology, we tend to make use of the proportional hazards hypothesis.

    * Cox Regression

* When we want the absolute risk:

    * Breslow estimator

    * Parametric models

# Motivations for a new method

* Julien and Hanley found that survival analysis rarely produces prognostic functions, even though the software is widely available in cox regression packages. [2]

* They believe the stepwise nature is the reason, as it reduces interpretability. [2]

* A streamlined approach for reaching a **smooth absolute risk** curve. [2]

# Dr. Cox's perspective

![](coxquote.png)

# Itinerary

1. SUPPORT study

2. Casebase sampling

3. Penalized logistic regression on survival data

4. Maximum likelihood with regularization

5. Absolute risk comparison

6. Conclusion

7. Future work

8. References

# SUPPORT dataset [4]

* **Study to Understand Prognoses and Preferences for Outcomes and Risks Treatments**

* Design: Prospective cohort study.

* Setting: 5 academic care centers in the United States.

* Participants: 9105 hospitalized.

* Follow-up-time: 5.56 years.

* 68% incidence rate.

# SUPPORT manual imputation [4]

* Notorious for missing data

\begin{table}[]

\centering

\resizebox{\textwidth}{!}{%

\begin{tabular}{ll}

\hline

Baseline Variable      & Normal Fill-in Value \\ \hline

Bilirubin              & 1.01                 \\

BUN                    & 6.51                 \\

Creatinine             & 1.01                 \\

PaO2/FiO2 ratio (pafi) & 333.3                \\

Serum albumin          & 3.5                  \\

Urine output           & 2502                 \\

White blood count      & 9 (thousands)        \\ \hline

\end{tabular}%

}

\caption{Suggested imputation values. {[}3{]}}

\label{tab:my-table}

\end{table}

# SUPPORT automated imputation

* Mice imputation package (R)

1. PMM (Predictive Mean Matching)  – For numeric variables

2. logreg(Logistic Regression) – Binary Variables 

3. polyreg(Bayesian polytomous regression) Factor Variables

# Removed variables [4]

* **Response variables**

  * Hospital Charges.

  * Patient ratio of costs to charges.

  * Patient Micro-costs.

* **Ordinal covariates**

  * functional disability.

  * Income.

* **Sparse covariate**

  * Surrogate activities of daily living.

* Previous model results. (6)

# Variable overview [4]

* **Response variables**

  * follow-up time, death.

* **Covariates**

  * Age, sex, race, education (6)

  * Disease group/class, comorbidities. (3)

  * Coma score, Therapeutic Intervention Scoring System (2)

  * Physiological variables. (11)

  * Activities of daily living. (2)

# Original SUPPORT analysis [4]

* Determined SUPPORT prognostic model on phase I patients.

* Tested on Phase II.

* Both on the scale of 180 days.

# Original SUPPORT analysis [4]

SUPPORT physiology score (SPS) was developed.

```{r, out.width = "80%",include=TRUE}

knitr::include_graphics("SPS.png")

```

# SUPPORT: my question

* How does their SPS perform over 5.56 years?

* How does the Apache III physiology score (APS) perform over 5.56 years?

* How does a model with all the covariates, excluding SPS and APS, perform?

# Analysis Process

1. Impute

2. Compare SPS and APS over ~5.56 years using absolute risk.

3. Compare to Kaplan-Meier curve (unadjusted model).

4. Compare to full model (excluding SPS and APS).

* All models is trained on 80% of the observations.

* Remaining observations are used to generate comparative absolute risk curves.

  * The absolute risk curve for each individual is averaged.

  * Averaged curve is expected to approach Kaplan-meier curve.

```{r data processing,include=FALSE}

#DATA retrieved from http://biostat.mc.vanderbilt.edu/wiki/pub/Main/DataSets/support2csv.zip

ew <- read_csv("support2.csv")

data<-read.csv('support2.csv', col.names=colnames(ew), fill=TRUE, header=TRUE)

rm(ew)

```

```{r cleaningImputation,include=FALSE}

#manual imputation

data$alb[is.na(data$alb)]=3.5

data$pafi[is.na(data$pafi)]=333.3

data$bili[is.na(data$bili)]=1.01

data$crea[is.na(data$crea)]=1.01

data$bun[is.na(data$bun)]=6.51

data$wblc[is.na(data$wblc)]=9

data$urine[is.na(data$urine)]=2502

# automated imputation

previousModelPositions<-c(18,19,20,21,22,23,24,25,26,27,28,29)

supportMain<-data[,-previousModelPositions]

supportPrevious<-data[,previousModelPositions]

#keep 1,2 after imputation from previous usable models.

#missingPattern<-md.pattern(supportMain)

supportMainImpute <-mice(supportMain[,-c(11,34,13,14,15)], m = 1,printFlag=FALSE) #impute while removing ordinal variables and other response variables.

imputedData<-cbind(complete(supportMainImpute,1),supportPrevious[,c(1,2)])

#adls is missing 33%, and is removed accordingly and hospitaldeath

#md.pattern(completeData)

completeData<-na.omit(imputedData[,-c(4,29)])

workingCompleteData=model.matrix(death~ .-d.time,data=completeData)[,-c(1)]#remove intercept

#Create u and xCox, which will be used when fitting Cox with glmnet.

newData=workingCompleteData

x=workingCompleteData

#x and y will be used to fit the casebase model under glmnet.

#x=as.matrix(sparse.model.matrix(death~ .-d.time+0,data=workingData))

y=data.matrix(completeData[,c(2,5)])

```

```{r CoxExplorationModels,include=FALSE, eval = TRUE}

fullDataCox<-as.data.frame(cbind(as.numeric(y[,2]),as.numeric(y[,1]),x))

train_index <- sample(1:nrow(fullDataCox), 0.8 * nrow(fullDataCox))

test_index <- setdiff(1:nrow(fullDataCox), train_index)

train<-fullDataCox[train_index,]

test<-fullDataCox[test_index,]

coxSPS<- survival::coxph(Surv(time=V1,event=V2) ~ sps, data = train)

abcoxSPS<-survival::survfit(coxSPS,newdata = test)

coxAPS<- survival::coxph(Surv(time=V1,event=V2) ~ aps, data = train)

abcoxAPS<-survival::survfit(coxAPS,newdata = test)

coxFull<- survival::coxph(Surv(time=V1,event=V2) ~ .-sps -aps, data = train)

abcoxFull<-survival::survfit(coxAPS,newdata = test)

KM<- survival::coxph(Surv(time=V1,event=V2) ~ 1, data = test)

abKM<-survival::survfit(KM)

abKMcompare<-cbind(abKM$time[abKM$time<=1631],1-abKM$surv[abKM$time<=1631])

colnames(abKMcompare)<-c("Time","KM")

baseCoxComparisonData<-cbind(abcoxSPS$time[abcoxSPS$time<=1631],1-rowMeans(abcoxSPS$surv[abcoxSPS$time<=1631,]),1-rowMeans(abcoxAPS$surv[abcoxSPS$time<=1631,]),1-abKM$surv,1-rowMeans(abcoxFull$surv[abcoxSPS$time<=1631,]))

colnames(baseCoxComparisonData)<-c("Time","SPS","APS","Unadjusted","Full")

myColors <- brewer.pal(5,"Dark2")

FigureWidth='90%'

```

# SPS vs APS

```{r SPSvsAPS,include=TRUE, out.width = FigureWidth}

#Full

ggplot(as.data.frame(baseCoxComparisonData), aes(Time)) + 

  geom_line(aes(y = SPS, colour = "SPS"),lwd=2) +

  geom_line(aes(y = APS, colour = "APS"),lwd=2) +

  scale_colour_manual(name = "Models",values = myColors[c(3,4)])+

  labs(y= "Probability of Death", x = "Survival time (Days)")

```

# SPS vs. Kaplan-Meier

```{r SPSvsKM,include=TRUE, out.width = FigureWidth}

ggplot(as.data.frame(baseCoxComparisonData), aes(Time)) + 

  geom_line(aes(y = SPS, colour = "SPS"),lwd=2) +

  geom_line(aes(y = APS, colour = "APS"),lwd=2) +

  geom_step(data=as.data.frame(abKMcompare),mapping=aes(x=Time,y=KM,color="Unadjusted"),lwd=2)+

  scale_colour_manual(name = "Models",values = myColors[c(3,4,2)])+

  labs(y= "Probability of Death", x = "Survival time (Days)")

```

# All covariates vs. physiology scores vs unadjusted

```{r FullvsKM,include=TRUE, out.width = FigureWidth}

ggplot(as.data.frame(baseCoxComparisonData), aes(Time)) + 

  geom_line(aes(y = SPS, colour = "SPS"),lwd=2) + 

  geom_line(aes(y = APS, colour = "APS"),lwd=2) +

  geom_step(data=as.data.frame(abKMcompare),mapping=aes(x=Time,y=KM,color="Unadjusted"),lwd=2)+

  geom_line(aes(y = Full, col = "All Covariates without SPS APS"),lwd=2)+

  scale_colour_manual(name = "Models",values = myColors[c(1,3,4,2)])+

  labs(y= "Probability of Death", x = "Survival time (Days)")

```

# Chosen absolute risk comparisons

```{r ,include=TRUE, out.width = FigureWidth}

ggplot(as.data.frame(baseCoxComparisonData), aes(Time)) + 

  geom_step(data=as.data.frame(abKMcompare),mapping=aes(x=Time,y=KM,color="Unadjusted"),lwd=2)+

  geom_line(aes(y = Full, colour = "Full"),lwd=2)+

  scale_colour_manual(name = "Models",values = myColors[c(1,2)])+

  labs(y= "Probability of Death", x = "Survival time (Days)")

```

# Chosen absolute risk comparisons: conclusion

* Linear associations without physiology scores perform similarly to SPS and APS alone.

* We choose the linear associations without physiology scores as the model of choice (Full model). 

# Casebase sampling overview

1. Clever sampling.

2. Implicitly deals with censoring.

3. Allows a parametric fit using *logistic regression*.

* Casebase is parametric, and allows different parametric fits by incorporation of the time component.

* Package contains an implementation for generating *population-time* plots.

# Casebase: Sampling

```{r DataModel, echo=FALSE, warning=FALSE,message=FALSE,include=TRUE, out.width = FigureWidth}

nobs <- nrow(data)

ftime <- data$d.time

ord <- order(ftime, decreasing = FALSE)

#ordSample=order(sample(ord,1000))

yCoords <- cbind(cumsum(data[ord, "death"] == 1),

cumsum(data[ord, "death"] == 0))

yCoords <- cbind(yCoords, nobs - rowSums(yCoords))

aspectRatio <- 0.75

height <- 8.5 * aspectRatio; width <- 11 * aspectRatio

cases <- data$death == 1

comps <-data$death == 2

# randomly move the cases vertically

moved_cases <- yCoords[cases[ord], 3] * runif(sum(cases))

moved_comps <- yCoords[comps[ord], 3] * runif(sum(comps))

moved_bases <- yCoords[cases[ord], 3] * runif(sum(cases))

plot(0, type = 'n', xlim = c(0, max(ftime)), ylim = c(0, nobs),

xlab = 'Follow-up time (days)', ylab = 'Population',cex.lab=1.5, cex.axis=1.5, cex.main=1.5, cex.sub=1.5)

ptsz=0.5

```

# Casebase: Sampling

```{r , echo=FALSE, warning=FALSE,message=FALSE,include=TRUE, out.width = FigureWidth}

plot(0, type = 'n', xlim = c(0, max(ftime)), ylim = c(0, nobs),

xlab = 'Follow-up time (days)', ylab = 'Population',cex.lab=1.5, cex.axis=1.5, cex.main=1.5, cex.sub=1.5)

polygon(c(0, 0, ftime[ord], max(ftime), 0),

c(0, nobs, yCoords[,3], 0, 0),

col = "grey60")

legend("topright", legend = c("Base"),

col = c("grey60"),

pch = 19,cex=1.5)

```

# Casebase: Sampling

```{r , echo=FALSE, warning=FALSE,message=FALSE,include=TRUE, out.width = FigureWidth}

plot(0, type = 'n', xlim = c(0, max(ftime)), ylim = c(0, nobs),

xlab = 'Follow-up time (days)', ylab = 'Population',cex.lab=1.5, cex.axis=1.5, cex.main=1.5, cex.sub=1.5)

polygon(c(0, 0, ftime[ord], max(ftime), 0),

c(0, nobs, yCoords[,3], 0, 0),

col = "grey60")

points((ftime[ord])[cases[ord]], moved_bases, pch = 19,

col = "blue", cex = ptsz)

legend("topright", legend = c("Base","base-series"),

col = c("grey60","blue"),

pch = 19,cex=1.5)

```

# Casebase: Sampling

```{r , echo=FALSE, warning=FALSE,message=FALSE,include=TRUE, out.width = FigureWidth}

plot(0, type = 'n', xlim = c(0, max(ftime)), ylim = c(0, nobs),

xlab = 'Follow-up time (days)', ylab = 'Population',cex.lab=1.5, cex.axis=1.5, cex.main=1.5, cex.sub=1.5)

polygon(c(0, 0, ftime[ord], max(ftime), 0),

c(0, nobs, yCoords[,3], 0, 0),

col = "grey60")

points((ftime[ord])[cases[ord]], moved_bases, pch = 19,

col = "blue", cex = ptsz)

points((ftime[ord])[cases[ord]], yCoords[cases[ord],3], pch = 19,

col = "firebrick3", cex = ptsz)

legend("topright", legend = c("Base","base-series","case-series (Death)"),

col = c("grey60","blue","firebrick3"),

pch = 19,cex=1.5)

```

# Casebase: Sampling

```{r , echo=FALSE, warning=FALSE,message=FALSE,include=TRUE, out.width = FigureWidth}

plot(0, type = 'n', xlim = c(0, max(ftime)), ylim = c(0, nobs),

xlab = 'Follow-up time (days)', ylab = 'Population',cex.lab=1.5, cex.axis=1.5, cex.main=1.5, cex.sub=1.5)

polygon(c(0, 0, ftime[ord], max(ftime), 0),

c(0, nobs, yCoords[,3], 0, 0),

col = "grey60")

points((ftime[ord])[cases[ord]], moved_bases, pch = 19,

col = "blue", cex = ptsz)

points((ftime[ord])[cases[ord]], moved_cases, pch = 19,

col = "firebrick3", cex = ptsz)

#points((ftime[ord])[bases[ord]], moved_bases, pch = 19,

#col = "blue", cex = ptsz)

legend("topright", legend = c("Base","base-series","case-series (Death)"),

col = c("grey60", "blue","firebrick3"),

pch = 19,cex=1.5)

```

# Casebase: Sampling

$$e^{L}=\frac{Pr(Y=1|x,t)}{Pr(Y=0|x,t)}=\frac{h(x,t)*B(x,t)}{b[B(x,t)/B]}=\frac{h(x,t)*B}{b}$$

* $L=\beta X$

* $b$ = base-series.

* $B$ = Base.

* *B(x,t)* =  Risk-set for survival time $t$.

# Casebase: Sampling

$e^{L}=\frac{Pr(Y=1|x,t)}{Pr(Y=0|x,t)}=\frac{h(x,t)*B(x,t)}{b[B(x,t)/B]}=\frac{h(x,t)*B}{b}$

```{r , echo=FALSE, warning=FALSE,message=FALSE,include=TRUE, out.width = FigureWidth}

plot(0, type = 'n', xlim = c(0, max(ftime)), ylim = c(0, nobs),

xlab = 'Follow-up time (days)', ylab = 'Population',cex.lab=1.5, cex.axis=1.5, cex.main=1.5, cex.sub=1.5)

polygon(c(0, 0, ftime[ord], max(ftime), 0),

c(0, nobs, yCoords[,3], 0, 0),

col = "grey60")

points((ftime[ord])[cases[ord]], moved_bases, pch = 19,

col = "blue", cex = ptsz)

points((ftime[ord])[cases[ord]], moved_cases, pch = 19,

col = "firebrick3", cex = ptsz)

#points((ftime[ord])[bases[ord]], moved_bases, pch = 19,

#col = "blue", cex = ptsz)

abline(v=990,lwd=2)

abline(v=1010,lwd=2)

legend("topright", legend = c("Base","base-series","case-series (Death)"),

col = c("grey60", "blue","firebrick3"),

pch = 19,cex=1.5)

```

# log-odds = log hazard

 $$e^{L}=\frac{\hat{h}(x,t)*B}{b} $$

 $$ \hat{h}(x,t)= \frac{b*e^{L}}{B}$$

 $$log(\hat{h}(x,t))= L+log(\frac{b}{B})$$

# Maximum log-likelihood [1]

 $$  log(l(\beta_{0},\beta))=\frac{1}{N}\sum^{N}_{i=1}\{y_{i}(\beta_{0}+x_{i}^{T}\beta)-log(1+e^{\beta_{0}+x_{i}^{T}\beta})\}  $$

# Maximum log-likelihood, with offset

 $$  log(l(\textcolor{blue}{log(\frac{b}{B})},\beta))=\frac{1}{N}\sum^{N}_{i=1}\{y_{i}(\textcolor{blue}{log(\frac{b}{B})}+x_{i}^{T}\beta)-log(1+e^{\textcolor{blue}{log(\frac{b}{B})}+x_{i}^{T}\beta})\}  $$

# Maximum log-likelihood, with offset and lasso

$\frac{1}{N}\sum^{N}_{i=1}\{y_{i}(\textcolor{blue}{log(\frac{b}{B})}+x_{i}^{T}\beta)-log(1+e^{\textcolor{blue}{log(\frac{b}{B})}+x_{i}^{T}\beta})\}  \textcolor{red}{-\lambda||\beta||}$

# Casebase: Parametric families

* We can now fit models of the form:

$$\ log(h(t;\alpha,\beta))=g(t;\alpha)+\beta X$$

* By changing the function $g(t;\alpha)$, we can model different parametric families easily:

# Casebase: Parametric models

*Exponential*: $g(t;\alpha)$ is equal to a constant

```{r,echo=TRUE,eval=FALSE}

casebase::fitSmoothHazard(status ~ X1 + X2)

```

*Gompertz*: $g(t;\alpha)=\alpha t$

```{r,echo=TRUE,eval=FALSE}

casebase::fitSmoothHazard(status ~ time + X1 + X2)

```

*Weibull*: $g(t;\alpha) = \alpha log(t)$

```{r,echo=TRUE,eval=FALSE}

casebase::fitSmoothHazard(status ~ log(time) + X1 + X2)

```

# Absolute Risk

* We have a bunch of different parametric hazard models now.

* To get the absolute risk, we need to evaluate the following equation in relation to the hazard:

$$CI(x,t)=1-e^{-\int^{t}_{0}h(x,u)du}$$

* *CI(x,t)*= Cumulative Incidence (Absolute Risk)

* *h(x,u)*= Hazard function

* Lets use the weibull hazard.

# Models to be compared

* Recall: Case study to demonstrate regularization using our method.

* **unadjusted**: (death, time) ~ 1

* **Cox**: (death, time) ~  $\beta X$

* **Coxnet**: (death, time) ~  $\beta X$ <--Lasso

* **Penalized logistic**: death ~ log(time) + $\beta X$ <--Lasso 

```{r supportSetUp, echo = FALSE, eval = TRUE}

ratio=10  #universal ratio for TRUE

#keep one individual out, at random, for which we will develop an 

#absolute risk curve.

#sam=sample(1:nrow(completeData), )

#Set workingData to the remaining individuals in the data

#workingData=completeData[-c(sam),]

#x and y will be used to fit the casebase model under glmnet.

#x=as.matrix(sparse.model.matrix(death~ .-d.time+0,data=workingData))

y=data.matrix(train[,c(2,1)])

colnames(y)<-c("death","d.time")

x=data.matrix(train[,-c(2,1)])

newData=data.matrix(test[,-c(2,1)])

```

```{r coxHazAbsolute, echo = FALSE }

#Create u and xCox, which will be used when fitting Cox with glmnet.

u=survival::Surv(time = as.numeric(y[,2]), event = as.factor(y[,1]))

xCox=as.matrix(cbind(data.frame("(Intercept)"=rep(1,each=length(x[,1]))),x))

colnames(xCox)[2]="(Intercept)"

#hazard for cox using glmnet

coxNet=glmnet::cv.glmnet(x=x,y=u, family="cox",alpha=1,standardize = TRUE)

#convergence demonstrated in plot

#plot(coxNet)

#taking the coefficient estimates for later use

nonzero_covariate_cox <- predict(coxNet, type = "nonzero", s = "lambda.1se")

nonzero_coef_cox <- coef(coxNet, s = "lambda.1se")

#creating a new dataset that only contains the covariates chosen through glmnet.

#cleanCoxData<- as.data.frame(cbind(as.numeric(workingData$d.time),as.factor(workingData$death), xCox[,nonzero_covariate_cox$X1]))

cleanCoxData<-as.data.frame(cbind(d.time=as.numeric(y[,2]),death=as.numeric(y[,1]),x[,nonzero_covariate_cox$X1]))

#newDataCox<-xCox[sam,nonzero_covariate_cox$X1]

#fitting a cox model using regular estimation, however we will not keep it.

#this is used more as an object place holder.

coxNet <- survival::coxph(Surv(time=d.time,event=death) ~ ., data = cleanCoxData)

#The coefficients of this object will be replaced with the estimates from coxNet.

#Doing so makes it so that everything is invalid aside from the coefficients.

#In this case, all we need to estimate the absolute risk is the coefficients.

#Std. error would be incorrect here, if we were to draw error bars.

coxNet$coefficients<-nonzero_coef_cox@x

#Fitting absolute risk curve for cox+glmnet. -1 to remove intercept

#posCox=nonzero_covariate_cox$X1

#newDataCoxI=cleanCoxDataI[,posCox]

abcoxNet<-survival::survfit(coxNet,type="breslow",newdata=as.data.frame(newData[,nonzero_covariate_cox$X1]))

abCN=as.data.frame(cbind(abcoxNet$time, 1-rowMeans(abcoxNet$surv)))

colnames(abCN)<-c("Time","Absolute risk")

abCN$Model="Coxnet"

abCN=abCN[abCN$Time<=max(abKMcompare[,1]),]

#plot(abCN[,1],1-abCN[,2])

```

```{r cbHazard}, echo = FALSE, eval = TRUE}

cbFull=casebase::fitSmoothHazard.fit(x,y,family="glmnet",time="d.time",event="death",formula_time = ~log(d.time),alpha=1,ratio=ratio,standardize = TRUE)

```

```{r CBabsolute, echo = FALSE, eval = TRUE}

#Estimating the absolute risk curve using the newData parameter.

abCbFull=casebase::absoluteRisk(cbFull,time = seq(1,max(abKMcompare[,1]), 25),newdata = newData , s="lambda.1se",method=c("numerical"))

abCB=data.frame(abCbFull[,1],rowMeans(abCbFull[,-c(1)]))

colnames(abCB)<-c("Time","Absolute risk")

abCB$Model="Penalized logistic"

ab=rbind(abCB,abCN)

```

# Survival comparison

```{r include=TRUE}

ggplot(as.data.frame(ab), aes(Time)) + 

    geom_line(data=as.data.frame(baseCoxComparisonData),mapping=aes(y = Full, colour = "Cox"),lwd=2)+

  geom_line(aes(y = `Absolute risk`, colour = Model),lwd=2)+

  geom_step(data=as.data.frame(abKMcompare),mapping=aes(x=Time,y=KM,color="Unadjusted"),lwd=2)+

  scale_colour_manual(name = "Models",values = myColors[c(4,1,3,2)])+

  labs(y= "Probability of Death", x = "Survival time (Days)")

```

# Conclusions / take homes

* Classical survival analysis requires methods to encorporate censorship in our data.

* Case-base sampling is a techinique that implicitly encorporates censorship implicitly.

* Logistic regression on SUPPORT dataset had slightly different results near the end of follow-up time. 

# Future work

* Comparative measure.

* Survival GWAS.

# References 1

1.Czepiel, S. A. (2002). Maximum likelihood estimation of logistic regression models: theory and implementation. Available at czep. net/stat/mlelr. pdf, 1825252548-1564645290.

2.Hanley, James A, and Olli S Miettinen. 2009. "Fitting Smooth-in-Time Prognostic Risk Functions via Logistic Regression." *The International Journal of Biostatistics 5 (1)*.

3.Harrell, F. (2020). SupportDesc < Main < Vanderbilt Biostatistics Wiki. [online] Biostat.mc.vanderbilt.edu. Available at: http://biostat.mc.vanderbilt.edu/wiki/Main/SupportDesc [Accessed 25 Feb. 2020].

# References 2

4.Knaus, W. A., Harrell, F. E., Lynn, J., Goldman, L., Phillips, R. S., Connors, A. F., ... & Layde, P. (1995). The SUPPORT prognostic model: Objective estimates of survival for seriously ill hospitalized adults. Annals of internal medicine, 122(3), 191-203.

5.Saarela, Olli, and Elja Arjas. 2015. "Non-Parametric Bayesian Hazard Regression for Chronic Disease Risk Assessment." Scandinavian Journal of Statistics 42 (2). Wiley Online Library: 609–26. 

6.Saarela, Olli. 2015. "A Case-Base Sampling Method for Estimating Recurrent Event Intensities." *Lifetime Data Analysis*. Springer, 1–17

# References 3

7.Scrucca L, Santucci A, Aversa F. Competing risk analysis using R: an easy guide for clinicians. *Bone Marrow Transplant*. 2007 Aug;40(4):381-7. doi: 10.1038/sj.bmt.1705727.

8.Turgeon, M. (2017, June 10). Retrieved May 05, 2019, from https://www.maxturgeon.ca/slides/MTurgeon-2017-Student-Conference.pdf

# Tutorial and slides

<center>

**Tutorial:**

http://sahirbhatnagar.com/casebase/

**Slides:**

https://github.com/Jesse-Islam/ATGC_survival_presentation_Feb.27.2020

**Questions?**

</center>