# Kentaro Fukumoto rm(list=ls()) ######################################### # Survival Analysis: Weibull distribution sum.log.failure.weibull <- function(parameter,time,x){ # Box and Jones pp. 26-7 # lambda = exp(xb) xb <- x%*%parameter[-1] p <- exp(parameter[1]) sum((p*xb) + ((p-1)*log(time)) - ((exp(xb)*time)^p) + parameter[1]) } sum.log.Survival.weibull <- function(parameter,time,x){ sum(-((exp(x%*%parameter[-1])* time)^exp(parameter[1]))) } log.likelihood.weibull <- function(parameter,end.time,start.time,event,x){ sum.log.failure.weibull(parameter,end.time[event==1],x[event==1,]) + # failure sum.log.Survival.weibull(parameter,end.time[event==0],x[event==0,]) - # right censoring sum.log.Survival.weibull(parameter,start.time,x) # left truncation } EHA.weibull <- function(formula,Event,Start.time){ model <- model.frame(formula) End.time <- as.matrix(model.extract(model, "response")) X <- as.matrix(model.matrix(formula)) k <- ncol(X) parameter.init <- rep(0,(k+1)) optim.out <- optim(parameter.init,log.likelihood.weibull,method = "BFGS", hessian = TRUE, control = list(fnscale = -1),end.time=End.time,x=X,event=Event,start.time=Start.time) names(optim.out$par) <- c("ancillary",colnames(X)) optim.out } summarize <- function(MLE.out,digit){ point.estimate <- MLE.out$par standard.error <- sqrt(diag(solve(-MLE.out$hessian))) p.value <- pnorm(-abs(MLE.out$par)/standard.error) round(cbind(point.estimate,standard.error,p.value),digit) } AIC <- function(optim.out){ -2*(length(optim.out$par)+optim.out$value) } # application library(survival) data() attach(jasa1) ?jasa1 survreg.out <- survreg(Surv(stop, event) ~ transplant + surgery + age) summary.survreg.out <- round(summary(survreg.out)$table,3) summary.survreg.out EHA.weibull.out <- EHA.weibull(stop ~ transplant + surgery + age, event,0) summarize.EHA.weibull.out <- summarize(EHA.weibull.out,3) summarize.EHA.weibull.out AIC(EHA.weibull.out) # left truncation EHA.weibull.out.2 <- EHA.weibull(stop ~ transplant + surgery + age, event,start) summarize.EHA.weibull.out.2 <- summarize(EHA.weibull.out.2,3) summarize.EHA.weibull.out.2 AIC(EHA.weibull.out.2) ################## # split population logistic <- function(x){1/(1+exp(-x))} sum.log.failure.weibull.splitpop <- function(parameter,time,x){ xb <- x%*%parameter[-c(1,2)] p <- exp(parameter[1]) sum((p*xb) + ((p-1)*log(time)) - ((exp(xb)*time)^p) + parameter[1] + log(logistic(parameter[2]))) # splitpop } Survival.weibull <- function(parameter,time,x){ exp(-((exp(x%*%parameter[-c(1,2)])* time)^exp(parameter[1]))) } log.likelihood.weibull.splitpop <- function(parameter,end.time,start.time,event,x){ delta <- logistic(parameter[2]) sum.log.failure.weibull.splitpop(parameter,end.time[event==1],x[event==1,])- sum(log(Survival.weibull(parameter,start.time[event==1],x[event==1,]))) + sum(log(1-delta+(delta*(Survival.weibull(parameter,end.time[event==0],x[event==0,])/ Survival.weibull(parameter,start.time[event==0],x[event==0,]))))) } EHA.weibull.splitpop <- function(formula,Event,Start.time){ model <- model.frame(formula) End.time <- as.matrix(model.extract(model, "response")) X <- as.matrix(model.matrix(formula)) k <- ncol(X) parameter.init <- rep(0,(k+2)) # splitpop optim.out <- optim(parameter.init,log.likelihood.weibull.splitpop,method = "BFGS", hessian = TRUE, control = list(fnscale = -1),end.time=End.time,x=X,event=Event,start.time=Start.time) names(optim.out$par) <- c("ancillary","splitpop",colnames(X)) optim.out } EHA.weibull.out.splitpop <- EHA.weibull.splitpop(stop ~ transplant + surgery + age, event,start) summarize.EHA.weibull.out.splitpop <- summarize(EHA.weibull.out.splitpop,3) summarize.EHA.weibull.out.splitpop AIC(EHA.weibull.out.splitpop)