---
title: "Smooth modeling of right-censored survival data with {icpack}"
author: "Hein Putter & Paul Eilers"
date: "`r format(Sys.time(), '%d %B, %Y')`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Smooth modeling of right-censored survival data with {icpack}}
  %\VignetteEngine{knitr::knitr}
  %\VignetteEncoding[UTF-8]{inputenc}
---


```{r global_options, include=FALSE}
knitr::opts_chunk$set(fig.width = 6, fig.height = 4,
                      echo = TRUE, warning = FALSE, message = FALSE)
```

# Introduction

The purpose of this vignette is to illustrate the use of {`icpack`}, using the right-censored `Ova` data, and to contrast the smooth modeling using {`icpack`} with the more standard non- and semi-parametric approaches of the `survival` package.

# Ovarian cancer data

The ovarian cancer data are available after loading the package.

```{r ova}
library(icpack)
```

## Preliminaries

There are three covariates (see the help file). Here is a summary of their frequencies and distribution.

```{r dusersummary}
table(Ova$Diameter)
table(Ova$FIGO)
table(Ova$Karnofsky)
```

`Diameter` is the diameter of the tumor, `FIGO` is the FIGO classification, where in fact `FIGO=0` corresponds to a FIGO III and `FIGO=1` to a FIGO IV classification, and `Karnofsky` is the often used Karnofsky score divided by 10, so for instance `Karnofsky=9` corresponds to a Karnofksy score of 90.

## Using icfit

The time-to-event data are right-censored. Both `Diameter` and `Karnofsky` are ordinal, so we could use them both as categorical and as continuous variables. Let's first use them both as categorical, and see whether it makes sense to use them as continuous variables by inspecting their (log) hazard ratios.

```{r icf1}
ic_fit <- icfit(Surv(time, d) ~ factor(Diameter) + FIGO + factor(Karnofsky),
                data = Ova)
ic_fit
```

It looks like both for `Diameter` and `Karnofsky` the increments of the log hazard ratios are reasonably equal for neighboring factor levels, which could argue for using both as continuous covariates in the future. Before we do that, however, let's compare the result with that from `coxph` in the `survival` package.

```{r coxph1}
cox_fit <- coxph(Surv(time, d) ~ factor(Diameter) + FIGO + factor(Karnofsky), data = Ova)
cox_fit
```

The results, in terms of the hazard ratios and their standard errors, are very similar indeed.

Let's refit the smooth and semi-parametric model, now using `Diameter` and `Karnofsky` as continuous variables.

```{r icfcoxph}
ic_fit <- icfit(Surv(time, d) ~ Diameter + FIGO + Karnofsky, data = Ova)
cox_fit <- coxph(Surv(time, d) ~ Diameter + FIGO + Karnofsky, data = Ova)
ic_fit
cox_fit
```

Without claiming in any way that using `Diameter` and `Karnofsky` as continuous variables is better, we continue with these last two models, mainly for simplicity.

## Plotting

The main difference between the smooth and semi-parametric fits of course resides in the estimates of the baseline hazard. We illustrate this by extracting the baseline hazards and plotting the corresponding baseline survival curves.

The `icpack` package uses `ggplot2` for plotting. The default plot is the (baseline) hazard; by setting `type="survival"` we get the (baseline) survival function. Default option is to add shaded pointwise confidence intervals; they have been switched off here.

```{r ploticf}
plot(ic_fit, type="survival", conf.int = FALSE,
     title = "Baseline survival function from icfit")
```

Here is the corresponding baseline survival function from `coxph`, with the smooth baseline survival curve from `icfit` added to it in blue.

```{r baselinecoxph}
bh <- basehaz(cox_fit, centered=FALSE)
bh$surv <- exp(-bh$hazard)
plot(c(0, bh$time), c(1, bh$surv), type="s",
     xlab="Time", ylab="Survival", main="Baseline survival functions from coxph and icfit")
nd_baseline <- data.frame(Diameter = 0, FIGO = 0, Karnofsky = 0)
pred_baseline <- predict(ic_fit, newdata = nd_baseline)[[1]]
lines(pred_baseline$time, pred_baseline$Surv, col = "blue")
legend("topright", c("coxph", "icfit"), lwd=1, col=c("black", "blue"), bty="n")
```

