R/AFTlnormal.R

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

########################################
# Negative log-likelihood of Log-normal AFT
###########################################
lnaft.llik<-function(init,st,status,xx){
p<-ncol(xx)
beta<-init[1:p]
kappa<-exp(init[p+1])
rho<-exp(init[p+2])
pred<-as.matrix(xx)%*%beta
u<-st*exp(-pred)
log.st<-plnorm(u,(-log(rho)),(1/kappa),lower.tail = FALSE,log.p=TRUE)
log.ft<-dlnorm(u,(-log(rho)),(1/kappa),log=TRUE)-pred
log.ht<-status*(log.ft-log.st)
Ht<--log.st
lf<-log.ht-Ht
return(-sum(lf[is.finite(lf)],na.rm=TRUE))
}
#######################################################
# Fit of the log-normal AFT.
#################################################
lnaft.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,2)),c(rep(0,p),rep(-1.5,2)),c(rep(0,p),rep(-2,2)),
   c(rep(0,p),rep(0,2)),c(rep(0,p),rep(-3,2)),c(rep(0,p),rep(-4,2)),
   c(rep(0,p),rep(0.5,2)),c(rep(0,p),rep(1,2)))
}
code<-1
k<-0
repeat{
k<-k+1
options(warn=-1)
fit1<-tryCatch({nlminb(start=init[k,],objective=lnaft.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=lnaft.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({optim(par=init[k,],fn=lnaft.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(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(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("Log-normal AFT","rho = rate parameter", "kappa = shape parameter"))
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="lnaft",code=0)
}
return(results)
}
##################################################
# Residual plot for log-normal
#################################################
lnaft.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])
rho<-exp(fit$Fit[(p+2),1])
pred<-as.matrix(xx)%*%beta
u<-st*exp(-pred)
Rhat<--plnorm(u,(-log(rho)),(1/kappa),lower.tail = FALSE,log.p=TRUE)
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=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 log-normal AFT
####################################################
Rlnaft.llik<-function(init,tstart,tstop,status,xx,precBits=NULL){
p<-ncol(xx)
beta<-init[1:p]
kappa<-exp(init[p+1])
rho<-exp(init[p+2])
pred<-as.matrix(xx)%*%beta
u0<-tstart*exp(-pred)
u1<-tstop*exp(-pred)
log.st0<-plnorm(u0,(-log(rho)),(1/kappa),lower.tail = FALSE,log.p=TRUE)
log.st1<-plnorm(u1,(-log(rho)),(1/kappa),lower.tail = FALSE,log.p=TRUE)
log.ft1<-dlnorm(u1,(-log(rho)),(1/kappa),log=TRUE)-pred
log.ht1<-status*(log.ft1-log.st1)
Ht<-log.st0-log.st1
lf<-log.ht1-Ht
return(-sum(lf[is.finite(lf)],na.rm=TRUE))
}
#####################################################
# Recurrent: Fit of the log-logistic AFT.
#################################################
Rlnaft.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(-0.5,2)),c(rep(0,p),rep(-1,2)),c(rep(0,p),rep(0,2)),
   c(rep(0,p),rep(-1.5,2)),c(rep(0,p),rep(-2,2)),c(rep(0,p),rep(-3,2)),c(rep(0,p),rep(-4,2)),
   c(rep(0,p),rep(0.5,2)),c(rep(0,p),rep(1,2)))
}
code<-1
k<-0
repeat{
k<-k+1
options(warn=-1)
fit1<-tryCatch({nlminb(start=init[k,],objective=Rlnaft.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=Rlnaft.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({optim(par=init[k,],fn=Rlnaft.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(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(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("Log-normal AFT","rho = rate parameter", "kappa = shape parameter"))
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="Rlnaft",code=0)
}
return(results)
}
##################################################
# Recurrent: Residual plot for log-normal
#################################################
Rlnaft.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])
rho<-exp(fit$Fit[(p+2),1])
pred<-as.matrix(xx)%*%beta
u0<-tstart*exp(-pred)
u1<-tstop*exp(-pred)
log.st0<-plnorm(u0,(-log(rho)),(1/kappa),lower.tail = FALSE,log.p=TRUE)
log.st1<-plnorm(u1,(-log(rho)),(1/kappa),lower.tail = FALSE,log.p=TRUE)
Rhat1<--log.st0
Rhat2<--log.st1
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)))
}