---
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 10. Examining the proportional hazards hypothesis (localised melanoma) #

-----------

Load the diet data using time-on-study as the timescale with a maximum of 10 years follow-up.

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 melanoma data and explore it.

```{webr-r}
#| autorun: true
library(biostat3)
localised <- subset(biostat3::melanoma, stage == "Localised") |>
    transform(## censoring for > 120 months
        death_cancer = ifelse( status == "Dead: cancer" & surv_mm <= 120, 1, 0),
        ## scale to years and truncate to 10 years        
        trunc_yy = pmin(surv_mm/12,10))

```

## (a) ## 
For the localised melanoma data with 10 years follow-up, plot the instantaneous cause-specific hazard for each calendar period.


```{webr-r}
#| autorun: true
library(bshazard)
library(tinyplot)
hazDiaDate <- lapply(levels(localised$year8594),
                     function(level) bshazard(Surv(trunc_yy,death_cancer)~1,
                                              data=subset(localised, year8594==level)) |>
                                     with(data.frame(time,hazard,lower.ci,upper.ci,
                                                     year8594=level))) |>
    do.call(what=rbind)
with(hazDiaDate, plt(hazard~time|year8594, ymin=lower.ci, ymax=upper.ci,
                     group=year8594,
                     type="ribbon",
                     xlab="Time since diagnosis (years)", 
                     ylab="Hazard"))
```

```{webr-r}
library(ggplot2)
# or using ggplot2
ggplot(hazDiaDate, aes(x=time, y=hazard, fill=year8594,
                       ymin=lower.ci, ymax=upper.ci)) +
    geom_line(aes(color=year8594)) +
    geom_ribbon(alpha=0.3) +
    xlab("Time since diagnosis (years)") +
    ylab("Hazard")

```

## (b) ## 
Now plot the instantaneous cause-specific hazard for each calendar period using a log scale for the y axis (use the option yscale(log)). What would you expect to see if a proportional hazards assumption was appropriate? Do you see it?

```{webr-r}
with(hazDiaDate, plt(hazard~time|year8594, ymin=lower.ci, ymax=upper.ci,
                     group=year8594,
                     type="ribbon",
                     log="y",
                     xlab="Time since diagnosis (years)", 
                     ylab="Hazard"))
```

## (c) ## 
Another graphical way of checking the proportional hazards assumption is to plot the log cumulative cause specific hazard function for each calendar period. These plots were not given extensive coverage in the lectures, so attempt this if you like or continue to part (d).

```{webr-r}
survfit1 <- survfit(Surv(trunc_yy,death_cancer)~year8594, data=localised)
plot(survfit1, col=1:2, fun=function(S) -log(-log(S)), log="x",
     xlab="log(time)", ylab="-log(H)")
legend("topright",legend=levels(localised$year8594),col=1:2,lty=1)

## or we can use
biostat3::survPHplot(Surv(trunc_yy,death_cancer)~year8594, data=localised)
```

## (d) ## 
Compare your estimated hazard ratio from part (a) with the one from a fitted Cox model with calendar period as the only explanatory variable. Are they similar?

```{webr-r}
cox1 <- coxph(Surv(trunc_yy, death_cancer==1) ~ year8594, data=localised)
summary(cox1)

```

## (e) ## 
Now fit a more complex model and use graphical methods to explore the assumption of proportional hazards by calendar period.

```{webr-r}
cox2 <- coxph(Surv(trunc_yy, death_cancer==1) ~ sex + year8594 + agegrp, data=localised)
summary(cox2)

cox2.phtest <- cox.zph(cox2, terms=FALSE) # for separate plots
print(cox2.phtest)
plot(cox2.phtest,var=2,resid=TRUE, se=TRUE, main="Schoenfeld residuals", ylim=c(-4,4))
```

What do you conclude?

## (f) ## 
Do plot the Schoenfeld residual plots for the variable agegrp. What are your
conclusions regarding the assumption of proportional hazards by age group?

```{webr-r}
cox2.phtest <- coxph(Surv(trunc_yy, death_cancer==1) ~ sex + year8594 + agegrp, 
                     data=localised) |>
    cox.zph(terms=FALSE)
par(mfrow=c(2,2))
for (var in 3:5)
    plot(cox2.phtest,var=var,resid=TRUE, se=TRUE, main="Schoenfeld residuals", ylim=c(-4,4))
```

## (g) ## 
Now formally test the assumption of proportional hazards.

```{webr-r}
coxph(Surv(trunc_yy, death_cancer==1) ~ sex + year8594 + agegrp, data=localised) |>
    cox.zph()
```

Are your conclusions from the test coherent with your conclusions from the graphical assessments?

## (h) ## 
Estimate separate age effects for the first two years of follow-up (and separate estimates for the remainder of the follow-up) while controlling for sex and period. Do the estimates for the effect of age differ between the two periods of follow-up? Write out the regression equation.

```{webr-r}
subset(localised, id %in% c(1,14), select=c(id, trunc_yy, death_cancer))
melanoma2p8Split <- survSplit(localised, cut=c(2), end="trunc_yy", start="start",
                              event="death_cancer", episode="fu") |>
    transform(fu = as.factor(fu))

##Tabulate ageband including risk_time

melanoma2p8Split |> subset(id %in% c(1,14), select=c(id, start, trunc_yy, death_cancer))

head(melanoma2p8Split)

cox2p8Split1 <- coxph(Surv(start, trunc_yy, death_cancer) ~ sex + year8594 + agegrp*fu,
                      data=melanoma2p8Split)
summary(cox2p8Split1)

cox2p8Split1b <- coxph(Surv(start, trunc_yy, death_cancer) ~ sex + year8594 + agegrp +
                           I(agegrp=="45-59" & fu=="2") + I(agegrp=="60-74" & fu=="2") +
                           I(agegrp=="75+" & fu=="2"), data=melanoma2p8Split)
summary(cox2p8Split1b)
```

