R/compareSurveysLRT.R

In fitPS: Fit Zeta Distributions to Forensic Data

#' Compare two or more surveys on the basis of their shape parameters using a Likelihood Ratio Test
#'
#' @param \ldots two or more objects of class \code{"psData"}---see \code{\link{readData}}.
#'
#' @details
#' This function **only** works for the zeta distribution. The function carries out a likelihood ratio test (LRT) to test the null hypothesis
#' \deqn{H_0: \alpha_1 = \alpha_2 = \ldots = \alpha_K}{H_0: alpha_1 = alpha_2 = \ldots = alpha_K} versus the alternative
#' \deqn{H_1: \alpha_i \neq \alpha_j \mbox{ for some } i \neq j \in \left\{1, \ldots, K\right\},}{H_1: alpha_i != alpha_j for some  i != j in {1, \ldots, K},}
#' where \eqn{\alpha_i}{alpha_i} is the shape parameter for the zeta distribution of the \eqn{i^\mathrm{th}}{ith} survey.
#'
#' @return The function returns a \code{list} of class \code{"htest"} with the following elements:
#' \describe{
#'  \item{\code{statistic}}{ -- the test statistic.}
#'  \item{\code{parameter}}{ -- the degrees of freedom for the test}
#'  \item{\code{p.value}}{ -- the P-value associated with the estimate.}
#'  \item{\code{method}}{ -- a character string describing the method hypothesis.}
#'  \item{\code{data.name}}{ -- the names of the data sets used in the test}
#' }
#'
#' @examples
#' data(Psurveys)
#' lau = Psurveys$lau
#' jackson = Psurveys$jackson
#' compareSurveysLRT(lau, jackson)
#'
#' ## Example with three surveys
#' roux = Psurveys$roux
#' compareSurveysLRT(lau, jackson, roux)
#'
#' @importFrom stats pchisq
#'
#' @export
compareSurveysLRT = function(...){
  Surveys = list(...)

  if(length(Surveys) < 2){
    stop("You need to supply at least two surveys to carry out this test")
  }

  if(any(sapply(Surveys, function(x){
    !is(x, "psData")
  }))){
    stop("All the arguments supplied to this function must be of class psData")
  }

  types = sapply(Surveys, function(x){
    x$type
  })

  if(length(unique(types)) > 1){
    stop("All the surveys must be of the same type - either P or S")
  }

  fits = lapply(Surveys, fitDist)

  ll.mle = sum(sapply(fits, function(x){
    -x$fit$value
  }))

  ll.H0 = -fitDist(combineSurveys(...))$fit$value

  LRT.stat = 2 * (ll.mle - ll.H0)
  names(LRT.stat) = "X-squared"
  parameter = length(Surveys) - 1
  names(parameter) = "df"

  P = pchisq(LRT.stat, df = length(Surveys) - 1, lower.tail = FALSE)
  dname = paste(unlist(match.call(expand.dots = FALSE)$...), collapse = ", ")

  rval = list(statistic = LRT.stat,
              parameter = length(Surveys) - 1,
              p.value = P,
              method = "Likelihood Ratio Test",
              data.name = dname
  )
  class(rval) = "htest"
  return(rval)
}