R/AFTlrtest.R

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

####################################################
# LR test for weibull as a submodel of eweibull or ggamma
###################################################
#' Likelihood ratio test for Weibull as a submodel of the exponentiated Weibull or 
#'       generalized gamma model
#' @description Likelihood ratio test for Weibull as a submodel of the exponentiated Weibull or
#'      generalized gamma model. 
#' @keywords Exponentiated Weibull, Generalized gamma, Likelihood ratio test, Weibull
#' @param fit a list of two \code{survreg.aft} fits: one must be the Weibull fit and the other one is the exponentiated
#'      Weibull or generalized gamma fit.
#' @details The log-likelihood values from the two fits are used for the likelihood ratio test.
#' @return \code{chi.sq:} chi-square test statistic.
#' @return \code{p value:} p value.
#' @author Shahedul Khan <khan@math.usask.ca>
#' @seealso \code{\link{survreg.aft}}, \code{\link{rresid.plot}}.
#' @examples 
#'   library(frailtypack)
#'   data(readmission)
#' # Recoding sex and chemo
#'   sex<-factor(2-as.numeric(readmission$sex),
#'        labels=c("Female","Male"))
#'   chemo<-factor(as.numeric(readmission$chemo)-1,
#'          labels=c("NonTreated","Treated"))
#'   Data<-readmission
#'   Data$sex<-sex
#'   Data$chemo<-chemo
#'   fit.gg<-survreg.aft(c(t.start,t.stop,event)~sex+chemo+charlson,
#'       Data=Data,model="ggamma")
#'   fit.ew<-survreg.aft(c(t.start,t.stop,event)~sex+chemo+charlson,
#'       Data=Data,model="eweibull")
#'   fit.w<-survreg.aft(c(t.start,t.stop,event)~sex+chemo+charlson,
#'       Data=Data,model="weibull")
#'   LR.test(list(fit.w,fit.ew))
#'   LR.test(list(fit.w,fit.gg))
#' @export

LR.test<-function(fit){
for(i in 1:length(fit)){
if(fit[[i]]$surv.model=="Rwaft" || fit[[i]]$surv.model=="waft") {L2<-fit[[i]][[6]][1]}
if(fit[[i]]$surv.model=="Rewaft" || fit[[i]]$surv.model=="ewaft" ||
   fit[[i]]$surv.model=="Rggaft" || fit[[i]]$surv.model=="ggaft") {L1<-fit[[i]][[6]][1]}
}
stat<-2*(L1-L2)
p.val<-pchisq(stat,1,lower.tail=FALSE)
results<-cbind(chi.sq=stat,"p value"=p.val)
rownames(results)<-""
return(results)
}