The plot does not show pointwise confidence intervals; it is possible to obtain those using `predict`, which is illustrated in the next section. Obviously, the estimated baseline survival curve form `coxph` is less smooth; their values however are very similar.

The estimate of the baseline hazard from the Cox model fit is not smooth; plotting it does not produce anything remotely useful.

```{r coxbasehaz}
plot(bh$time, diff(c(0, bh$hazard)),
     xlab="Time", ylab="Hazard", main="Baseline hazard function from coxph")
```

The baseline hazard obtained from `icfit` is smooth and does give a useful illustration of the hazard, along with its pointwise confidence intervals.

```{r ploticfit2}
plot(ic_fit, type = 'hazard', title = 'Baseline hazard function from icfit')
```

The cumulative hazard can also be shown using `type="cumhazard"`.

## Predicting

Actually the baseline hazards and survival function are not really that interesting, in the sense that they correspond to a subject with `Diameter=0`, `FIGO=0` and `Karnofsky=0`. Such persons (at least `Karnofsky=0`) are not in the data and are even biologically impossible.

Modeled after the `survival` package, the function `predict` with a `newdata` object can be used to obtain hazards and survival functions for subjects with given characteristics. Let us illustrate this with subject A, who has `Diameter=0`, `FIGO=0` and `Karnofsky=10`, the most optimistic possibility, and  subject B, who has `Diameter=4`, `FIGO=0` and `Karnofsky=6`, the most pessimistic possibility, apart from FIGO. This is simply entered into a data frame with the same names as the original data.

```{r newd}
newd <- data.frame(subject = c("A", "B"),
                   Diameter = c(0, 4),
                   FIGO = c(0, 0),
                   Karnofsky = c(10, 6))
newd
```

```{r predict}
picf <- predict(ic_fit, newdata=newd)
```

The result is an object of class `predict.icfit`, which is a list with 2 items (one for each subject), each containing a data frame with time points and hazards, cumulative hazards, survival functions at these time points with associated standard errors and confidence bounds. The `summary` method can be used to extract the hazards etc at specified time points.

```{r summarypredicticfit}
summary(picf, times=(0:7) * 365.25)
```

Finally, the hazards or cumulative hazards or survival functions can be plotted.

```{r icfitAB}
library(gridExtra)
plot(picf, type="surv", ylim=c(0, 1), xlab="Days since randomisation",
     title = c("Subject A", "Subject B"), nrow=1)
```

The results from `coxph` for these patients are shown in the plots below as black lines, in blue the results from `icpack`.

```{r coxAB}
sfAB <- survfit(cox_fit, newdata=newd)
sfA <- sfAB[1]
sfB <- sfAB[2]
# Retrieve separate plots from icfit for comparison
picf1 <- picf[[1]]
class(picf1) <- c("predict.icfit", "data.frame")
plt1 <- plot(picf1, type='surv', ylim=c(0, 1),
             xlab = "Days since randomisation",
             title = "Subject A")
picf2 <- picf[[2]]
class(picf2) <- c("predict.icfit", "data.frame")
plt2 <- plot(picf2, type='surv', ylim=c(0, 1),
             xlab = "Days since randomisation",
             title = "Subject B")

oldpar <- par(mfrow=c(1, 2))
plot(sfA, xlab="Days since randomisation", ylab="Survival", main="Subject A")
lines(plt1$data$time, plt1$data$Surv, col = "blue")
lines(plt1$data$time, plt1$data$lowerSurv, col = "blue")
lines(plt1$data$time, plt1$data$upperSurv, col = "blue")
legend("topright", c("coxph", "icfit"), lwd=1, col=c("black", "blue"), bty="n")
plot(sfB, xlab="Days since randomisation", ylab="Survival", main="Subject B")
lines(plt2$data$time, plt2$data$Surv, col = "blue")
lines(plt2$data$time, plt2$data$lowerSurv, col = "blue")
lines(plt2$data$time, plt2$data$upperSurv, col = "blue")
legend("topright", c("coxph", "icfit"), lwd=1, col=c("black", "blue"), bty="n")

par(oldpar)
```

Again, the estimates and their standard errors are very similar.