Simulated example from a Proportional Excess Hazards model

Keywords: General Hazard Structure, Accelerated Failure Time, Accelerated Hazards, Excess Hazard, Exponentiated Weibull distribution, Net Survival, Proportional Hazards.

In cancer epidemiology using population-based registry data, the quantity of interest is usually the excess mortality hazard instead of the overall mortality hazard. The basic idea behind excess mortality hazard models consists of decomposing the mortality hazard function associated to an individual, \(h(t;{\bf x})\), as the sum of an excess mortality hazard, \(h_{\text{E}}(t;{\bf x})\), and the general population mortality hazard \(h_{\text{P}}(t;{\bf z})\). The general population mortality hazard represents the other-causes hazard in our population of interest, assuming that the contribution in the general population mortality hazard of the disease of interest is small compared to all others. This excess mortality hazard is interpreted as the mortality hazard that could be due, directly or indirectly, to the cancer under study. This is [3], \[\begin{eqnarray} \text{observed hazard = population hazard + excess hazard.} \end{eqnarray}\] Formally, this can be written as: \[\begin{eqnarray}\label{HazardAddModel} h(t;{\bf x}) = h_{\text{P}}(\text{age}+t;\text{year}+t;{\bf z}) + h_{\text{E}}(t;{\bf x}), \end{eqnarray}\]

