Around cure models

Setup

Method: Mixture cure models with age and year as continuous covariables

In the following models, we assumed a non-linear effect of age at diagnosis both on the survival probability for the uncured and on the probability of cure. However, for the year of diagnosis, we fitted 4 models, which vary depending on linear or non-linear effect assumed for the year of diagnosis on either the survival probability of the uncured or the probability of cure.

We fit the following models, with S(t;x)=π(x)+(1π(x))Su(t;x), where π(x) is the probability of cure and Su(t;x) is the survival of the uncured.

We assumed a Weibull distribution for the uncured, Su(t;x)=exp(λ(x)tκ) , where λ is the scale and κ is the shape parameter, and we assumed a logit parametrisation for the cure probability, π(x)=exp(xβ)1+exp(xβ). The matrix x represents the age at diagnosis and the year of diagnosis.

Mod.Fixed.mixcure_2X <- optim(par=rep(0.1, 15), fn=Tot.LogLik.mixcure_2X, hessian=T, method = "BFGS",
                              Time=mydat$timesurv15y, Status = mydat$vstat15y, varyear=mydat$ydiagcr, 
                              varage=mydat$agediagcr, control=list(maxit=10000))
Mod.Fixed.mixcure_2X$convergence
## [1] 0
Mod.Fixed.mixcure_2X.SuYLin <- optim(par=rep(0.1, 13), fn=Tot.LogLik.mixcure_2X, hessian=T, method = "BFGS",
                              Time=mydat$timesurv15y, Status = mydat$vstat15y, varyear=mydat$ydiagcr, 
                              varage=mydat$agediagcr, control=list(maxit=10000), ModForm="SuYLin")
Mod.Fixed.mixcure_2X.SuYLin$convergence
## [1] 0
# Akaike Information Criteria

print("Both Non linear:")
## [1] "Both Non linear:"
2*Mod.Fixed.mixcure_2X$value+2*length(Mod.Fixed.mixcure_2X$par)
## [1] 1912.64
print("Su with linear effect of year")
## [1] "Su with linear effect of year"
2*Mod.Fixed.mixcure_2X.SuYLin$value+2*length(Mod.Fixed.mixcure_2X.SuYLin$par)
## [1] 1910.75
ModFinal <- Mod.Fixed.mixcure_2X.SuYLin

Results

expit <- function(x){exp(x)/(1+exp(x))}
theagediag=seq(0,24,by=0.1); myagediagcr=(theagediag-agediagk)/10 # corresponds to mydat$agediagc/10
linpredage <- ModFinal$par[c(2,6:8)]%*%rbind(1,myagediagcr, myagediagcr^2, (myagediagcr^2)*(myagediagcr>0))
myt <- seq(0,15,by=0.1)

kap     <- ModFinal$par[1]
alpha0   <- ModFinal$par[2]
alpha1   <- ModFinal$par[3]
alpha2   <- ModFinal$par[4]
alpha3   <- ModFinal$par[5]
alpha4   <- ModFinal$par[6]
alpha5   <- ModFinal$par[7]
alpha6   <- ModFinal$par[8]
beta0    <- ModFinal$par[9]
beta1    <- ModFinal$par[10]
beta4    <- ModFinal$par[11]
beta5    <- ModFinal$par[12]
beta6    <- ModFinal$par[13]

myvarcov <- solve(ModFinal$hessian)
stdErr <- sqrt(diag(myvarcov))

Related