Disclaimer

The artificial data provided here and the numerical results represent numerical illustrations of the methods proposed in Rubio, Putter, and Belot (2023). These data and results should not be used to make epidemiological claims.

The excess hazard regression model

We refer the reader to Rubio et al. (2019) and Eletti et al. (2022) for a review of excess hazard regression models and the assumptions behind the relative survival framework. The model is formulated in terms of the excess hazard function multiplied by an individual frailty: \[ \begin{eqnarray} \lambda h_E(t; {\bf x}, \alpha, \beta, \theta) &=& \lambda h_0(t \exp\{ {\tilde{\bf x}}^{\top}\alpha\}; \theta) \exp\{{\bf x}^{\top}\beta\}, \\ \lambda &\sim & G, \end{eqnarray} \] where \(h_0(;\theta)\) is a parametric baseline hazard. This model contains the Proportional hazards (\(\alpha=0\)), accelerated hazards (\(\beta=0\)), and accelerated failure time (\(\alpha = \beta\), \(\tilde{\bf x} = {\bf x}\)) models as particular cases.

We compare this model with the Proportional Hazards (PH) and Accelerated Failure Time (AFT) models.

See also: Rubio et al. (2019) and Eletti et al. (2022).

The simulated data set

The data set used here is based on the Simulacrum data set, which contains artificial patient-like cancer data to help researchers gain insights. We consider the cohort of female lung cancer patients diagnosed in 2013, with follow-up until 2019. Background mortality rates were obtained for each cancer patient from population life tables for England defined for each calendar year in 2010-2015, and stratified by single year, age, sex, and deprivation category.

We consider the time dependent effect of the variable age (\(\tilde{\bf x}\), scaled), and the hazard-level effects (\({\bf x}\)) of age (scaled) and deprivation (categorical I-V). We explore the effect of missing the variable deprivation.

See also: GHSurv.

Data Preparation

rm(list=ls())

############################################################################################################
# Data preparation
############################################################################################################

# Required packages
library(survival)
library(knitr)
library(mvtnorm)
library(numDeriv)
library(knitr)

#library(devtools)
#install_github("FJRubio67/IFNS")
library(IFNS)

# Load the simulated data set
load('data_lung.Rda')

dim(df)
## [1] 16828    14
# A snipet of the data frame
head(df)
##          time year dep dep.1 dep.2 dep.3 dep.4 dep.5 age  age.out sex status
## 1 5.776865000 2013   2     0     1     0     0     0  18 23.77687   1      0
## 2 0.084873371 2013   2     0     1     0     0     0  20 20.08487   1      1
## 3 0.388774810 2013   1     1     0     0     0     0  20 20.38877   1      1
## 4 5.771389500 2013   4     0     0     0     1     0  20 25.77139   1      0
## 5 0.002737851 2013   5     0     0     0     0     1  21 21.00274   1      1
## 6 5.264886900 2013   2     0     1     0     0     0  22 27.26489   1      0
##          hrate      agec
## 1 0.0001792260 -5.094741
## 2 0.0001898829 -4.912291
## 3 0.0001898829 -4.912291
## 4 0.0001898829 -4.912291
## 5 0.0001944859 -4.821066
## 6 0.0002008885 -4.729842
# Summaries

# Deprivation level
table(df$dep)
## 
##    1    2    3    4    5 
## 2340 3043 3431 3747 4267
# Age at diagnosis
summary(df$age)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   18.00   67.00   74.00   73.85   82.00   99.00
# Proportion of observed times to event
mean(df$status)
## [1] 0.7747801
#******************
# Data preparation
#******************

n.ind <- dim(df)[1]  # number of individuals
# Design matrices
X <- as.matrix(cbind( df$agec, df$dep.2, df$dep.3, df$dep.4, df$dep.5 ))
colnames(X) <- cbind("std age", "dep.2", "dep.3", "dep.4", "dep.5")
Xt <- as.matrix(df$agec)
q <- dim(Xt)[2]
p <- dim(X)[2]

