# Biostatistics III in R

#### Code by Annika Tillander, Johan Zetterqvist and Mark Clements

# Exercise 9. Localised melanoma: modelling cause-specific mortality using Cox regression

In exercise 7 we modelled the cause-specific mortality of patients diagnosed with localised melanoma using Poisson regression. We will now model cause-specific mortality using Cox regression and compare the results to those we obtained using the Poisson regression model.

To fit a Cox proportional hazards model (for cause-specific survival) with calendar period as the only explanatory variable, the following commands can be used. Note that we are censoring all survival times at 120 months (10 years) in order to facilitate comparisons with the Poisson regression model in exercise 7.

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(splines)   # ns (recommended package)

    ## Read melanoma data, select subcohorts and create a death indicator
    melanoma.l <- biostat3::melanoma %>%
        filter(stage=="Localised") %>%
        mutate(death_cancer = as.numeric(status=="Dead: cancer"))

    melanoma.l2 <-    mutate(melanoma.l,
                             ## Create a death indicator (only count deaths within 120 months)
                             death_cancer = death_cancer * as.numeric( surv_mm <= 120),
                             ## Create a new time variable
                             surv_mm = pmin(surv_mm, 120))

## (a)

Interpret the estimated hazard ratio, including a comment on statistical significance. Write the linear predictor using pen and paper. Draw a graph showing the shape of the fitted hazard rates for the two calendar periods, also indicating the distance between the curves (HINT: draw on the log hazard scale if you find it easier). Compare this to how you drew the graph in exercise 7h.

    summary( coxfit9a <- coxph(Surv(surv_mm, death_cancer) ~ year8594,
                               data = melanoma.l2) )

## (b)

(This part is more theoretical and is not required in order to understand the remaining parts.)

R reports a Wald test of the null hypothesis that survival is independent of calendar period. The test statistic (and associated P-value), is reported in the table of parameter estimates. Under the null hypothesis, the test statistic has a standard normal (Z) distribution, so the square of the test statistic will have a chi square distribution with one degree of freedom.

R also reports a likelihood ratio test statistic of the null hypothesis that none of the parameters in the model are associated with survival. In general, this test statistic will have a chi-square distribution with degrees of freedom equal to the number of parameters in the model. For the current model, with only one parameter, the test statistic has a chi square distribution with one degree of freedom.

Compare these two test statistics with each other and with the log rank test statistic (which also has a chi square distribution) calculated in question 3c (you should, however, recalculate the log rank test since we have restricted follow-up to the first 10 years in this exercise). Would you expect these test statistics to be similar? Consider the null and alternative hypotheses of each test and the assumptions involved with each test.

    summary( coxfit3c <- coxph(Surv(surv_mm, death_cancer) ~ year8594,
                               data = melanoma.l) )

## (c)

Now include sex and age (in categories) in the model. Write the linear predictor using pen and paper.

1. Interpret the estimated hazard ratio for the parameter labelled agegrp45-59, including a comment on its statistical significance.
2. Is the effect of calendar period strongly confounded by age and sex? That is, does the inclusion of sex and age in the model change the estimate for the effect of calendar period?
3. Perform a Wald test of the overall effect of age and interpret the results.

    summary( coxfit9c <- coxph(Surv(surv_mm, death_cancer) ~ year8594 + sex + agegrp,
                               data = melanoma.l2) )

    ## Test if the effect of age is significant
    ## Use a Wald test with the car package
    library(car)
    linearHypothesis(coxfit9c,c("agegrp45-59 = 0","agegrp60-74 = 0","agegrp75+ = 0"))

## (d)

Perform a likelihood ratio test of the overall effect of age and interpret the results.

    ## Test if the effect of age is significant
    ## Use an LR test with the anova function
    ## The model without agegrp
    summary( coxfit9d <- coxph(Surv(surv_mm, death_cancer) ~ year8594 + sex,
                               data = melanoma.l2) )
    ## LR test
    anova(coxfit9c, coxfit9d)

Compare your findings to those obtained using the Wald test. Are the findings similar? Would you expect them to be similar?

## (e)

The model estimated in question 9c is similar to the model estimated in question 7i.

1. Both models adjust for sex, year8594, and i.agegrp but the Poisson regression model in question 7i appears to adjust for an additional variable (`fu`). Is the Poisson regression model adjusting for an additional factor? Explain.
2. Would you expect the parameter estimate for sex, period, and age to be similar for the two models? Are they similar?
3. Do both models assume proportional hazards? Explain.

    summary( coxfit9e <- coxph(Surv(surv_mm, death_cancer) ~ year8594 + sex + agegrp,
                               data = melanoma.l2) )

    ## Split follow up by year
    melanoma.spl <- survSplit(melanoma.l2, cut=12*(0:10), end="surv_mm", start="start",
                               event="death_cancer")

    ## Calculate persontime and
    ## recode start time as a factor
    melanoma.spl <- mutate(melanoma.spl,
                           pt = surv_mm - start,
                           fu = as.factor(start) )

    ## Run Poisson regression
    summary(poisson9e <- glm( death_cancer ~ fu + year8594 + sex + agegrp + offset( log(pt) ),
                             family = poisson,
                             data = melanoma.spl ))

    eform(poisson9e)
    summary(coxfit9e)

