# Biostatistics III in R

#### Code by Xing-Rong Liu and Mark Clements

## Exercise 28. Flexible parametric survival models in R

Flexible parametric survival models use a smooth function to model for a transformation of survival. These models are also called Royston-Parmar models. For our model extensions, we also use the term generalized survival models.

Note that there are some differences in the Stata implementation (`stpm2`) and our implementation in R (`rstpm2::stpm2` and `rstpm::pstpm2`). For the parametric models, Stata uses restricted cubic splines using truncated power basis, while we use B-splines. Moreover, we allow for most smooth functions to be used for the baseline function.

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(rstpm2)  # for the flexible parametric model
    library(dplyr)   # for data manipulation

We start by reading the data and then define a 1/0 varible for the events that we are interested in.

    data(melanoma)

    ## Extract subset and create 0/1 outcome variables
    melanoma0 <- melanoma %>% filter(stage=="Localised") %>%
                 mutate(event = ifelse(status=="Dead: cancer" & surv_mm<120, 1, 0),
                        time = pmin(120, surv_mm)/12,
                        agegrp1 = (agegrp=="0-44")+0,  # used by time-dependent effect
                        agegrp2 = (agegrp=="45-59")+0, # used by time-dependent effect
                        agegrp3 = (agegrp=="60-74")+0, # used by time-dependent effect
                        agegrp4 = (agegrp=="75+")+0)   # used by time-dependent effect

## (a)

Flexible parametric model with df=4.

    ## (a) Flexible parametric model with df=4

    fpma <- stpm2(Surv(time,event) ~ year8594, data=melanoma0, df=4)
    exp(cbind(IRR=coef(fpma),confint(fpma)))

    ## Profiling...

    ## Warning in .local(object, parm, level, ...): non-monotonic spline fit
    ## to profile (nsx(log(time), df = 4)2): reverting from spline to linear
    ## approximation

    summary(fpma)

Making a cox model for comparison.

    cox <- coxph(Surv(time, event) ~ year8594,
                 data=melanoma0) # year8594 is a categorical variable
    summary(cox)

The hazard ratio, 95% confidence interval and statistical significance are very similar to the Cox model.

## (b)

Prediction and plots of survival and hazard by calendar period.

    ## (b) Prediction and plots of survival and hazard by calendar period
    years <- levels(melanoma0$year8594)
    alegend <- function() legend("topright", legend=years, col=1:2, lty=1, bty="n")

    plot(fpma,newdata=data.frame(year8594=years[1]),
         xlab="Time since diagnosis (years)")
    lines(fpma,newdata=data.frame(year8594=years[1]), ci=FALSE) # repeated
    alegend()

Localised skin melanoma. Predicted survival functions from a exible parametric model.

    plot(fpma,newdata=data.frame(year8594=years[1]), type="haz",
         xlab="Time since diagnosis (years)", ylab="Hazard")
    lines(fpma,newdata=data.frame(year8594=years[2]), type="haz", col=2)
    alegend()

Localised skin melanoma. Predicted hazard functions from a exible parametric model.

## (c)

Plotting the hazards on the log scale. There is a constant difference as the predictions are from a proportional hazards model and a multiplicative effect becomes additive on the log scale.

    ## (c) hazards on log scale, adding log="y"
    plot(fpma,newdata=data.frame(year8594=years[1]), type="haz",
         xlab="Time since diagnosis (years)", log="y", ci=FALSE,
         ylab="Hazard (log scale)")
    lines(fpma,newdata=data.frame(year8594=years[2]), type="haz", col=2)
    alegend()

Localised skin melanoma. Predicted hazard functions on log scale from a exible parametric model.

## (d)

Investigate DF for the baseline with criteria AIC and BIC.

    summary(fpma)
    aicc <- bicc <- beta <- se <- rep(NULL,6)
    BIC <- function(object, nknots){
        -2 * as.numeric(bbmle::logLik(object)) + log(sum(melanoma0$event)) * (1 + nknots)
    }

    for (i in 1:6 ) {
      fitaic <- stpm2(Surv(time, event) ~ year8594, data=melanoma0, df=i)
      aicc[i] <- AIC(fitaic)
      bicc[i] <- BIC(fitaic, i)
      beta[i] <- as.numeric(coef(fitaic)[2])
      se[i] <- as.numeric(sqrt(fitaic@vcov[2,2]))
    }
    rbind(beta=beta, se=se, aicc=aicc, bicc=bicc)

The log hazard ratios (and hence the hazard ratios) from 2 df and up are similar and for 3 df they are very similar. The main difference is for 1 df, which is equivalent to a Weibull model. The Weibull model enforces a monotonic hazard function and as the hazard function in the melanoma data has a turning point it is clearly inappropriate.