# Histogram of survival times
hist(df$time, breaks = 50, probability = T, 
     xlab = "Time (years)", ylab = "Density", main = "Simulacrum Lung Cancer (F)")
box()

# Required quantities
status <- as.vector(df$status)
times <- as.vector(df$time) # in years
n.obs <- sum(status)
# Population hazard rates for all individuals
hp.rates <- df$hrate
# Population hazard rates for the uncensored individuals
hp.obs <- df$hrate[as.logical(status)]

Model fitting using the IFNS package

The IFNS R package implements several excess hazard models with and without frailties . The R command PGWMLE can be used to fit Excess Hazard Models with a PGW baseline and several hazard structures: PGW without covariates (“PGW”), PGW with AFT model (PGWAFT) PGW with PH model (PGWPH) PGW with AH model (PGWAH) PGW with GH model (PGWGH). The R command IFPGWMLE can be used to fit frailty Excess Hazard Models with a PGW baseline and several hazard structures: PGW without covariates (“IFPGW”), PGW with AFT model (IFPGWAFT) PGW with PH model (IFPGWPH) PGW with AH model (IFPGWAH) PGW with GH model (IFPGWGH). The corresponding models with LogNormal baseline hazard can be fitted using the command LNMLE and IFLNMLE.

Models with all covariates

The following R code shows how to fit the classical (no frailty) and the frailty excess hazards models including all covariates available.

# Initial points
(modinit <- PGWMLE(init=c(0,0,2,rep(0,ncol(X))), hstr = "PGWPH", times = times, 
                   status = status, rates = hp.rates, des = X, maxit = 10000))
## $OPT
## $OPT$par
## [1] -1.73089387  0.02579858  1.90280331 -0.21810282  0.69136786  0.65910209
## [7]  0.76428032  0.80343345
## 
## $OPT$value
## [1] 19616.31
## 
## $OPT$counts
## function gradient 
##     1915       NA 
## 
## $OPT$convergence
## [1] 0
## 
## $OPT$message
## NULL
## 
## 
## $loglik
## function (par) 
## {
##     ae0 <- exp(par[1])
##     be0 <- exp(par[2])
##     ce0 <- exp(par[3])
##     beta <- par[4:(3 + p)]
##     exp.x.beta <- as.vector(exp(des %*% beta))
##     exp.x.beta.obs <- exp.x.beta[status]
##     haz0 <- hp.as.obs + hpgw(times.obs, ae0, be0, ce0) * exp.x.beta.obs
##     val <- -sum(log(haz0)) + sum(chpgw(times, ae0, be0, ce0) * 
##         exp.x.beta)
##     return(val)
## }
## <bytecode: 0x1412b2238>
## <environment: 0x1412b9860>
vect.c0 <- c(modinit$OPT$par[1:3], rep(0,ncol(Xt)), modinit$OPT$par[-c(1:3)])


#################################################################################
# MODEL FITTING
#################################################################################

#################################################################################
# All covariates
#################################################################################


#--------------------------------------------------------------------------------
# Classical model
#--------------------------------------------------------------------------------

OPTC01 <- try(PGWMLE(vect.c0,
                     times = times, status = status, rates = hp.rates,
                     hstr = "PGWGH", des_t = Xt, des = X, method = "Nelder-Mead", maxit = 10000), silent=T)
if (class(OPTC01)=="try-error") OPTC01$value<-Inf
OPTC02 <- try(PGWMLE(vect.c0,
                     times = times, status = status, rates = hp.rates,
                     hstr = "PGWGH", des_t = Xt, des = X, method = "nlminb", maxit = 10000), silent=T)
if (class(OPTC02)=="try-error") OPTC02$objective<-Inf
OPTC03 <- try(PGWMLE(vect.c0*0,
                     times = times, status = status, rates = hp.rates,
                     hstr = "PGWGH", des_t = Xt, des = X, method = "Nelder-Mead", maxit = 10000), silent=T)