## (f)

ADVANCED: By splitting at each failure time we can estimate a Poisson regression model that is identical to the Cox model.

    sfit9f <- survfit(Surv(surv_mm, event=death_cancer) ~ 1,
                      data = melanoma.l2 )
    ## Have a look at the fitted object
    str(sfit9f, 1)

    head(sfit9f$time)

    ## Split follow up by year
    melanoma.spl2 <- survSplit(melanoma.l2, cut=sfit9f$time, end="surv_mm", start="start",
                               event="death_cancer")

    ## Calculate persontime and
    ## recode start time as a factor
    melanoma.spl2 <- mutate(melanoma.spl2,
                            pt = surv_mm - start,
                            fu = as.factor(start) )

    ## Collapse
    melanoma.spl3 <- melanoma.spl2 %>%
        group_by(fu,year8594,sex,agegrp) %>%
        summarise(pt=sum(pt), death_cancer=sum(death_cancer)) %>%
        data.frame

    ## Run Poisson regression
    poisson9f <- glm( death_cancer ~ fu + year8594 + sex + agegrp + offset( log(pt) ),
                     family = poisson,
                     data = melanoma.spl3 )

    ## IRR
    coef9f <- eform(poisson9f)
    ## We are not interested in nuisance parameters fu1, fu2, etc
    npar <- nrow(coef9f)
    pars <- (npar-4):npar
    coef9f[pars,]

    summary(coxfit9e)

## (g)

ADVANCED: Split the data finely (e.g., 3-month intervals) and model the effect of time using a restricted cubic spline with Poisson regression. Also plot the rate and survival.

    ## split and collapse
    cuts.splines <- seq(0,max(sfit9f$time),by=3)
    mid.splines <- cuts.splines + 1.5
    melanoma.spl4 <-
        survSplit(Surv(surv_mm,death_cancer)~., data=melanoma.l2, cut=cuts.splines,
                  start="tstart", end="tstop") %>%
        mutate(cut=cut(tstop,cuts.splines),
               mid=mid.splines[unclass(cut)]+1.5) %>%
        group_by(mid,year8594,sex,agegrp) %>%
        summarise(pt=sum(tstop-tstart), death_cancer=sum(death_cancer))

    poisson9g <- glm( death_cancer ~ ns(mid,df=3) + year8594 + sex + agegrp + offset( log(pt) ),
                     family = poisson,
                     data = melanoma.spl4 )
    summary(poisson9g)
    eform(poisson9g)

    ## quick approach: use the effects package
    library(effects)

    ## Loading required package: carData

    ##
    ## Attaching package: 'carData'

    ## The following objects are masked from 'package:car':
    ##
    ##     Guyer, UN, Vocab

    ## lattice theme set by effectsTheme()
    ## See ?effectsTheme for details.

    plot(Effect("mid", poisson9g))

    ## 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
    times <- seq(0,max(cuts.splines),length=1000)
    delta <- diff(times)[1]
    newdata <- data.frame(mid=times+delta/2, year8594="Diagnosed 85-94",
                          sex="Male", agegrp="45-59",
                          pt=1)
    ## predict rates and 95% CI
    pred <- predict(poisson9g, newdata=newdata, se.fit=TRUE)
    predrate <- exp(pred$fit)
    ci <- with(pred, exp(cbind(fit-1.96*se.fit, fit+1.96*se.fit)))
    ## plot
    matplot(newdata$mid, ci, type="n", xlab="Time since diagnosis (months)",
            ylab="Rate", main="Males aged 45-59 years diagnosed 1985-94")
    polygon.ci(newdata$mid, ci)
    lines(newdata$mid, predrate)

    ## predict survival and 95% CI
    library(rstpm2)

    ##
    ## Attaching package: 'rstpm2'

    ## The following objects are masked from 'package:biostat3':
    ##
    ##     colon, eform, popmort

    ## The following object is masked from 'package:survival':
    ##
    ##     colon

    logcumhaz <- rstpm2::predictnl(poisson9g,
              fun=function(fit,newdata) log(cumsum(delta*predict(fit, newdata, type="response"))),
              newdata=newdata)
    surv <- exp(-exp(logcumhaz$Estimate))
    ci <- with(logcumhaz, exp(-exp(cbind(Estimate-1.96*SE, Estimate+1.96*SE))))

    ## plot
    matplot(newdata$mid, ci, type="n", xlab="Time since diagnosis (months)",
            ylab="Survival", main="Males aged 45-59 years diagnosed 1985-94")
    polygon.ci(newdata$mid, ci)
    lines(newdata$mid, surv)