The lowest AIC is for the model with 5 df and for the BIC it is the model with 2 df. The penalty term in the AIC is twice the number of parameters (2 x k ) whereas in the BIC it is ln( D ) > k where D is the number of events. Since ln( D ) > k the BIC penalizes extra parameters much more strongly than AIC. Since we have a large data set and there are no disadvantages to including extra parameters we would use 5df for the baseline hazard.

## (e)

Plots of predicted baseline survivals and hazards with df=1:6.

    ## Baseline survival
    fitaic0 <- stpm2(Surv(time, event) ~ year8594, data=melanoma0, df=6)
    plot(fitaic0,newdata=data.frame(year8594=years[1]), lty=6, ci=F,
         xlab="Time since diagnosis (years)")

    dfs <- c("s_df1","s_df2","s_df3","s_df4","s_df5","s_df6")
    for (i in 1:5 ) {
      fitaic <- stpm2(Surv(time, event) ~ year8594, data=melanoma0, df=i)
      plot(fitaic,newdata=data.frame(year8594=years[1]), add=TRUE,lty=i,
           xlab="Time since diagnosis (years)")
    }
    legend("topright", legend=dfs[1:6], lty=1:6)

    ## Baseline hazard
    fitaic1 <- stpm2(Surv(time, event) ~ year8594, data=melanoma0, df=6)
    plot(fitaic1,newdata=data.frame(year8594=years[1]), lty=6, type="haz",
         ci=F, xlab="Time since diagnosis (years)", ylab="Hazard")

    for (i in 1:5 ) {
      fitaic <- stpm2(Surv(time, event) ~ year8594, data=melanoma0, df=i)
      plot(fitaic,newdata=data.frame(year8594=years[1]), add=TRUE, lty=i,
           type="haz", xlab="Time since diagnosis (years)", ylab="Hazard")
    }
    dfs2 <- c("h_df1","h_df2","h_df3","h_df4","h_df5","h_df6")
    legend("topright", legend=dfs2[1:6], lty=1:6)

With the exception of 1 df (the Weibull model), the survival and hazard functions show similar shapes, so as long we have enough knots our conclusions would be very similar.

## (f)

Adjust for sex and age (in categories).

    fpmf <- stpm2(Surv(time, event) ~ sex + year8594 + agegrp,
                  data=melanoma0, df=4)

    beta <- coef(fpmf)[2:6] ## log(HR)
    se <- sqrt(diag(fpmf@vcov[2:6,2:6]))            # standard errors of beta
    z <- beta/se                                    # z-scores
    confin <- cbind(beta - 1.96*se, beta + 1.96*se) # 95% confidence interval of beta=log(HR)
    p_value <- 1-pchisq((beta/se)^2, 1)             # Wald-type test of beta

    ## Summaries of Beta and HR
    cbind(beta = beta, se = se, z = beta/se, p_value = p_value) # summaries of beta
    cbind(HR = exp(beta), lower.95 = exp(confin)[,1], upper.95 =  exp(confin)[,2])

    ## To test the joint significance of categorized age with Wald-type test
    ## Joint H0:
    ##          beta_(agegrp45-59) = 0
    ##          beta_(agegrp60-74) = 0
    ##          beta_(agegrp75+) = 0
    statistic <- t(coef(fpmf)[4:6]) %*% solve(vcov(fpmf)[4:6,4:6]) %*% coef(fpmf)[4:6]
    p_value <- 1 - pchisq(statistic, 3)
    cbind(chi2_statistics =statistic, P_value = p_value)

    ## To test the overall effect of age with LR test
    fpmf2 <- stpm2(Surv(time,event) ~ sex + year8594 , data=melanoma0, df=4)
    anova(fpmf, fpmf2)

The estimates are similar to those obtained from the Cox model. The Wald test yields a very highly significant result, which is similar to that obtained from the comparable test for the Cox model.

## (g)

The estimates are so similar because very similar models are being fitted with exactly the same covariates. The two models differ only in the manner in which they account for the baseline hazard. In the Cox model it is assumed arbitrary and not directly estimated. In the flexible parametric model the baseline hazard is modelled using splines. The 5 df spline allows sufficient flexibility to make the model estimates virtually identical.

## (h)