if (class(OPTC03)=="try-error") OPTC03$value<-Inf
OPTC04 <- try(PGWMLE(vect.c0*0,
                     times = times, status = status, rates = hp.rates,
                     hstr = "PGWGH", des_t = Xt, des = X, method = "nlminb", maxit = 10000), silent=T)
if (class(OPTC04)=="try-error") OPTC04$objective<-Inf


idc0 <- which.min(c(OPTC01$OPT$value, OPTC02$OPT$objective, OPTC03$OPT$value, OPTC04$OPT$objective))
if(idc0 == 1) OPTC0 = OPTC01
if(idc0 == 2) OPTC0 = OPTC02
if(idc0 == 3) OPTC0 = OPTC03
if(idc0 == 4) OPTC0 = OPTC04
(MLEC0<- c(exp(OPTC0$OPT$par[1:3]),OPTC0$OPT$par[-c(1:3)]))
## [1]  0.10073483  1.15936186  6.01470158  0.38982360 -0.03574526  0.16875930
## [7]  0.11142544  0.14610973  0.22164112
loglikc0 <- c(OPTC01$OPT$value, OPTC02$OPT$objective, OPTC03$OPT$value, OPTC04$OPT$objective)[idc0]


#--------------------------------------------------------------------------------
# Frailty model
#--------------------------------------------------------------------------------
vect.f0 <- c(OPTC0$OPT$par,0)

OPTF01 <- try(IFPGWMLE(vect.f0,
                       times = times, status = status, rates = hp.rates,
                       hstr = "IFPGWGH", des_t = Xt, des = X, method = "Nelder-Mead", maxit = 10000), silent=T)
if (class(OPTF01)=="try-error") OPTF01$value<-Inf
OPTF02 <- try(IFPGWMLE(vect.f0,
                       times = times, status = status, rates = hp.rates,
                       hstr = "IFPGWGH", des_t = Xt, des = X, method = "nlminb", maxit = 10000), silent=T)
if (class(OPTF02)=="try-error") OPTF02$value<-Inf
OPTF03 <- try(IFPGWMLE(vect.f0*0,
                       times = times, status = status, rates = hp.rates,
                       hstr = "IFPGWGH", des_t = Xt, des = X, method = "Nelder-Mead", maxit = 10000), silent=T)
if (class(OPTF03)=="try-error") OPTF03$value<-Inf
OPTF04 <- try(IFPGWMLE(vect.f0*0,
                       times = times, status = status, rates = hp.rates,
                       hstr = "IFPGWGH", des_t = Xt, des = X, method = "nlminb", maxit = 10000), silent=T)
if (class(OPTF04)=="try-error") OPTF04$objective<-Inf

idf0 <- which.min(c(OPTF01$OPT$value, OPTF02$OPT$objective, OPTF03$OPT$value, OPTF04$OPT$objective))
if(idf0 == 1) OPTF0 = OPTF01
if(idf0 == 2) OPTF0 = OPTF02
if(idf0 == 3) OPTF0 = OPTF03
if(idf0 == 4) OPTF0 = OPTF04
(MLEF0 <- c(exp(OPTF0$OPT$par[1:3]),OPTF0$OPT$par[-c(1:3,length(OPTF0$OPT$par))],
            exp(OPTF0$OPT$par[length(OPTF0$OPT$par)])))
##  [1]  0.3145450  1.0789340  2.1679828  1.5640598 -0.1161342  0.2409230
##  [7]  0.1762914  0.2040043  0.3193708  1.3057667
loglikf0 <- c(OPTF01$OPT$value, OPTF02$OPT$objective, OPTF03$OPT$value, OPTF04$OPT$objective)[idf0]

Models after removing deprivation

#################################################################################
# NO deprivation
#################################################################################

# Design matrix without stage
X1 <- as.matrix(X[,-which(colnames(X) %in% c("dep.2","dep.3","dep.4","dep.5" ))])

# Initial points
(modinit <- PGWMLE(init=c(0,0,2,rep(0,ncol(X1))), hstr = "PGWPH", times = times, 
                   status = status, rates = hp.rates, des = X1, maxit = 10000))
