# Biostatistics III in R

#### Code by Annika Tillander, Andreas Karlsson and Mark Clements

# 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.

    library(biostat3)
    library(dplyr)    # for data manipulation
    library(ggplot2)

Load melanoma data and explore it.

    localised <- dplyr::filter(melanoma, stage == "Localised") %>%
        dplyr::mutate(death_cancer = ifelse( status == "Dead: cancer" & surv_mm <= 120, 1, 0), #censoring for > 120 months
               trunc_yy = pmin(surv_mm/12,10))  #scale to years and truncate to 10 years

## (a)

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

    # Using muhaz2 to smooth the Kaplan-Meier hazards by strata
    hazDiaDate <- muhaz2(Surv(trunc_yy,death_cancer)~year8594, data=localised)
    hazDiaDateDf <- as.data.frame(hazDiaDate)

    ## Comparing hazards
    plot(hazDiaDate, haz.scale=1000,
         xlab="Time since diagnosis (years)",
         ylab="Hazard per 1000 person-years")
    # or using ggplot2
    ggplot(hazDiaDateDf, aes(x=x, y=y*1000, colour= strata)) + geom_line() +
        xlab("Time since diagnosis (years)") +
        ylab("Hazard per 1000 person-years")

## (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?

    ## Comparing hazards on a log scales
    plot(hazDiaDate, log="y")

## (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).

    ## Calculating -log cumulative hazards per strata
    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?

    # Cox regression with time-since-entry as the timescale
    # Note that R uses the Efron method for approximating the likelihood in the presence of ties
    # whereas Stata (and some other software) use the Breslow method
    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.

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

    ## Plot of the scaled Schoenfeld residuals for calendar period 1985–94.
    ## The smooth line shows the estimated log hazard ratio as a function of time.
    cox2.phtest <- cox.zph(cox2, transform="identity") #Stata appears to be using 'identity'
    plot(cox2.phtest[2],resid=TRUE, se=TRUE, main="Schoenfeld residuals", ylim=c(-4,4))

What do you conclude?

## (f)

Do part (a)–(e) but now for the variable agegrp. What are your conclusions regarding the assumption of proportional hazards? No written solutions for this part. No written solutions for this part.

## (g)

Now formally test the assumption of proportional hazards.

    ## The results from the previous proportional hazards assumption test
    print(cox2.phtest)

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?

    melanoma2p8Split <- survSplit(localised, cut=c(2), end="trunc_yy", start="start",
                                  event="death_cancer", episode="fu") %>%
        mutate(fu = as.factor(fu))

    ##Tabulate ageband including risk_time
    melanoma2p8Split %>% select(id, start, trunc_yy) %>% filter(id<=3) %>% arrange(id, trunc_yy)

    head(melanoma2p8Split)

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

    ## Warning in coxph(Surv(start, trunc_yy, death_cancer) ~ sex + year8594 + : X
    ## matrix deemed to be singular; variable 6

    summary(cox2p8Split1)

    cox2p8Split1b <- coxph(Surv(start, trunc_yy, death_cancer) ~ sex + year8594 + agegrp + 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+ | – | –

    cox2p8Split2 <- coxph(Surv(start, trunc_yy, death_cancer) ~ sex + year8594 + fu + fu:agegrp, data=melanoma2p8Split)

    ## Warning in coxph(Surv(start, trunc_yy, death_cancer) ~ sex + year8594 + : X
    ## matrix deemed to be singular; variable 3

    summary(cox2p8Split2)

    ## Alternative approach using tt():
    ## http://cran.r-project.org/web/packages/survival/vignettes/timedep.pdf
    cox2p8tvc2 <- coxph(Surv(trunc_yy, death_cancer) ~ sex + year8594 + agegrp + tt(agegrp), data=localised,
                       tt = function(x, t, ...) (x=="75+")*(t>=2))
    summary(cox2p8tvc2)

    cox2p8tvct <- coxph(Surv(trunc_yy, death_cancer) ~ sex + year8594 + agegrp + tt(agegrp), data=localised,
                       tt = function(x, t, ...) (x=="75+")*t)
    summary(cox2p8tvct)
    lincom(cox2p8tvct, "agegrp75+ + tt(agegrp)",eform=TRUE)

    cox2p8tvclogt <- coxph(Surv(trunc_yy, death_cancer) ~ sex + year8594 + agegrp + tt(agegrp), data=localised,
                        tt = function(x, t, ...) (x=="75+")*log(t))
    summary(cox2p8tvclogt)

## (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.

    library(splines)
    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)) %>%
        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)

    ## utility function to draw a confidence interval
    polygon.ci <- function(time, interval, col="lightgrey")
        polygon(c(time,rev(time)), c(interval[,1],rev(interval[,2])), col=col, border=col)
    ## define exposures
    newdata <- data.frame(mid=seq(0,max(time.cuts),length=100), year8594="Diagnosed 85-94",
                          sex="Male", agegrp="75+", age75=1, pt=1)

    library(rstpm2)
    logirr <- rstpm2::predictnl(poisson2p8tvc,
              fun=function(fit,newdata) predict(fit, newdata) -
                            predict(fit, transform(newdata, agegrp='0-44', age75=0)),
              newdata=newdata)
    pred <- exp(logirr$fit)
    ci <- exp(confint(logirr))
    ## plot
    matplot(newdata$mid, ci, type="n", xlab="Time since diagnosis (months)",
            ylab="Rate ratio", main="Ages 75+ compared with ages 0-44 years")
    polygon.ci(newdata$mid, ci)
    lines(newdata$mid, pred)