R/lrt.R

#' Compute LRT for rate models
#'
#' Performs a likelihood ratio test for two models 
#' @param h0 rateModel object that represents the null hypothesis
#' @param hA rateModel object that represents the alternative hypothesis
#' @name lrt
#' @include rateModel-class.R
#' @rdname lrt
#' @return a data table including D-statistic, the p-value, and the degrees of freedom
#' @examples
#' 
#' @export
methods::setGeneric("lrt", function(h0,hA) {
  standardGeneric("lrt")
})

#' @name lrt
#' @rdname lrt
methods::setMethod("lrt", signature(h0 = "rateModel",hA = "rateModel"), function(h0,hA) {
  d=-2*(logLikelihood(obj = h0)-logLikelihood(obj = hA))
  deltaDF=sum(!hA@fixed)-sum(!h0@fixed)
  p=pchisq(d,lower.tail = FALSE,df = deltaDF)
  return(data.table::data.table(D=d,p=p,df=deltaDF))
})
ndukler/epiAllele documentation built on May 5, 2019, 4:50 p.m.