## $OPT
## $OPT$par
## [1] -2.3398345  0.1724236  1.7319940 -0.2445874
## 
## $OPT$value
## [1] 19434.84
## 
## $OPT$counts
## function gradient 
##      339       NA 
## 
## $OPT$convergence
## [1] 0
## 
## $OPT$message
## NULL
## 
## 
## $loglik
## function (par) 
## {
##     ae0 <- exp(par[1])
##     be0 <- exp(par[2])
##     ce0 <- exp(par[3])
##     beta <- par[4:(3 + p)]
##     exp.x.beta <- as.vector(exp(des %*% beta))
##     exp.x.beta.obs <- exp.x.beta[status]
##     haz0 <- hp.as.obs + hpgw(times.obs, ae0, be0, ce0) * exp.x.beta.obs
##     val <- -sum(log(haz0)) + sum(chpgw(times, ae0, be0, ce0) * 
##         exp.x.beta)
##     return(val)
## }
## <bytecode: 0x1412b2238>
## <environment: 0x1424c6240>
vect.c1 <- c(modinit$OPT$par[1:3], rep(0,ncol(Xt)), modinit$OPT$par[-c(1:3)])

#--------------------------------------------------------------------------------
# Classical model
#--------------------------------------------------------------------------------

OPTC11 <- try(PGWMLE(vect.c1,
                     times = times, status = status, rates = hp.rates,
                     hstr = "PGWGH", des_t = Xt, des = X1, method = "Nelder-Mead", maxit = 10000), silent=T)
if (class(OPTC11)=="try-error") OPTC11$value<-Inf
OPTC12 <- try(PGWMLE(vect.c1,
                     times = times, status = status, rates = hp.rates,
                     hstr = "PGWGH", des_t = Xt, des = X1, method = "nlminb", maxit = 10000), silent=T)
if (class(OPTC12)=="try-error") OPTC12$objective<-Inf
OPTC13 <- try(PGWMLE(vect.c1*0,
                     times = times, status = status, rates = hp.rates,
                     hstr = "PGWGH", des_t = Xt, des = X1, method = "Nelder-Mead", maxit = 10000), silent=T)
if (class(OPTC13)=="try-error") OPTC13$value<-Inf
OPTC14 <- try(PGWMLE(vect.c1*0,
                     times = times, status = status, rates = hp.rates,
                     hstr = "PGWGH", des_t = Xt, des = X1, method = "nlminb", maxit = 10000), silent=T)
if (class(OPTC14)=="try-error") OPTC14$objective<-Inf


idc1 <- which.min(c(OPTC11$OPT$value, OPTC12$OPT$objective, OPTC13$OPT$value, OPTC14$OPT$objective))
if(idc1 == 1) OPTC1 = OPTC11
if(idc1 == 2) OPTC1 = OPTC12
if(idc1 == 3) OPTC1 = OPTC13
if(idc1 == 4) OPTC1 = OPTC14
(MLEC1 <- c(exp(OPTC1$OPT$par[1:3]),OPTC1$OPT$par[-c(1:3)]))
## [1]  0.09246080  1.16991677  5.61848823  0.38947868 -0.03740248
loglikc1 <- c(OPTC11$OPT$value, OPTC12$OPT$objective, OPTC13$OPT$value, OPTC14$OPT$objective)[idc1]


#--------------------------------------------------------------------------------
# Frailty model
#--------------------------------------------------------------------------------
vect.init <- c(OPTC1$OPT$par,0)

OPTF11 <- try(IFPGWMLE(vect.init,
                       times = times, status = status, rates = hp.rates,
                       hstr = "IFPGWGH", des_t = Xt, des = X1, method = "Nelder-Mead", maxit = 10000), silent=T)
if (class(OPTF11)=="try-error") OPTF11$value<-Inf
OPTF12 <- try(IFPGWMLE(vect.init,
                       times = times, status = status, rates = hp.rates,
                       hstr = "IFPGWGH", des_t = Xt, des = X1, method = "nlminb", maxit = 10000), silent=T)