Change to time-varying effect of agegrp2:4.

    ## (h) Change to time-varying effect of agegrp2:4
    ## NB: including main effect of agerp
    fpmh <- stpm2(Surv(time,event) ~ sex + year8594 + agegrp2 + agegrp3 + agegrp4,
                  data=melanoma0, tvc=list(agegrp2 = 2, agegrp3 = 2, agegrp4 = 2),
                  smooth.formula=~ nsx(log(time),df=4)  )
    ## NB: no main effect of agegrp
    fpmh1 <- stpm2(Surv(time,event) ~ sex + year8594 ,
                   data=melanoma0, tvc=list(agegrp2 = 2, agegrp3 = 2, agegrp4 = 2),
                   smooth.formula=~ nsx(log(time),df=4)  )

    beta <- coef(fpmh)[2:6] ## log(HR)
    se <- sqrt(diag(fpmh@vcov[2:6,2:6])) ## standard errors of beta
    z <- beta/se ## z-scores
    confin <- cbind(beta - 1.96*se, beta + 1.96*se) ## 95% confidence interval of beta=log(HR)
    p_value <- 1-pchisq((beta/se)^2, 1) ## Wald-type test of beta

    ## Summaries of Beta and HR
    cbind(beta = beta, se = se, z = beta/se, p_value = p_value) ## summaries of beta
    cbind(HR = exp(beta), lower.95 = exp(confin)[,1], upper.95 =  exp(confin)[,2])

    ## LR test comparing fpmh (non-PH for agegrp2:4) with fpmf(PH for agegrp2:4)
    anova(fpmh, fpmf)
    anova(fpmh1, fpmf)

    ## I investigated the non-proportional effect of age with penalized models
    ## here, sp is the optimal smoothing parameters estimated from models without sp argument
    pfit0 <- pstpm2(Surv(time,event) ~ sex + year8594 + agegrp2 + agegrp3 + agegrp4,
                    smooth.formula=~s(log(time)), data=melanoma0, sp=0.1359685)

    ## The time-dependent effects including linear forms of age groups
    pfit1 <- pstpm2(Surv(time,event) ~ sex + year8594,
                    smooth.formula=~s(log(time)) + s(log(time),by=agegrp2) +
                                    s(log(time),by=agegrp3) + s(log(time),by=agegrp4),
                    data=melanoma0, sp=c( 0.1429949, 1.6133966, 1.3183117, 1.9958815))
    anova(pfit1, pfit0)## the results also suggest there is strong evidence

There is strong evidence of a non-proportional effect of age.

## (i)

Plot of baseline hazard with fpmh. This baseline is for the youngest age group who are male and diagnosed in 1975-1984, i.e, when all the covariates are equal to zero.

    ## (i) Plot of baseline hazard with fpmh
    sexs <- levels(melanoma$sex)
    years <- levels(melanoma$year8594)
    agegrps <- levels(melanoma$agegrp)
    par(mfrow=c(1,1))
    newdata1 <- data.frame(sex=sexs[1],year8594=years[1],agegrp2=0, agegrp3=0, agegrp4=0)
    plot(fpmh,newdata=newdata1, xlab="Time since diagnosis (years)", ylab="Hazard",
         type="haz", ci=FALSE, ylim=c(0,0.051), rug=FALSE)
    legend("topright", legend=c(sexs[1], years[1], agegrps[1]), lty=1, bty="n")

## (j)

Hazard ratios between age groups.

    plot(fpmh, newdata=data.frame(sex=sexs[1],year8594=years[1],agegrp2=0, agegrp3=0, agegrp4=0),
         type="hr", rug=FALSE, line.col=1, ci=FALSE,log="y", ylim=c(1,500), lty=1,
         exposed=function(data) transform(data,sex=sexs[1],year8594=years[1],agegrp2=1, agegrp3=0, agegrp4=0),
         xlab="Time since diagnosis (years)",ylab="Hazards ratio")
    plot(fpmh, newdata=data.frame(sex=sexs[1],year8594=years[1],agegrp2=0, agegrp3=0, agegrp4=0),
         type="hr", rug=FALSE, line.col=1, ci=FALSE, add=T,log="y", lty=2,
         exposed=function(data) transform(data,sex=sexs[1],year8594=years[1],agegrp2=0, agegrp3=1, agegrp4=0),
         xlab="Time since diagnosis (years)",ylab="Hazards ratio")
    plot(fpmh, newdata=data.frame(sex=sexs[1],year8594=years[1],agegrp2=0, agegrp3=0, agegrp4=0),
         type="hr", rug=FALSE, line.col=1, ci=FALSE,add=T,log="y", lty=3,
         exposed=function(data) transform(data,sex=sexs[1],year8594=years[1],agegrp2=0, agegrp3=0, agegrp4=1),
         xlab="Time since diagnosis (years)",ylab="Hazards ratio")
    legend("topright", legend=agegrps[2:4], lty=1:3)

