#' @title Fitting one or a set of failure time models
#' @description Code for fitting a particular tag-life model or models
#'
#'
#' @param time numeric vector of failure times
#' @param model character vector of specified model(s)
#' @param rc.value right-censoring cutoff value; only observations with times > rc.value are considered to have failed prematurely.
#' @param rt.value right-truncation cutoff value; only observations with times < rt.value will be included in the model fitting
#'
#' @return A failure model object, if one model specified OR a failure model list object if multiple models are specified
#'
#' @details
#' Model fitting routine used to fit one or a set of failure time models
#'
#' If a single model is specified a "failmod_obj" is created, which can be
#' used to adjust a CJS model in the "ATLAS" package.
#'
#' If multiple models are specified a "failmod_list" is created containing
#' output from all model fits. This object serves as an input in the
#' fail_rank() function, which ranks the performance of the model using
#' the \href{http://animalbiotelemetry.biomedcentral.com/articles/10.1186/s40317-020-00213-z}{Skalski and Whitlock (2020)} GOF measure.
#'
#' (named list)
#' "mod_choice" = character object of model name or vector of model names,
#' "fit_vals" = data frame with fitted values for each time and 95% confidence interval
#' "mod_objs" = model objects from various packages: kaplan-meiser ("surival"), vitality models ("vitality")
#' "par_tab" = dataframe of one or more model's parameter estimates
#'
#' @importFrom survival Surv
#' @importFrom flexsurv flexsurvreg
#' @importFrom vitality vitality.ku
#' @importFrom vitality vitality.4p
#' @export fail_fit
#'
#'
fail_fit=function(time,model,rc.value=NULL,rt.value=NULL,...){
rc=FALSE #temp def
rt=FALSE #temp def
if(!is.vector(time)|!is.numeric(time)){stop("A numeric vector is expected for the 'time' argument")}
if(model[1]=="all"){model=c("weibull",'weibull3', "gompertz", "gamma", "lognormal", "llogis", "gengamma","vitality.ku","vitality.4p")
message("Fitting all default parametric survival models")}
# WARNINGS AND VALIDATION
if(any(is.na(time))){
message(paste(length(which(is.na(time))),"NA times removed"))}
test=!model %in% c("weibull",'weibull3', "gompertz", "gamma", "lognormal", "llogis", "gengamma","vitality.ku","vitality.4p","kaplan-meier")
if(any(test)){
message(paste("Model names not recognized:",paste(model[which(test)],collapse=";")," \n Default model names = {'weibull','weibull3','gompertz','gamma','lognormal','llogis','gengamma','vitality.ku','vitality.4p','kaplan-meier'}",sep=""))
model=model[!test]}
y=sort(time) # sorted data necessary for Vitality package functions
y_sfrac=sapply(y,function(x){1-length(which(y<=x))/length(y)}) # survival fraction calc
non_cen=rep(1,length(y))
if(!is.null(rt.value)){
non_trunc=time<rt.value # vector used by "flexsurv"
y=y[non_trunc]
non_cen=non_cen[non_trunc]
# recalculate the survival fraction after truncating values beyond the threshold
y_sfrac=sapply(y,function(x){1-length(which(y<=x))/length(y)}) # survival fraction calc
}
if(!is.null(rc.value)){
rc=TRUE # change this value for later if statement
non_cen=ifelse(time<rc.value,1,0) # vector used by "flexsurv"
y=sort(y) # sorted data necessary for Vitality package functions
y_sfrac=sapply(y,function(x){1-length(which(y<=x))/length(y)}) # survival fraction calc
# For vitality model
y_cen=y[y<rc.value]
y_cen_sfrac=y_sfrac[y<rc.value]
n_cen=length(y)
}
KM_mod=survival::survfit(survival::Surv(time=y,event=non_cen)~1)
KM_sls=summary(KM_mod)
KM_DF=rbind(data.frame(model="kaplan-meier",time=0,est=1,lcl=1,ucl=1),
data.frame(model="kaplan-meier",time=KM_sls$time,est=KM_sls$surv,lcl=KM_sls$lower,ucl=KM_sls$upper))
KM_DF[nrow(KM_DF),c("lcl","ucl")]=c(0,0) # replacing last two NAs wih 0s
if("kaplan-meier" %in% model & length(model)>1){
model=model[model!="kaplan-meier"]
message("'kaplan-meier' cannot be included in a failmod_list with parametric models, and will be removed")
}
fit=list()
fit_vals=NULL
for (i in 1:length(model)){
# FITTING DISTRIBUTIONS IN THE FLEXSURV PACKAGE
if(model[i] %in% c("weibull", "gompertz", "gamma", "lognormal", "llogis", "gengamma")){
fit[[model[i]]] <- flexsurv::flexsurvreg(survival::Surv(time=y,event=non_cen) ~ 1, dist = model[i])
preds=summary(fit[[i]])[[1]]
# preds[match(time,preds$time),]
fit_vals=rbind(fit_vals,
data.frame(model=model[i],time=0,est=1,lcl=1,ucl=1),
data.frame(model=model[i],preds[match(y,preds$time),]))
}
else{
if(model[i]=="vitality.ku"){
if(rc){
dTmp=vitality::dataPrep(c(0,y_cen),(n_cen:(n_cen-length(y_cen)))/n_cen,datatype="CUM",rc.data=(n_cen>length(y_cen)))
fit[[model[i]]]=vitality::vitality.ku(dTmp[,"time"],sdata = dTmp[,"sfract"],rc.data = T,pplot =F,lplot = F,
se=T,...) # added ... here to get the init.params argument to pass through
pars_tmp=fit[[model[i]]][,"params"]
fit_vals=rbind(fit_vals,
data.frame(model="vitality.ku",
time=c(0,y), # survival proportions
est=c(1,vitality::SurvFn.ku(y,pars_tmp[1],pars_tmp[2],pars_tmp[3],pars_tmp[4])),
lcl=0,ucl=0))
}
else{
fit[[model[i]]] = vitality::vitality.ku(time = sort(y),sdata = y_sfrac,se=T,pplot =F,lplot = F, silent = T)
pars_tmp=fit[[model[i]]][,"params"]
fit_vals=rbind(fit_vals,
data.frame(model="vitality.ku",
time=c(0,y), # survival proportions
est=c(1,vitality::SurvFn.ku(y,pars_tmp[1],pars_tmp[2],pars_tmp[3],pars_tmp[4])),
lcl=0,ucl=0))
}
}
if(model[i]=="vitality.4p"){
if(rc){
dTmp=vitality::dataPrep(c(0,y_cen),(n_cen:(n_cen-length(y_cen)))/n_cen,datatype="CUM",rc.data=(n_cen>length(y_cen)))
fit[[model[i]]]=vitality::vitality.4p(dTmp[,"time"],sdata = dTmp[,"sfract"],rc.data = T,pplot =F,
se=T,...) # added ... here to get the init.params argument to pass through
pars_tmp=fit[[model[i]]][,"params"]
fit_vals=rbind(fit_vals,
data.frame(model="vitality.4p",
time=c(0,y), # survival proportions
est=c(1,vitality::SurvFn.4p(y,pars_tmp[1],pars_tmp[2],pars_tmp[3],pars_tmp[4])),
lcl=0,ucl=0))
}
else{
fit[[model[i]]] = vitality::vitality.4p(time = sort(y),sdata = y_sfrac,se=T,pplot =F,silent = T)
pars_tmp=fit[[model[i]]][,"params"]
fit_vals=rbind(fit_vals,
data.frame(model="vitality.4p",
time=c(0,y), # survival proportions
est=c(1,vitality::SurvFn.4p(y,pars_tmp[1],pars_tmp[2],pars_tmp[3],pars_tmp[4])),
lcl=0,ucl=0))
}
}
if(model[i]=="weibull3"){
tmp=taglife.fn_weib3(sort(y),model.in = "weibull")
fit[[model[i]]] = tmp$mod_obj
pars_tmp=tmp$par_tab
fit_vals=rbind(fit_vals,tmp$fit_vals)
}
}
}
# Combined table of parameter estimates for all fitted models (Associated with failmod_list)
if(length(model)>1){
par_ls=sapply(seq_along(fit),function(i){
tmp_mod_nm=names(fit)[i]
if(class(fit[[i]])[1]=="flexsurvreg"){tmp=fit[[i]]$res[,c("est","se")]}
else{tmp=fit[[i]][,c("params","std")]
pnms=as.character(1:dim(tmp)[1]) # number parameters by default
# substitute vitality parameter names
if(tmp_mod_nm=="vitality.ku"){pnms=c("r","s","k","u")}
if(tmp_mod_nm=="vitality.4p"){pnms=c("r","s","lambda","beta")}
if(tmp_mod_nm=="weibull3"){pnms=c("shape","thrsh","scale")}
dimnames(tmp)=list(pnms,c("est","se"))
tmp=data.frame(tmp)
}
return(tmp)
},simplify = FALSE)
tmp_tab=do.call(rbind,par_ls)
param=unlist(sapply(par_ls,rownames,simplify = F))
par_tab=data.frame(model=rep(model,sapply(par_ls,nrow)),param,tmp_tab,row.names = NULL)
out_ls=list("mod_choice"=model,
"times"=data.frame(time=y,surv_frac=y_sfrac),
"fit_vals"=fit_vals,
"mod_objs"=fit,
"par_tab"=par_tab,
"KM_DF"=KM_DF)
out=structure(out_ls,class="failmod_list")
}
# Table of parameter estimates for the single-model case (Associated with failmod_obj)
else{
# if only K-M specified
if(model=="kaplan-meier"){
par_tab=fit_vals=NULL
fit[[1]]=KM_mod
fit_vals=KM_DF
}
# Otherwise
else{
if(model=="vitality.ku" | model=="vitality.4p" | model =="weibull3"){
par_tab=fit[[1]][,c("params","std")]
pnms=as.character(1:dim(par_tab)[1]) # number parameters by default
if(model=="vitality.ku"){ pnms=c("r","s","k","u")}
if(model=="vitality.4p"){pnms=c("r","s","lambda","beta")}
dimnames(par_tab)=list(pnms,c("est","se"))
}
else{
par_tab=fit[[1]]$res[,c("est","se")]
}
}
out_ls=list("mod_choice"=model,
"times"=data.frame(time=y,surv_frac=y_sfrac),
"fit_vals"=fit_vals,
"mod_objs"=fit[[1]],
"par_tab"=par_tab,
"KM_DF"=KM_DF)
out=structure(out_ls,class="failmod_obj")
}
return(out)
# DEV NOTES
# can use a tapply() to execute function for multiple LotIDs
# Consider trying alternative optimization methods
}
# function that limits the amount of output displayed when
# a failmod_obj or a failmod_list is called
#' @title print.failmod_list()
#'
#' @description Printed output for class "failmod_list"
#'
#' @param x an object of class "failmod_list"
#' @param ... ignored
#'
#' @return description of list of models
#'
#' @export
#'
print.failmod_list <- function(x,...){
cat("Failure model list object\n\n")
cat("Contains the following",length(x[["mod_choice"]]),"models: \n",paste(x[["mod_choice"]],collapse = " ; "))
if(is.null(x$"GOF_tab")){ cat("\n\n*this object can used to compare model fit using the function: fail_rank()\n")}
if(!is.null(x$"GOF_tab")){cat("\n\nRanked list\n")
print(x$"GOF_tab")}
invisible(x)
}
#' @title print.failmod_obj()
#'
#' @description Printed output for class "failmod_obj"
#'
#' @param x an object of class "failmod_obj"
#' @param ... ignored
#'
#' @export
#'
print.failmod_obj <- function(x,...){
if(x[["mod_choice"]]=="kaplan-meier"){
cat("Kaplan-Meier estimates for increments between failure times\n")
print(x[["KM_DF"]][,-c(1)])
}
else{
cat(paste(x[["mod_choice"]],"failure model object \n\n"))
cat("Parameter estimates:\n")
print(x[["par_tab"]])
}
cat("\n*This object can be used to adjust survival estimates using the 'ATLAS' package\n")
invisible(x)
}
#' @title Summary of failmod object
#' @details Detailed summary of model objects in failmod object
#'
#' @param object object of class failmod_obj
#' @param ... ignore
#'
#' @return Summary of failmod_obj model fit calls
#' @export
#'
summary.failmod_obj <- function(object,...){
cat("Summary of",paste(object[["mod_choice"]],"failure model object \n\n"))
print(object$"mod_obj")
cat("\n*This object can be used to adjust survival estimates using the 'ATLAS' package\n")
invisible(object)
}
#' @title Sumary of failmod list
#' @details Detailed summary of model objects in failmod list object
#'
#' @param object object of class failmod_list
#' @param ... ignore
#'
#' @return Summary of model fitting calls and GOF rankings (if available)
#' @export
summary.failmod_list <- function(object,...){
cat("Summary failure model list \n\n")
cat("Contains the following",length(object[["mod_choice"]]),"models: \n",paste(object[["mod_choice"]],collapse = " ; "),"\n\n")
print(object$"mod_obj")
if(is.null(object$"GOF_tab")){ cat("\n\n*this object can used to compare model fit using the function: fail_rank()\n")}
if(!is.null(object$"GOF_tab")){cat("\n\nRanked list\n")
print(object$"GOF_tab")}
invisible(object)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.