if (class(OPTF12)=="try-error") OPTF12$value<-Inf
OPTF13 <- try(IFPGWMLE(vect.init*0,
                       times = times, status = status, rates = hp.rates,
                       hstr = "IFPGWGH", des_t = Xt, des = X1, method = "Nelder-Mead", maxit = 10000), silent=T)
if (class(OPTF13)=="try-error") OPTF13$value<-Inf
OPTF14 <- try(IFPGWMLE(vect.init*0,
                       times = times, status = status, rates = hp.rates,
                       hstr = "IFPGWGH", des_t = Xt, des = X1, method = "nlminb", maxit = 10000), silent=T)
if (class(OPTF14)=="try-error") OPTF14$objective<-Inf

idf1 <- which.min(c(OPTF11$OPT$value, OPTF12$OPT$objective, OPTF13$OPT$value, OPTF14$OPT$objective))
if(idf1 == 1) OPTF1 = OPTF11
if(idf1 == 2) OPTF1 = OPTF12
if(idf1 == 3) OPTF1 = OPTF13
if(idf1 == 4) OPTF1 = OPTF14
(MLEF1 <- c(exp(OPTF1$OPT$par[1:3]),OPTF1$OPT$par[-c(1:3,length(OPTF1$OPT$par))],
            exp(OPTF1$OPT$par[length(OPTF1$OPT$par)])))
## [1]  0.2669244  1.0843595  2.0530040  1.6020597 -0.1176661  1.3227951
loglikf1<- c(OPTF11$OPT$value, OPTF12$OPT$objective, OPTF13$OPT$value, OPTF14$OPT$objective)[idf1]

Models with no covariates

#################################################################################
# No covariates
#################################################################################

#--------------------------------------------------------------------------------
# Classical model
#--------------------------------------------------------------------------------

OPTC3 <- PGWMLE(init=c(0,0,2), hstr = "PGW", times = times, method = "nlminb", 
                status = status, rates = hp.rates, maxit = 10000)
MLEC3 <- c(exp(OPTC3$OPT$par[1:3]))
loglikc3 <- OPTC3$OPT$objective


#--------------------------------------------------------------------------------
# Frailty model
#--------------------------------------------------------------------------------
vect.f3 <- c(OPTC3$OPT$par,2)

OPTF3 <- IFPGWMLE(vect.f3,
                  times = times, status = status, rates = hp.rates,
                  hstr = "IFPGW", method = "nlminb", maxit = 10000)

MLEF3 <- c(exp(OPTF3$OPT$par[1:4]))
loglikf3 <- OPTF3$OPT$objective

Model comparison

# Model comparison using AIC
AICC0 <- 2*loglikc0 + 2*length(MLEC0)
AICF0 <- 2*loglikf0 + 2*length(MLEF0)
AICC1 <- 2*loglikc1 + 2*length(MLEC1)
AICF1 <- 2*loglikf1 + 2*length(MLEF1)
AICC3 <- 2*loglikc3 + 2*length(MLEC3)
AICF3 <- 2*loglikf3 + 2*length(MLEF3)

AICS <- c(AICC0,AICF0,AICC1,AICF1,AICC3,AICF3)
which.min(AICS)
## [1] 2
# Using AIC, we would select the model with all covariates
# followed by the model with age, dep, and comorbidities
AICC <- c(AICC0,AICC1,AICC3)
which.min(AICC)
## [1] 1
order(AICC)
## [1] 1 2 3
# Using AIC, we would select the model with all covariates
# followed by the model with age, dep, and comorbidities
AICF <- c(AICF0,AICF1,AICF3)
which.min(AICF)
## [1] 1
order(AICF)
## [1] 1 2 3
# MLEs classical models
MLECs <- cbind(MLEC0,c(MLEC1,rep(NA,4)),c(MLEC3,rep(NA,6)))
colnames(MLECs) <- c("All vars", "No dep", "No vars")
kable(MLECs, digits = 3)
All vars No dep No vars
0.101 0.092 0.091
1.159 1.170 1.198
6.015 5.618 5.711
0.390 0.389 NA
-0.036 -0.037 NA
0.169 NA NA
0.111 NA NA
0.146 NA NA
0.222 NA NA 
# MLEs frailty models