The hazard ratios decrease as a function of follow-up time. The hazard ratio is so high during the early years of follow-up because the hazard in the reference group is close to zero.

    plot(fpmh, newdata=data.frame(sex=sexs[1],year8594=years[1],agegrp2=0, agegrp3=0, agegrp4=0),
         type="hr", rug=FALSE, line.col=1, ci=T, log="y",
         exposed=function(data) transform(data,sex=sexs[1],year8594=years[1],agegrp2=0, agegrp3=0, agegrp4=1),
         xlab="Time since diagnosis (years)",ylab="Hazards ratio")

    ## Warning in xy.coords(x, y, xlabel, ylabel, log = log): 9 y values <= 0
    ## omitted from logarithmic plot

The hazard ratio for the oldest age group with a 95% confidence interval.

## (k)

The hazard difference for the oldest group.

    plot(fpmh,newdata=data.frame(sex=sexs[1],year8594=years[1],agegrp2=0, agegrp3=0, agegrp4=0),
         type="hdiff", rug=FALSE, line.col=1, ci=T,
         exposed=function(data) transform(data,sex=sexs[1],year8594=years[1],agegrp2=0, agegrp3=0, agegrp4=1),
         xlab="Time since diagnosis (years)",ylab="Hazards difference")

Localised skin melanoma. Predicted difference in hazard rates (with a 95% confidence interval) for oldest age group from a fexible parametric model (other covariates are set to zero).

The hazard difference is small early on, despite the hazard ratio being large, because the underlying hazard is so low.

## (l)

The survival difference for the oldest group.

    plot(fpmh,newdata=data.frame(sex=sexs[1],year8594=years[1],agegrp2=0, agegrp3=0, agegrp4=0),
         type="sdiff", rug=FALSE, line.col=1, ci=T,
         exposed=function(data) transform(data, sex=sexs[1],year8594=years[1],agegrp2=0, agegrp3=0, agegrp4=1),
         xlab="Time since diagnosis (years)",ylab="Survivals difference")

Localised skin melanoma. Predicted difference in survival functions (with a 95% confidence interval) for oldest age group from a flexible parametric model (for females diagnosed in 1984-1994).

## (m)

Fit models with 1, 2 and 3 df for time-dependent effect of age.

    aicc <- bicc <- beta <- se <- rep(1,3)
    nevent <- sum(melanoma0$event)
    BIC <- function(object, nknots){
      -2 * as.numeric(bbmle::logLik(object)) + log(nevent) * (1 + nknots)
    }
    for (i in 1:3 ) {
      fitdf <- stpm2(Surv(time,event) ~ sex + year8594 + agegrp2 + agegrp3 + agegrp4,
                      data=melanoma0, tvc=list(agegrp2 = i, agegrp3 = i, agegrp4 = i),
                      smooth.formula=~ nsx(log(time),df=4)  )
      aicc[i] <- AIC(fitaic)
      bicc[i] <- BIC(fitaic, i)
    }
    cbind(aicc=aicc, bicc=bicc)

    ## PLots with different df
    fitdf1 <- stpm2(Surv(time,event) ~ sex + year8594 + agegrp2 + agegrp3 + agegrp4, data=melanoma0,
                    tvc = list(agegrp2 = 1, agegrp3 = 1, agegrp4 = 1), smooth.formula=~ nsx(log(time),df=4))

    plot(fitdf1, newdata = data.frame(sex=sexs[1], year8594=years[1], agegrp2=0, agegrp3=0, agegrp4=0),
         type="hr", rug=FALSE, line.col=1, ci=F, log="y", lty=1,
         exposed=function(data) transform(data, sex = sexs[1],year8594=years[1],
                                            agegrp2=0, agegrp3 = 0, agegrp4=1),
         xlab="Time since diagnosis (years)", ylab="Hazards ratio")

    for (i in 2:3 ) {
        fitdf <- stpm2(Surv(time,event) ~ sex + year8594 + agegrp2 + agegrp3 + agegrp4,
                       data=melanoma0, tvc=list(agegrp2 = i, agegrp3 = i, agegrp4 = i),
                       smooth.formula=~ nsx(log(time),df=4))

        plot(fitdf, newdata=data.frame(sex=sexs[1],year8594=years[1],agegrp2=0, agegrp3=0, agegrp4=0),
             type="hr", rug=FALSE, line.col=1, ci=F, log="y", add=T, lty=i,
             exposed=function(data) transform(data, sex=sexs[1],year8594=years[1],agegrp2=0, agegrp3=0, agegrp4=1),
             xlab="Time since diagnosis (years)",ylab="Hazards ratio")
    }
    legend("topright", legend=c("1 df", "2 df", "3 df"), lty=1:3)

AIC selects 2 df and BIC selects 1 df. As discussed in the previous part, the BIC imposes a stronger penalty on aditional parameters. The fitted time-dependent effects are similar. We would suggest 2 df for the time-varying effect of age.