---
title: "Biostatistics III in R"
author:
- Code by Annika Tillander, Andreas Karlsson and Mark Clements
format:
  html:
    minimal: true
filters:
  - webr
execute:
  echo: true
  message: false
  cache: false
  fig-width: 7
  fig-height: 6
---


# Exercise 2. Comparing survival proportions and mortality rates by stage for cause-specific and all-cause survival #
	
The purpose of this exercise is to study survival of the patients
using two alternative measures --- survival proportions and
mortality rates. A second purpose is to study the difference between
cause-specific and all-cause survival.

You may have to install the required packages the first time you use them. You can install a package by `install.packages("package_of_interest")` for each package you require.

Load dependencies and define two 1/0 variables for the events that we are interested in:

```{webr-r}
#| autorun: true
#| results: "hide"
library(biostat3) # melanoma (depends on survival)
library(bshazard) # bshazard
library(knitr)    # kable
as.data.frame.bshazard <- function(x, ...)
    with(x, data.frame(time,hazard,lower.ci,upper.ci))
melanoma <- 
    transform(biostat3::melanoma,
              death_cancer = ifelse(status == "Dead: cancer", 1, 0),
              death_all = ifelse(status %in% c("Dead: cancer", "Dead: other"), 1, 0))
```

We then list the first few observations to get an idea about the data.

```{webr-r}
#| results: 'asis'
head(biostat3::melanoma) |> knitr::kable("html")
```

## (a) Plot estimates of the survivor function and hazard function by stage ##

We now tabulate the distribution of the melanoma patients by cancer stage at diagnosis.

```{webr-r}
#| results: 'asis'
Freq <- xtabs(~stage, data=melanoma)
cbind(Freq, Prop=proportions(Freq)) |> kable("html")
```

We then plot the survival and hazards by stage. Does it appear that stage is associated with patient survival?

```{webr-r}
library(dplyr)
hazards <- group_by(melanoma, stage) |> 
    do(as.data.frame(bshazard(Surv(surv_mm/12, death_cancer)~1, data=., verbose=FALSE))) |>
    ungroup()
with(hazards, plt(hazard~time|stage, ymin=lower.ci, ymax=upper.ci,
                  type="ribbon",
                  xlab="Time since diagnosis (years)",
                  ylab="Hazard",
                  col=1:4, lty=1, xlim=c(0,20), ylim=0:1))
```


## (b) Estimate the mortality rates for each stage using, for example, the `survRate` command ##

```{webr-r}
#| results: "asis"
survRate(Surv(surv_mm/12, death_cancer) ~ stage, data=melanoma) |> kable("html")
```

What are the units of the estimated rates? The `survRate` function, as the name suggests, is used to estimates rates. Look at the help pages if you are not familiar with the function (e.g. `?survRate` or `help(survRate)`).


We can also do this using more general tools:

```{webr-r}
#| results: "asis"
library(dplyr)
melanoma |>
    group_by(stage) |>
    summarise(D = sum(death_cancer), M = sum(surv_mm/12), Rate = D/M,
              lower.ci = stats::poisson.test(D,M)$conf.int[1],
              upper.ci = stats::poisson.test(D,M)$conf.int[2]) |>
    kable("html")
```

## (c) ##
If you haven’t already done so, estimate the mortality rates for each stage per 1000 person-years of follow-up.

```{webr-r}
#| results: "asis"
survRate(Surv(surv_mm/12/1000, death_cancer) ~ stage, data=melanoma) |> kable("html")
```

## (d) ##

Study whether survival is different for males and females (both by plotting the survivor function and by tabulating mortality rates). Is there a difference in survival between males and females? If yes, is the difference present throughout the follow up?

