R/AFTeweibull.R

In sa4khan/AFTjmr: Recurrent event analysis and joint modeling using accelerated failure time (AFT) models

###################################################
# Negative log-likelihood of Exp Weibull AFT
####################################################
ewaft.llik<-function(init,st,status,xx,precBits=NULL){
p<-ncol(xx)
beta<-init[1:p]
kappa<-exp(init[p+1])
gam<-exp(init[p+2])
rho<-exp(init[p+3])
pred<-as.matrix(xx)%*%beta
u<-rho*st*exp(-pred)
if(is.null(precBits)){
r1<-log1p(-exp(gam*log1p(-exp(-u^kappa))))
s1<-log1p(-exp(-u^kappa))
} else{
r1<-as.numeric(log1p(-exp(gam*log1p(-exp(-mpfr(u^kappa,precBits))))))
s1<-as.numeric(log1p(-exp(-mpfr(u^kappa,precBits))))
}
lf<-status*(log(kappa)+log(gam)+log(rho)+(kappa-1)*log(u)+(gam-1)*s1-u^kappa-pred-r1)+r1
return(-sum(lf[is.finite(lf)],na.rm=TRUE))
}
#######################################################
# Fit of the Exp weibull AFT.
#################################################
ewaft.fit<-function(Formula,init=NULL,data,iter.max=150){
variables<-all.vars(Formula)
data0<-data[,variables]
if(any(is.na(data0))){
stop("missing values in the data")
}
Formula1<-as.formula(paste("~",as.character(Formula)[3],collapse=""))
xx0<-model.matrix(Formula1,data)
cnames<-colnames(xx0)[-1]
xx<-as.matrix(xx0[,-1])
colnames(xx)<-cnames
st<-c(as.matrix(data[variables[1]]))
status<-c(as.matrix(data[variables[2]]))
p<-ncol(xx)
if(is.null(init)){
init<-rbind(c(rep(0,p),rep(-1,3)),c(rep(0,p),rep(-1.5,3)),c(rep(0,p),rep(-2,3)),
   c(rep(0,p),rep(0,3)),c(rep(0,p),rep(-3,3)),c(rep(0,p),rep(-4,3)),
   c(rep(0,p),rep(0.5,3)),c(rep(0,p),rep(1,3)))
}
code<-1
k<-0
repeat{
k<-k+1
options(warn=-1)
fit1<-tryCatch({nlminb(start=init[k,],objective=ewaft.llik,st=st,
  status=status,xx=xx,control=list(iter.max=iter.max))},error = function(e){NULL})
options(warn=0)
if(!is.null(fit1)){
if(fit1$message=="relative convergence (4)" & fit1$convergence==0 & is.finite(fit1$objective) & 
  fit1$objective!=0 & !is.na(fit1$objective)){
options(warn=-1)
hess<-hessian(func=ewaft.llik,x=fit1$par,st=st,
  status=status,xx=xx)
options(warn=0)
if(all(is.finite(hess))){
if(!any(eigen(hess)$values<=0)){
logL<--fit1$objective
code<-0
optimizer<-"nlminb"}}}} 
if(k==nrow(init) || code==0) break}
if(code==1){
k<-0
repeat{
k<-k+1
options(warn=-1)
fit1<-tryCatch({nlminb(start=init[k,],objective=ewaft.llik,st=st,
  status=status,xx=xx,precBits=120,control=list(iter.max=iter.max))},error = function(e){NULL})
options(warn=0)
if(!is.null(fit1)){
if(fit1$message=="relative convergence (4)" & fit1$convergence==0 & is.finite(fit1$objective) & 
  fit1$objective!=0 & !is.na(fit1$objective)){
options(warn=-1)
hess<-hessian(func=ewaft.llik,x=fit1$par,st=st,
  status=status,xx=xx,precBits=120)
options(warn=0)
if(all(is.finite(hess))){
if(!any(eigen(hess)$values<=0)){
logL<--fit1$objective
code<-0
optimizer<-"nlminb"}}}} 
if(k==nrow(init) || code==0) break}
}
if(code==1){
k<-0
repeat{
k<-k+1
options(warn=-1)
fit1<-tryCatch({optim(par=init[k,],fn=ewaft.llik,st=st,
  status=status,xx=xx,method="BFGS",hessian=TRUE,control=list(maxit=iter.max))},
  error = function(e){NULL})
options(warn=0)
if(!is.null(fit1)){
hess<-fit1$hessian
if(all(is.finite(hess)) & fit1$convergence==0 & is.finite(fit1$value) & fit1$value!=0 & !is.na(fit1$value)){
if(!any(eigen(hess)$values<=0)){
code<-0
logL<--fit1$value
optimizer<-"optim"}}}
if(k==nrow(init) || code==0) break}}
if(code==1){
return(list(code=1,message="non-convergence; try with a different set of initial values"))
} else{
est<-fit1$par
cov.mat<-solve(hess)
rownames(cov.mat)<-colnames(cov.mat)<-c(colnames(xx),"log(kappa)","log(gamma)","log(rho)")
se<-sqrt(diag(cov.mat))
z<-est/se
ci.lo<-est-se*qnorm(0.975)
ci.up<-est+se*qnorm(0.975)
pval<-2*pnorm(-abs(z))
result1<-data.frame(est=est,se=se,z=z,
  "p value"=ifelse(pval<0.001,"<0.001",round(pval,digits=3)),lower.95=ci.lo,upper.95=ci.up)
rownames(result1)<-c(colnames(xx),"log(kappa)","log(gamma)","log(rho)")
result2<-cbind("exp(est)"=exp(est)[1:p],"lower.95"=exp(ci.lo)[1:p],upper.95=exp(ci.up)[1:p])
rownames(result2)<-colnames(xx)
result3<-cbind("exp(-est)"=exp(-est)[1:p],"lower.95"=exp(-ci.up)[1:p],upper.95=exp(-ci.lo)[1:p])
rownames(result3)<-colnames(xx)
dev<--2*logL
aic<-dev+2*length(est)
result4<-cbind(logL=logL,deviance=dev,AIC=aic)
rownames(result4)<-""
dat.sum<-cbind(n=length(st),"events"=sum(status),predictors=p)
rownames(dat.sum)<-""
model<-noquote(rbind("Exponentiated Weibull AFT","rho = rate parameter", "kappa, gamma = shape parameters"))
colnames(model)<-""
rownames(model)<-rep("",3)
results<-list(Model=model,
  "Data summary"=dat.sum,Fit=result1,"exp(est)"=result2,
  "exp(-est)"=result3,"Fit criteria"=result4,optimizer=optimizer,cov=cov.mat,
  st=st,status=status,design.mat=xx,surv.model="ewaft",code=0)
}
return(results)
}
##################################################
# Residual plot for eweibull
#################################################
ewaft.resid<-function(fit,plot=FALSE,xlim=NULL,ylim=NULL,main=NULL){
st<-fit$st
status<-fit$status
xx<-as.matrix(fit$design.mat)
p<-ncol(xx)
beta<-fit$Fit[1:p,1]
kappa<-exp(fit$Fit[(p+1),1])
gam<-exp(fit$Fit[(p+2),1])
rho<-exp(fit$Fit[(p+3),1])
pred<-as.matrix(xx)%*%beta
u<-rho*st*exp(-pred)
Rhat<--as.numeric(log1p(-exp(gam*log1p(-exp(-mpfr(u^kappa,120))))))
km.Rhat<-survfit(Surv(Rhat,status)~1, type="kaplan-meier")
if(plot){
H<-cbind(km.Rhat$time,-log(km.Rhat$surv))
H[which(!is.finite(H))]<-NA
H<-na.omit(H)
if(is.null(xlim) || is.null(ylim)){
xlim<-ylim<-c(min(c(H)),max(c(H)))
}
plot(H[,1],H[,2],type="p",pch=20,xlab="Residuals",
  ylab="Estimated Cumulative Hazard Function",
  xlim=xlim,ylim=ylim,main="",
  cex.lab=1.25,cex.main=1.5)
abline(a=0,b=1)
}
return(list(residuals=as.vector(km.Rhat$time)))
}
###################################################
# Recurrent: Negative log-likelihood of Exp Weibull AFT
####################################################
Rewaft.llik<-function(init,tstart,tstop,status,xx,precBits=NULL){
p<-ncol(xx)
beta<-init[1:p]
kappa<-exp(init[p+1])
gam<-exp(init[p+2])
rho<-exp(init[p+3])
pred<-as.matrix(xx)%*%beta
u0<-rho*tstart*exp(-pred)
u1<-rho*tstop*exp(-pred)
if(is.null(precBits)){
r0<-log1p(-exp(gam*log1p(-exp(-u0^kappa))))
r1<-log1p(-exp(gam*log1p(-exp(-u1^kappa))))
s1<-log1p(-exp(-u1^kappa))
} else{
r0<-as.numeric(log1p(-exp(gam*log1p(-exp(-mpfr(u0^kappa,precBits))))))
r1<-as.numeric(log1p(-exp(gam*log1p(-exp(-mpfr(u1^kappa,precBits))))))
s1<-as.numeric(log1p(-exp(-mpfr(u1^kappa,precBits))))
}
lf<-status*(log(kappa)+log(gam)+log(rho)+(kappa-1)*log(u1)+(gam-1)*s1-u1^kappa-pred-r1)-(r0-r1)
return(-sum(lf[is.finite(lf)],na.rm=TRUE))
}
#######################################################
# Recurrent: Fit of the Exp weibull AFT.
#################################################
Rewaft.fit<-function(Formula,init=NULL,data,iter.max=150){
variables<-all.vars(Formula)
data0<-data[,variables]
if(any(is.na(data0))){
stop("missing values in the data")
}
Formula1<-as.formula(paste("~",as.character(Formula)[3],collapse=""))
xx0<-model.matrix(Formula1,data)
cnames<-colnames(xx0)[-1]
xx<-as.matrix(xx0[,-1])
colnames(xx)<-cnames
tstart<-c(as.matrix(data[variables[1]]))
tstop<-c(as.matrix(data[variables[2]]))
if(any(tstop<=tstart)){
stop("stop time must be greater than the start time")
}
status<-c(as.matrix(data[variables[3]]))
p<-ncol(xx)
if(is.null(init)){
init<-rbind(c(rep(0,p),rep(-1,3)),c(rep(0,p),rep(-1.5,3)),c(rep(0,p),rep(-2,3)),
   c(rep(0,p),rep(0,3)),c(rep(0,p),rep(-3,3)),c(rep(0,p),rep(-4,3)),
   c(rep(0,p),rep(0.5,3)),c(rep(0,p),rep(1,3)))
}
code<-1
k<-0
repeat{
k<-k+1
options(warn=-1)
fit1<-tryCatch({nlminb(start=init[k,],objective=Rewaft.llik,tstart=tstart,tstop=tstop,
  status=status,xx=xx,control=list(iter.max=iter.max))},error = function(e){NULL})
options(warn=0)
if(!is.null(fit1)){
if(fit1$message=="relative convergence (4)" & fit1$convergence==0 & is.finite(fit1$objective) & 
  fit1$objective!=0 & !is.na(fit1$objective)){
options(warn=-1)
hess<-hessian(func=Rewaft.llik,x=fit1$par,tstart=tstart,tstop=tstop,
  status=status,xx=xx)
options(warn=0)
if(all(is.finite(hess))){
if(!any(eigen(hess)$values<=0)){
logL<--fit1$objective
code<-0
optimizer<-"nlminb"}}}} 
if(k==nrow(init) || code==0) break}
if(code==1){
k<-0
repeat{
k<-k+1
options(warn=-1)
fit1<-tryCatch({nlminb(start=init[k,],objective=Rewaft.llik,tstart=tstart,tstop=tstop,
  status=status,xx=xx,precBits=120,control=list(iter.max=iter.max))},error = function(e){NULL})
options(warn=0)
if(!is.null(fit1)){
if(fit1$message=="relative convergence (4)" & fit1$convergence==0 & is.finite(fit1$objective) & 
  fit1$objective!=0 & !is.na(fit1$objective)){
options(warn=-1)
hess<-hessian(func=Rewaft.llik,x=fit1$par,tstart=tstart,tstop=tstop,
  status=status,xx=xx,precBits=120)
options(warn=0)
if(all(is.finite(hess))){
if(!any(eigen(hess)$values<=0)){
logL<--fit1$objective
code<-0
optimizer<-"nlminb"}}}} 
if(k==nrow(init) || code==0) break}
}
if(code==1){
k<-0
repeat{
k<-k+1
options(warn=-1)
fit1<-tryCatch({optim(par=init[k,],fn=Rewaft.llik,tstart=tstart,tstop=tstop,
  status=status,xx=xx,method="BFGS",hessian=TRUE,control=list(maxit=iter.max))},
  error = function(e){NULL})
options(warn=0)
if(!is.null(fit1)){
hess<-fit1$hessian
if(all(is.finite(hess)) & fit1$convergence==0 & is.finite(fit1$value) & fit1$value!=0 & !is.na(fit1$value)){
if(!any(eigen(hess)$values<=0)){
code<-0
logL<--fit1$value
optimizer<-"optim"}}}
if(k==nrow(init) || code==0) break}}
if(code==1){
return(list(code=1,message="non-convergence; try with a different set of initial values"))
} else{
est<-fit1$par
cov.mat<-solve(hess)
rownames(cov.mat)<-colnames(cov.mat)<-c(colnames(xx),"log(kappa)","log(gamma)","log(rho)")
se<-sqrt(diag(cov.mat))
z<-est/se
ci.lo<-est-se*qnorm(0.975)
ci.up<-est+se*qnorm(0.975)
pval<-2*pnorm(-abs(z))
result1<-data.frame(est=est,se=se,z=z,
  "p value"=ifelse(pval<0.001,"<0.001",round(pval,digits=3)),lower.95=ci.lo,upper.95=ci.up)
rownames(result1)<-c(colnames(xx),"log(kappa)","log(gamma)","log(rho)")
result2<-cbind("exp(est)"=exp(est)[1:p],"lower.95"=exp(ci.lo)[1:p],upper.95=exp(ci.up)[1:p])
rownames(result2)<-colnames(xx)
result3<-cbind("exp(-est)"=exp(-est)[1:p],"lower.95"=exp(-ci.up)[1:p],upper.95=exp(-ci.lo)[1:p])
rownames(result3)<-colnames(xx)
dev<--2*logL
aic<-dev+2*length(est)
result4<-cbind(logL=logL,deviance=dev,AIC=aic)
rownames(result4)<-""
dat.sum<-cbind(n=length(tstop),"events"=sum(status),predictors=p)
rownames(dat.sum)<-""
model<-noquote(rbind("Exponentiated Weibull AFT","rho = rate parameter", "kappa, gamma = shape parameters"))
colnames(model)<-""
rownames(model)<-rep("",3)
results<-list(Model=model,
  "Data summary"=dat.sum,Fit=result1,"exp(est)"=result2,
  "exp(-est)"=result3,"Fit criteria"=result4,optimizer=optimizer,cov=cov.mat,
  st=cbind(start=tstart,stop=tstop),status=status,design.mat=xx,
  surv.model="Rewaft",code=0)
}
return(results)
}
##################################################
# Recurrent: Residual plot for eweibull
#################################################
Rewaft.resid<-function(fit,plot=FALSE,xlim=NULL,ylim=NULL,main=NULL){
tstart<-fit$st[,1]
tstop<-fit$st[,2]
status<-fit$status
xx<-as.matrix(fit$design.mat)
p<-ncol(xx)
beta<-fit$Fit[1:p,1]
kappa<-exp(fit$Fit[(p+1),1])
gam<-exp(fit$Fit[(p+2),1])
rho<-exp(fit$Fit[(p+3),1])
pred<-as.matrix(xx)%*%beta
u0<-rho*tstart*exp(-pred)
u1<-rho*tstop*exp(-pred)
Rhat1<--as.numeric(log1p(-exp(gam*log1p(-exp(-mpfr(u0^kappa,120))))))
Rhat2<--as.numeric(log1p(-exp(gam*log1p(-exp(-mpfr(u1^kappa,120))))))
km.Rhat<-survfit(Surv(Rhat1,Rhat2,status)~1, type="kaplan-meier")
if(plot){
H<-cbind(km.Rhat$time,-log(km.Rhat$surv))
H[which(!is.finite(H))]<-NA
H<-na.omit(H)
if(is.null(xlim) || is.null(ylim)){
xlim<-ylim<-c(min(c(H)),max(c(H)))
}
plot(H[,1],H[,2],type="p",pch=20,xlab="Residuals",
  ylab="Estimated Cumulative Hazard",
  xlim=xlim,ylim=ylim,main=main,
  cex.lab=1.25,cex.main=1.5)
abline(a=0,b=1)
}
return(list(residuals=as.vector(km.Rhat$time)))
}