We see effects of age (i.e., the hazard ratios) for the period 0–2 years subsequent to diagnosis along with the interaction effects. An advantage of the default parameterisation is that one can easily test the statistical significance of the interaction effects. Before going further, test whether the age*follow-up interaction is statistically significant (using a Wald and/or LR test).

## (i) ## 
Often we wish to see the effects of exposure (age) for each level of the modifier (time since diagnosis). That is, we would like to complete the table below with relevant hazard ratios. To get the effects of age for the period 2+ years after diagnosis, using the default parametrization, we must multiply the hazard ratios for 0–2 years by the appropriate interaction effect. Now let’s reparameterise the model to directly estimate the effects of age for each level of time since diagnosis.

              0–2 years   2+ years
------------- ----------- ----------
  Agegrp0-44     1.00        1.00
  Agegrp45-59    --          --
  Agegrp60-74    --          -- 
  Agegrp75+      --          --


```{webr-r}
coxph(Surv(start, trunc_yy, death_cancer) ~ sex + year8594 + fu + fu:agegrp, data=melanoma2p8Split) |> summary()

```

We can also use the `tt` argument in `coxph` for modelling for time-varying effects:

```{webr-r}
coxph(Surv(trunc_yy, death_cancer) ~ sex + year8594 + agegrp +
          tt(agegrp=="45-59") + tt(agegrp=="60-74") + tt(agegrp=="75+"),
      data=localised,
      tt = function(x, t, ...) x*(t>=2)) |> summary()
## The tt labels do not play nicely with lincom:(

cox2p8tvct <- coxph(Surv(trunc_yy, death_cancer) ~ sex + year8594 + agegrp + tt(agegrp),
                    data=localised,
                    tt = function(x, t, ...) cbind(`45-59`=(x=="45-59")*t,
                                                   `60-74`=(x=="60-74")*t,
                                                   `75+`=(x=="75+")*t))
summary(cox2p8tvct)
if (requireNamespace("car")) {
    lincom(cox2p8tvct, "agegrp75+",eform=TRUE)                   # t=0
    lincom(cox2p8tvct, "agegrp75+ + tt(agegrp)75+",eform=TRUE)   # t=1
    lincom(cox2p8tvct, "agegrp75+ + 2*tt(agegrp)75+",eform=TRUE) # t=2
}
```

Write out the regression equations for models `cox2p8Split2` and `cox2p8tvct`. Based on model `cox2p8tvct`, write out a formula for the hazard ratio for those aged 75 years and over compared with those aged less than 45 years as a function of time.



## (j) ## 

ADVANCED: Fit an analogous Poisson regression model. Are the parameter estimates similar? HINT: You will need to split the data by time since diagnosis. Calculate and plot the rates.


```{webr-r}
library(splines) # ns
library(dplyr)   # mutate, group_by, summarise
library(rstpm2)  # predictnl
time.cuts <- seq(0,10,length=100)
delta <- diff(time.cuts)[1]
## split and collapse
melanoma2p8Split2 <- survSplit(Surv(trunc_yy,death_cancer)~., data=localised,
                               cut=time.cuts, end="tstop", start="tstart",
                               event="death_cancer") |>
    mutate(fu=cut(tstop,time.cuts),
           mid=time.cuts[unclass(fu)]+delta/2) |>
    group_by(mid,sex,year8594,agegrp) |>
    summarise(pt=sum(tstop-tstart), death_cancer=sum(death_cancer),
              .groups="keep") |>
    mutate(age75 = (agegrp=="75+")+0)
poisson2p8tvc <- glm(death_cancer ~ sex + year8594 + agegrp + ns(mid,df=3) +
                         age75:ns(mid,df=3) + offset(log(pt)),
                     data=melanoma2p8Split2, family=poisson)
newdata <- data.frame(mid=seq(0,max(time.cuts),length=100), year8594="Diagnosed 85-94",
                      sex="Male", agegrp="75+", age75=1, pt=1)
predictnl(poisson2p8tvc,
          fun=function(fit,newdata)
              predict(fit, newdata) -
              predict(fit, transform(newdata, agegrp='0-44', age75=0)),
          newdata=newdata) |>
    cbind(newdata) |>
    transform(RateRatio=exp(fit),
              lower.ci=exp(fit-1.96*se.fit),
              upper.ci=exp(fit+1.96*se.fit)) |>
    with(plt(mid, RateRatio, type="ribbon", xlab="Time since diagnosis (months)",
             ymin=lower.ci, ymax=upper.ci,
             ylab="Rate ratio", main="Ages 75+ compared with ages 0-44 years"))

newdata <- expand.grid(mid=c(0.5,2), year8594="Diagnosed 85-94",
                       sex="Male", agegrp=levels(localised$agegrp), pt=1) |>
    transform(age75=(agegrp=="75+")+0)
predictnl(poisson2p8tvc,
          fun=function(fit,newdata)
              predict(fit, newdata) -
              predict(fit, transform(newdata, agegrp='0-44', age75=0)),
          newdata=newdata) |>
    cbind(newdata) |>
    transform(RateRatio=exp(fit),
              lower.ci=exp(fit-1.96*se.fit),
              upper.ci=exp(fit+1.96*se.fit)) |>
    subset(select=c(agegrp, mid, RateRatio, lower.ci, upper.ci))
```