MLEFs <- cbind(MLEF0,c(MLEF1[1:5],rep(NA,4),MLEF1[6]),c(MLEF3[1:3],rep(NA,6),MLEF3[4]))
colnames(MLEFs) <- c("All vars", "No dep", "No vars")
kable(MLEFs, digits = 3)
All vars No dep No vars
0.315 0.267 0.349
1.079 1.084 1.085
2.168 2.053 1.647
1.564 1.602 NA
-0.116 -0.118 NA
0.241 NA NA
0.176 NA NA
0.204 NA NA
0.319 NA NA
1.306 1.323 1.741

Net survival curves: whole population

###############################################################################################
#  Fitted net survival
###############################################################################################

# Classical models
NS.fitc0 <- Vectorize(function(t) NSfit_PGWGH(t,MLEC0,Xt,X))
NS.fitc1 <- Vectorize(function(t) NSfit_PGWGH(t,MLEC1,Xt,X1))
NS.fitc3 <- Vectorize(function(t) spgw(t,MLEC3[1],MLEC3[2],MLEC3[3]))

# Frailty models
NS.fitf0 <- Vectorize(function(t) NSfit_IFPGWGH(t,MLEF0,Xt,X))
NS.fitf1 <- Vectorize(function(t) NSfit_IFPGWGH(t,MLEF1,Xt,X1))
NS.fitf3 <- Vectorize(function(t) 1/(1+MLEF3[4]*chpgw(t,MLEF3[1],MLEF3[2],MLEF3[3]))^(1/MLEF3[4]) )

curve(NS.fitc0,0,4,lwd=3,n=100,xlab="time",ylab="Net Survival",cex.axis=1.5,cex.lab=1.5,ylim=c(0,1), col="gray")
curve(NS.fitf0,0,4,lwd=3,n=100,lty=1,col="black",add=T)
curve(NS.fitc1,0,4,lwd=3,n=100,lty=2,col="gray",add=T)
curve(NS.fitf1,0,4,lwd=3,n=100,lty=2,col="black",add=T)
curve(NS.fitc3,0,4,lwd=3,n=100,lty=3,col="gray",add=T)
curve(NS.fitf3,0,4,lwd=3,n=100,lty=3,col="black",add=T)
legend("topright", c("I","II", "III","IV","V","VI"),
       text.col = c("black","black","black","black","black","black"), 
       col = c("gray","black","gray","black","gray","black"), 
       lty = c(1,1,2,2,3,3),lwd=3, merge = TRUE, bg = "white")

Net survival curves: dep-I vs dep-V

#########################################################################
# Comparison of net survival for Deprivation I-V
#########################################################################

indI <- which(X[,2]==0 & X[,3]==0 & X[,4] == 0 & X[,5] == 0)
indV <- which(X[,5]==1)

# Classical model
NS.fitc0I <- Vectorize(function(t) NSfit_PGWGH(t,MLEC0,Xt[indI],X[indI,]))
NS.fitc0V <- Vectorize(function(t) NSfit_PGWGH(t,MLEC0,Xt[indV],X[indV,]))

# Frailty model
NS.fitf0I <- Vectorize(function(t) NSfit_IFPGWGH(t,MLEF0,Xt[indI],X[indI,]))
NS.fitf0V <- Vectorize(function(t) NSfit_IFPGWGH(t,MLEF0,Xt[indV],X[indV,]))