where \({\bf x}\) and \({\bf z}\) are vectors of covariables, \({\bf z}\) typically corresponds to a subset of covariables of \({\bf x}\), and ‘age’ is the age at diagnosis and ‘year’ is the year of diagnosis. The population hazard \(h_{\text{P}}(\cdot;{\bf z})\) is typically obtained from national life tables based on the available sociodemographic characteristics (\({\bf z}\) in addition to age and year. A survival function derived from the excess hazard \(h_{\text{E}}(\cdot;{\bf x})\) is called the . Net survival represents a useful way of reporting the probability of survival of cancer patients since this allows for a fairer comparison of cancer survival between different populations or countries.

In practice, the excess hazard function is typically modelled using a proportional hazards structure (see below). [2] considered modelling the excess hazard function using more general hazard-based structures. In particular, they study the use of the following hazard structures:

Proportional Hazards model (PH). \[\begin{eqnarray}\label{PH} h_{\text{E}}^{\text{PH}}\left(t;{\bf x}_j\right) &=& h_0(t)\exp\left({\bf x}_j^T\beta\right)\\ H_{\text{E}}^{\text{PH}}\left(t;{\bf x}_j\right) &=& H_0(t)\exp\left({\bf x}_j^T\beta\right)\nonumber \end{eqnarray}\] Accelerated Hazards model (AH). \[\begin{eqnarray*}\label{AH} h_{\text{E}}^{\text{AH}}\left(t;{\bf x}_j\right) &=& h_0\left(t \exp({\bf x}_j^T\beta)\right)\\ H_{\text{E}}^{\text{AH}}\left(t;{\bf x}_j\right) &=& H_0\left(t \exp({\bf x}_j^T\beta)\right)\exp(-{\bf x}_j^T\beta)\nonumber \end{eqnarray*}\] Accelerated Failure Time model (AFT). \[\begin{eqnarray}\label{AFT} h_{\text{E}}^{\text{AFT}}\left(t;{\bf x}_j\right) &=& h_0\left(t \exp({\bf x}_j^T\beta)\right)\exp({\bf x}_j^T\beta),\\ H_{\text{E}}^{\text{AFT}}\left(t;{\bf x}_j\right) &=& H_0\left(t \exp({\bf x}_j^T\beta)\right).\nonumber \end{eqnarray}\] General Hazards model (GH). \[\begin{eqnarray}\label{GH} h_{\text{E}}^{\text{GH}}\left(t;{\bf x}_j\right) &=& h_0\left(t \exp({\bf x}_j^T\beta_1)\right)\exp({\bf x}_j^T\beta_2),\\ H_{\text{E}}^{\text{GH}}\left(t;{\bf x}_j\right) &=& H_0\left(t \exp({\bf x}_j^T\beta_1)\right)\exp(-{\bf x}_j^T\beta_1+{\bf x}_j^T\beta_2).\nonumber \end{eqnarray}\]

where the baseline hazard \(h_0\) is modelled using the Exponentiated Weibull distribution (see [1]). We refer the reader to [2] for more details on these models.

The following R code illustrates the use of these models in a simulated data set. The data set was simulated using the Proportional Hazards (PH) structure. The idea is to fit the parametric regression models with hazard structures PH, AH, AFT, and GH and select the one favoured by the Akaike Information Criterion (AIC). The simulated data set contains \(n=10,000\) observations with approximately \(30\%\) censoring. The user can fit subsets of this data set to evaluate the effect of sample size.

R Code and output

Importing the data, defining the useful functions (hazards, survival, likelihood) and optimisation

rm(list=ls())
# library to compile the functions
library(compiler)
enableJIT(3)
## [1] 3
library(knitr)
## Warning: package 'knitr' was built under R version 3.4.3
# True values of the parameters
ae <- 1.75; be <- 0.5; ce <- 2.5; beta <- c(0.035,0.2,0.3); rate0 <- 0; betaH = c(0,0,0)
true.par <- c(ae,be,ce,betaH,beta)
true.parc <- c(ae,be,ce,beta)

# read the data (available at: https://sites.google.com/site/fjavierrubio67/dataPH.txt?attredirects=0&d=1)
mydata <- read.table("dataPH.txt")
colnames(mydata)
## [1] "agec"      "sex"       "TTT"       "surv.time" "status"    "rate"
# create the variables for the models
x <- as.matrix(mydata[,1:3]) # covariates
sim <- mydata[,4] # survival time
ind.cens <- mydata[,5] # vital status
hp.as <- mydata[,6] # Population (expected) mortality hazard

# Summary of the data
summary(sim) # simulated survival times
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 0.000082 0.797842 2.187653 2.440057 4.225682 4.999801
apply(x,2,summary)  # covariates
##               agec    sex   TTT
## Min.    -39.991059 0.0000 0.000
## 1st Qu.  -4.999624 0.0000 0.000
## Median    2.199544 0.0000 0.000
## Mean     -1.672924 0.4999 0.498
## 3rd Qu.   8.760779 1.0000 1.000
## Max.     14.992899 1.0000 1.000
mean(ind.cens) # censoring rate
## [1] 0.6986
##############################################################################################################
# Hazard functions
##############################################################################################################

# Weibull Hazard
hw <- function(t,lambda,kappa){
  pdf0 <-  dweibull(t,scale=lambda,shape=kappa)
  cdf0 <- pweibull(t,scale=lambda,shape=kappa)
  return(pdf0/(1-cdf0))
}                                                                                      

# Weibull Cumulative hazard
Hw <- function(t,lambda,kappa){
  cdf <- pweibull(t,scale=lambda,shape=kappa)  
  return(-log(1-cdf))
}

# Exponentiated Weibull Distribution
# Hazard
hexpw <- function(t,lambda,kappa,alpha){
  pdf0 <-  alpha*dweibull(t,scale=lambda,shape=kappa)*pweibull(t,scale=lambda,shape=kappa)^(alpha-1) 
  cdf0 <- pweibull(t,scale=lambda,shape=kappa)^alpha
  cdf0 <- ifelse(cdf0==1,0.9999999,cdf0)
  return(pdf0/(1-cdf0))
}                                                                                      

# Cumulative hazard
Hexpw <- function(t,lambda,kappa,alpha,log.p=FALSE){
  cdf <- pweibull(t,scale=lambda,shape=kappa)^alpha  
  return(-log(1-cdf))
} 


###################################################################################################################
# Estimation
###################################################################################################################
# Implemented using nlminb(), alternatively, use optim()
#--------------------------------------------------------------------------------------------------------
# Proportional Hazards
#--------------------------------------------------------------------------------------------------------

#----------------------------------------------
# log-likelihood: Weibull model
#----------------------------------------------

log.likw1 <- cmpfun(function(par){
  if(anyNA(par)) return(Inf)
  ae0 <- par[1]; be0 <- par[2]; beta0 <- par[3:(3+dim(x)[2]-1)];
  if(par[1]>0 & par[2]>0 ){
    exp.x.beta <- exp(x%*%beta0)
    haz0 <- hp.as + hw(sim,ae0,be0)*exp.x.beta
    val <- - ind.cens*log(haz0) + Hw(sim,ae0,be0)*exp.x.beta
    return(sum(val))
  }
  else return(Inf)
})
# Alternatively use optim(c(ae,be,beta),log.likw1,control=list(maxit=10000))
OPTW1 <- nlminb(c(ae,be,beta),log.likw1,control=list(iter.max=10000))
MLEW1 <- OPTW1$par
AICW1 <- 2*OPTW1$objective + 2*length(MLEW1)

#----------------------------------------------
# log-likelihood: EW Model
#----------------------------------------------
log.likew1 <- cmpfun(function(par){
  if(anyNA(par)) return(Inf)
  ae0 <- par[1]; be0 <- par[2]; ce0 <- par[3]; beta0 <- par[4:(4+dim(x)[2]-1)];
  if(par[1]>0 & par[2]>0 & par[3]>0){
    exp.x.beta <- exp(x%*%beta0)
    haz0 <- hp.as + hexpw(sim,ae0,be0,ce0)*exp.x.beta
    CH0 <- Hexpw(sim,ae0,be0,ce0)*exp.x.beta
    if(anyNA(CH0)) return(Inf)
    if(anyNA(haz0)) return(Inf)
    else  val <- - ind.cens*log(haz0) + CH0
    return(sum(val))
  }
  else return(Inf)
})

OPTEW1 <- nlminb(c(ae,be,ce,beta),log.likew1,control=list(iter.max=10000))
MLEEW1 <- OPTEW1$par
AICEW1 <- 2*OPTEW1$objective + 2*length(MLEEW1)

#--------------------------------------------------------------------------------------------------------
# Accelerated Hazards
#--------------------------------------------------------------------------------------------------------

#----------------------------------------------
# log-likelihood: EW Model
#----------------------------------------------
log.likew2 <- cmpfun(function(par){
  if(anyNA(par)) return(Inf)
  ae0 <- par[1]; be0 <- par[2]; ce0 <- par[3]; beta0 <- par[4:(4+dim(x)[2]-1)];
  if(par[1]>0 & par[2]>0 & par[3]>0){
    exp.x.beta <- exp(x%*%beta0)
    haz0 <- hp.as + hexpw(sim*exp.x.beta,ae0,be0,ce0)
    CH0 <- Hexpw(sim*exp.x.beta,ae0,be0,ce0)/exp.x.beta
    if(anyNA(CH0)) return(Inf)
    if(anyNA(haz0)) return(Inf)
    else  val <- - ind.cens*log(haz0) + CH0
    return(sum(val))
  }
  else return(Inf)
})

OPTEW2 <- nlminb(c(ae,be,ce,beta),log.likew2,control=list(iter.max=10000))
MLEEW2 <- OPTEW2$par
AICEW2 <- 2*OPTEW2$objective + 2*length(MLEEW2)

#--------------------------------------------------------------------------------------------------------
# Accelerated Failure Time
#--------------------------------------------------------------------------------------------------------

#----------------------------------------------
# log-likelihood: EW Model
#----------------------------------------------
log.likew3 <- cmpfun(function(par){
  if(anyNA(par)) return(Inf)
  ae0 <- par[1]; be0 <- par[2]; ce0=par[3]; beta0 <- par[4:(4+dim(x)[2]-1)];
  if(par[1]>0 & par[2]>0 & par[3]>0){
    exp.x.beta <- exp(x%*%beta0)
    haz0 <- hp.as + hexpw(sim*exp.x.beta,ae0,be0,ce0)*exp.x.beta
    CH0 <- Hexpw(sim*exp.x.beta,ae0,be0,ce0)
    if(anyNA(CH0)) return(Inf)
    if(anyNA(haz0)) return(Inf)
    else  val <- - ind.cens*log(haz0) + CH0
    return(sum(val))
  }
  else return(Inf)
})

OPTEW3 <- nlminb(c(ae,be,ce,beta),log.likew3,control=list(iter.max=10000))
MLEEW3 <- OPTEW3$par
AICEW3 <- 2*OPTEW3$objective + 2*length(MLEEW3)

#--------------------------------------------------------------------------------------------------------
# General Model
#--------------------------------------------------------------------------------------------------------

#----------------------------------------------
# log-likelihood: EW Model
#----------------------------------------------
log.likewG <- cmpfun(function(par){
  if(anyNA(par)) return(Inf)
  ae0 <- par[1]; be0 <- par[2]; ce0=par[3]; beta0 <- par[4:(4+dim(x)[2]-1)]; beta1 <- tail(par,n = dim(x)[2]);
  if(par[1]>0 & par[2]>0 & par[3]>0){
    exp.x.beta0 <- exp(x%*%beta0)
    exp.x.beta1 <- exp(x%*%beta1)
    exp.x.dif <- exp(x%*%(beta1-beta0))
    haz0 <- hp.as + hexpw(sim*exp.x.beta0,ae0,be0,ce0)*exp.x.beta1
    CH0 <- Hexpw(sim*exp.x.beta0,ae0,be0,ce0)*exp.x.dif
    if(anyNA(CH0)) return(Inf)
    if(anyNA(haz0)) return(Inf)
    else  val <- - ind.cens*log(haz0) + CH0
    return(sum(val))
  }
  else return(Inf)
})

OPTEWG <- nlminb(c(ae,be,ce,betaH,beta), log.likewG,control=list(iter.max=10000))
MLEEWG <- OPTEWG$par
AICEWG <- 2*OPTEWG$objective + 2*length(MLEEWG)


###################################################################################################################
# Model Comparison
###################################################################################################################

AIC <- c(AICW1,AICEW1,AICEW2,AICEW3,AICEWG)
AIC.min <- which(AIC==min(AIC))
model.names <- c("W-PH","EW-PH","EW-AH","EW-AFT","EW-General")
names(AIC) <- model.names
AIC <- as.matrix(AIC);  colnames(AIC) <- c("AIC")
print(kable(AIC),digits=4)
## 
## 
##                    AIC
## -----------  ---------
## W-PH          26865.59
## EW-PH         26803.65
## EW-AH         26977.45
## EW-AFT        26816.40
## EW-General    26809.14
model.names[order(AIC)]
## [1] "EW-PH"      "EW-General" "EW-AFT"     "W-PH"       "EW-AH"
cat("Best Model based on AIC = ",model.names[AIC.min])
## Best Model based on AIC =  EW-PH
# True parameter values vs MLE of best model
print(kable(cbind(true.parc,MLEEW1),digits=4))
## 
## 
##  true.parc   MLEEW1
## ----------  -------
##      1.750   1.5221
##      0.500   0.4811
##      2.500   2.6847
##      0.035   0.0358
##      0.200   0.1965
##      0.300   0.2694

Derivation of the confidence intervals

###################################################################################################################
# Confidence intervals for best model
###################################################################################################################

# Required package
library(numDeriv)

###########################################################################################
# Function to calculate the normal confidence intervals
# The parameters indicated with "index" are transformed to the real line using log()
###########################################################################################
# FUN   : minus log-likelihood function to be used to calculate the confidence intervals
# MLE   : maximum likelihood estimator of the parameters of interest
# level : confidence level
# index : position of the positive parameters under the original parameterisation

Conf.Int <- function(FUN,MLE,level=0.95,index=NULL){
  sd.int <- abs(qnorm(0.5*(1-level)))
  tempf <- function(par){
    par[index] = exp(par[index])
    return(FUN( par ))
  }
  r.MLE <- MLE
  r.MLE[index] <- log(MLE[index])
  HESS <- hessian(tempf,x=r.MLE)
  Fisher.Info <- solve(HESS)
  Sigma <- sqrt(diag(Fisher.Info))
  U<- r.MLE + sd.int*Sigma
  L<- r.MLE - sd.int*Sigma
  C.I <- cbind(L,U,r.MLE, Sigma)
  names.row <- paste0("par", seq_along(1:length(MLE)))
  names.row[index] <- paste0("log.par", seq_along(index))
  rownames(C.I)<- names.row
  colnames(C.I)<- c("Lower","Upper","Transf MLE", "Std. Error")
  return(C.I)
}

CI <- Conf.Int(log.likew1,MLEEW1,level=0.95,index=c(1,2,3))
print(kable(CI,digits=4))
## 
## 
##               Lower     Upper   Transf MLE   Std. Error
## ---------  --------  --------  -----------  -----------
## log.par1    -0.0270    0.8671       0.4201       0.2281
## log.par2    -0.8843   -0.5790      -0.7317       0.0779
## log.par3     0.7248    1.2503       0.9876       0.1341
## par4         0.0334    0.0381       0.0358       0.0012
## par5         0.1406    0.2524       0.1965       0.0285
## par6         0.2133    0.3255       0.2694       0.0286

Results: Graphical comparisons and Net Survival estimates derived from the selected model

###################################################################################################################
# Fitted hazard vs True hazard
###################################################################################################################


hazard.fit <- Vectorize(function(t)  hexpw(t,MLEEW1[1],MLEEW1[2],MLEEW1[3]))
hazard.true <- Vectorize(function(t)  hexpw(t,ae,be,ce))

curve(hazard.true,0,5,lwd=3,n=10000,xlab="time",ylab="Excess Hazard",cex.axis=1.5,cex.lab=1.5,ylim=c(0,0.26), col="black")
curve(hazard.fit,0,5,lwd=3,n=10000,lty=2,col="black",add=T)
legend(3.5,0.24, c("True","Fitted"),
       text.col = c("black","black"), col = c("black","black"), lty = c(1, 2),lwd=3,
       merge = TRUE, bg = "gray90")

###################################################################################################################
# Fitted Net survival for a single patient with sex = 0, age = 70, Comorbidity = 0, 1 
###################################################################################################################

NS0 <- Vectorize(function(t)  exp(-Hexpw(t,MLEEW1[1],MLEEW1[2],MLEEW1[3])))
NS1 <- Vectorize(function(t)  exp(-Hexpw(t,MLEEW1[1],MLEEW1[2],MLEEW1[3])*exp(MLEEW1[6])))

curve(NS0,0,5,lwd=3,n=10000,xlab="time",ylab="Net Survival",cex.axis=1.5,cex.lab=1.5,ylim=c(0,1), col="blue")
curve(NS1,0,5,lwd=3,n=10000,lty=2,col="red",add=T)
legend(3.5,0.9, c("Comorb = 0","Comorb = 1"),
       text.col = c("blue","red"), col = c("blue","red"), lty = c(1, 2),lwd=3,
       merge = TRUE, bg = "gray90")

########################################################################################################
# Confidence intervals for the net survival based on Bootstrap
########################################################################################################
# Required package
library(mvtnorm)

# t0 : time at what the confidence interval will be calculated
# level : confidence level
# n.boot : number of bootstrap iterations
# x.int : subset of the design matrix associated to the group of interest

# The variance of the normal approximation to the distribution of the MLE
 tempf <- function(par){
    par[1:3] = exp(par[1:3])
    return(log.likew1( par ))
  }
  r.MLE <- MLEEW1
  r.MLE[1:3] <- log(MLEEW1[1:3])
  HESS <- hessian(tempf,x=r.MLE)
  Sigma <- solve(HESS)

  # The function 
conf.int.NS <- function(t0,level,n.boot,x.int){
  boot <- vector()
  S.par <- function(par) mean( exp(-Hexpw(t0,par[1],par[2],par[3])*exp(x.int%*%par[-c(1:3)])))

  for(i in 1:n.boot) {
    val <- rmvnorm(1,mean = r.MLE, sigma = Sigma)
    val[1:3] <- exp(val[1:3])
    boot[i] <- S.par(val)
  }
  
  L <- quantile(boot,(1-level)*0.5)
  U <- quantile(boot,(1+level)*0.5)
  
  M <- S.par(MLEEW1)
  
  return(c(L,M,U))
}

#-----------------------------------
# Net Survival for the whole Cohort
# at Comorbidity = 0 , 1 
# at years 1, 2, 3, 4, 4.9
#-----------------------------------

# times
times <- c(1,2,3,4,4.9)

CIS0 <- CIS1 <- matrix(0, ncol = 4, nrow = length(times))

for(k in 1:length(times)) CIS0[k,] <- c(times[k],conf.int.NS(times[k],0.95,10000,x[which(x[,3]==0),]))
for(k in 1:length(times)) CIS1[k,] <- c(times[k],conf.int.NS(times[k],0.95,10000,x[which(x[,3]==1),]))

colnames(CIS0) <- cbind("year","lower","net_survival","upper")
print(kable(CIS0,digits=4))
## 
## 
##  year    lower   net_survival    upper
## -----  -------  -------------  -------
##   1.0   0.7575         0.7660   0.7786
##   2.0   0.6020         0.6131   0.6302
##   3.0   0.4958         0.5079   0.5266
##   4.0   0.4185         0.4311   0.4507
##   4.9   0.3652         0.3779   0.3980
colnames(CIS1) <- cbind("year","lower","net_survival","upper")
print(kable(CIS1,digits=4))
## 
## 
##  year    lower   net_survival    upper
## -----  -------  -------------  -------
##   1.0   0.6992         0.7084   0.7229
##   2.0   0.5223         0.5332   0.5509
##   3.0   0.4094         0.4206   0.4398
##   4.0   0.3319         0.3430   0.3614
##   4.9   0.2806         0.2916   0.3106
#--------------------------------------------
# Net Survival for an age group: 60 - 70
# at Comorbidity = 0 , 1 
# sex = 1
#--------------------------------------------

CISA0 <- CISA1 <- matrix(0, ncol = 4, nrow = length(times))

for(k in 1:length(times)) CISA0[k,] <- c(times[k],conf.int.NS(times[k],0.95,10000,x[which(x[,3]==0 & x[,1]>=-10 & x[,1]<=0 & x[,2]==1),]))
for(k in 1:length(times)) CISA1[k,] <- c(times[k],conf.int.NS(times[k],0.95,10000,x[which(x[,3]==1 & x[,1]>=-10 & x[,1]<=0 & x[,2]==1),]))

colnames(CISA0) <- cbind("year","lower","net_survival","upper")
print(kable(CISA0,digits=4))
## 
## 
##  year    lower   net_survival    upper
## -----  -------  -------------  -------
##   1.0   0.7660         0.7764   0.7903
##   2.0   0.6086         0.6228   0.6425
##   3.0   0.4975         0.5130   0.5349
##   4.0   0.4144         0.4306   0.4544
##   4.9   0.3556         0.3723   0.3973
colnames(CISA1) <- cbind("year","lower","net_survival","upper")
print(kable(CISA1,digits=4))
## 
## 
##  year    lower   net_survival    upper
## -----  -------  -------------  -------
##   1.0   0.7046         0.7163   0.7315
##   2.0   0.5214         0.5357   0.5568
##   3.0   0.3995         0.4149   0.4373
##   4.0   0.3140         0.3294   0.3529
##   4.9   0.2564         0.2721   0.2951