```{webr-r}
library(dplyr) # slow:(
library(tinyplot) # fast:)
library(ggplot2) # slow:(
hazards <- group_by(melanoma, sex) |>
    do(as.data.frame(bshazard(Surv(surv_mm/12, death_cancer)~1, data=., verbose=FALSE))) |>
    ungroup()
with(hazards,
     plt(hazard~time|sex,ymin=lower.ci, ymax=upper.ci,type="ribbon",
         xlab='Follow-up Time', ylab='Hazard'))
ggplot(hazards,aes(x=time,y=hazard,group=sex)) + geom_line(aes(col=sex)) +
    geom_ribbon(aes(ymin=lower.ci, ymax=upper.ci, fill=sex), alpha=0.3) + ylim(0,0.01) +
    xlab('Follow-up Time') + ylab('Hazard')
```

## (e) ##

The plots you made above were based on cause-specific survival
(i.e., only deaths due to cancer are counted as events, deaths
due to other causes are censored). In the next part of this
question we will estimate all-cause survival (i.e., any death is
counted as an event). First, however, study the coding of vital
status and tabulate vital status by age group.

How many patients die of each cause? Does the distribution of
cause of death depend on age?


```{webr-r}
Freq <- xtabs(~status+agegrp, melanoma)
Freq
# proportions for status by agegrp
round(proportions(Freq,"agegrp")*100,1)
chisq.test(Freq[-4,])
```

## (f) ##

To get all-cause survival, specify all deaths (both cancer
and other) as events.

Now plot the survivor proportion for all-cause survival
by stage. We name the graph to be able to separate them in the
graph window. Is the survivor proportion different compared to
the cause-specific survival you estimated above? Why?

```{webr-r}
par(mfrow=c(1, 1))
afit <- survfit(Surv(surv_mm/12, death_all) ~ stage, data = melanoma)
plot(afit, col=1:4,
     xlab = "Time since diagnosis (years)",
     ylab = "Survival",
     main = "Kaplan-Meier survival estimates\nAll-cause")
legend("topright", levels(melanoma$stage), col=1:4, lty = 1)
```

## (g) ##

It is more common to die from a cause other than cancer in
older ages. How does this impact the survivor proportion for
different stages? Compare cause-specific and all-cause survival
by plotting the survivor proportion by stage for the oldest age
group (75+ years) for both cause-specific and
all-cause survival.

```{webr-r}
par(mfrow=c(1, 2))
mfit75 <- survfit(Surv(surv_mm/12, death_cancer) ~ stage, data = subset(melanoma,agegrp=="75+"))
plot(mfit75, col=1:4,
     xlab = "Since since diagnosis (years)",
     ylab = "Survival",
     main = "Kaplan-Meier survival estimates\nCancer | Age 75+")
legend("topright", levels(melanoma$stage), col=1:4, lty = 1)

afit75 <- survfit(Surv(surv_mm/12, death_all) ~ stage, data = subset(melanoma,agegrp=="75+"))
plot(afit75, col=1:4,
     xlab = "Time since diagnosis (years)",
     ylab = "Survival",
     main = "Kaplan-Meier survival estimates\nAll-cause | Age 75+")
legend("topright", levels(melanoma$stage), col=1:4, lty = 1)
```

## (h) ##

Now estimate both cancer-specific and all-cause survival for
        each age group.


```{r}
par(mfrow=c(1, 2))
survfit(Surv(surv_mm/12, death_cancer) ~ agegrp, data = melanoma) |>
	plot(col=1:4,
         xlab = "Time since diagnosis (years)",
         ylab = "Survival",
         main = "Kaplan-Meier estimates of\ncancer survival by age group")
legend("topright", levels(melanoma$agegrp), col=1:4, lty = 1)

survfit(Surv(surv_mm/12, death_all) ~ agegrp, data = melanoma) |>
    plot(col=1:4,
         xlab = "Time since diagnosis (years)",
         ylab = "Survival",
         main = "Kaplan-Meier estimates of\nall-cause survival by age group")
legend("topright", levels(melanoma$agegrp), col=1:4, lty = 1)
```