curve(NS.fitc0I,0,4,lwd=3,n=100,xlab="time",ylab="Net Survival",cex.axis=1.5,cex.lab=1.5,ylim=c(0,1), col="gray")
curve(NS.fitc0V,0,4,lwd=3,n=100,lty=2,col="gray",add=T)
curve(NS.fitf0I,0,4,lwd=3,n=100,lty=1,col="black",add=T)
curve(NS.fitf0V,0,4,lwd=3,n=100,lty=2,col="black",add=T)
legend("topright", c("Dep I", "Dep V"),
       text.col = c("black","black"), 
       col = c("black","black"), 
       lty = c(1,2),lwd=3, merge = TRUE, bg = "white")

Confidence intervals

######################################################
# Confidence intervals for Net survival at time t0s
######################################################

# dep I: Classical vs Frailty
t0s <- c(1,2,3,4)

C.depI <- F.depI <- matrix(0, ncol = 3, nrow = length(t0s))

for(i in 1:length(t0s)){
        set.seed(i)
        C.depI[i,] <- CI.NSfit_PGWGH(t0 = t0s[i], NMC = 1000, MLE = MLEC0, level = 0.95, times = times, status = status, rates = hp.rates, 
                                       des_t = Xt,des = X, subgroup = indI)
        
        F.depI[i,] <- CI.NSfit_IFPGWGH(t0 = t0s[i], NMC = 1000, MLE = MLEF0, level = 0.95, times = times, status = status, rates = hp.rates, 
                                         des_t = Xt,des = X, subgroup = indI)
}


M.depI <- cbind(C.depI, F.depI)
colnames(M.depI) <- c("C-lower","C-estimate","C-upper","F-lower","F-estimate","F-upper")
rownames(M.depI) <- c("t=1","t=2","t=3","t=4")

kable(M.depI, digits = 3)
C-lower C-estimate C-upper F-lower F-estimate F-upper
t=1 0.547 0.566 0.584 0.541 0.557 0.573
t=2 0.441 0.458 0.476 0.428 0.445 0.462
t=3 0.381 0.400 0.418 0.371 0.388 0.405
t=4 0.341 0.361 0.378 0.335 0.351 0.367 
# dep V: Classical vs Frailty
t0s <- c(1,2,3,4)

C.depV <- F.depV <- matrix(0, ncol = 3, nrow = length(t0s))

for(i in 1:length(t0s)){
        set.seed(i)
        C.depV[i,] <- CI.NSfit_PGWGH(t0 = t0s[i], NMC = 1000, MLE = MLEC0, level = 0.95, times = times, status = status, rates = hp.rates, 
                                        des_t = Xt,des = X, subgroup = indV)
        
        F.depV[i,] <- CI.NSfit_IFPGWGH(t0 = t0s[i], NMC = 1000, MLE = MLEF0, level = 0.95, times = times, status = status, rates = hp.rates, 
                                          des_t = Xt,des = X, subgroup = indV)
}

M.depV <- cbind(C.depV, F.depV)
colnames(M.depV) <- c("C-lower","C-estimate","C-upper","F-lower","F-estimate","F-upper")
rownames(M.depV) <- c("t=1","t=2","t=3","t=4")

kable(M.depV, digits = 3)
C-lower C-estimate C-upper F-lower F-estimate F-upper
t=1 0.474 0.488 0.501 0.468 0.481 0.494
t=2 0.362 0.375 0.388 0.360 0.373 0.386
t=3 0.302 0.316 0.330 0.309 0.320 0.332
t=4 0.265 0.278 0.290 0.277 0.287 0.299

References

Eletti, A., G. Marra, M. Quaresma, R. Radice, and F. J. Rubio. 2022. “A Unifying Framework for Flexible Excess Hazard Modelling with Applications in Cancer Epidemiology.” Journal of the Royal Statistical Society: Series C 71: 1044–62.
Rubio, F. J., H. Putter, and A. Belot. 2023. “Individual Frailty Excess Hazard Models in Cancer Epidemiology.” Statistics in Medicine 42: 1066–81.
Rubio, F. J., L. Remontet, N. P. Jewell, and A. Belot. 2019. “On a General Structure for Hazard-Based Regression Models: An Application to Population-Based Cancer Research.” Statistical Methods in Medical Research 28: 2404–17.