R/rss.ELR.test.R

In generalRSS: Statistical Tools for Balanced and Unbalanced Ranked Set Sampling

#' @details This function performs a one-sample empirical likelihood ratio (ELR) test on ranked set sample data using the method introduced by Ahn et al. (2024). Given a data frame of RSS data `data` with `rank` and `y` columns, the function calculates the empirical likelihood ratio test statistic, confidence interval, and p-value based on the hypothesized mean value `mu0`.
#' @title RSS empirical likelihood ratio (ELR) test for one-sample population mean
#' @name rss.ELR.test
#' @description The rss.ELR.test function conducts a one-sample empirical likelihood ratio test on ranked set sample data to assess the population mean.
#'
#' @param data A numeric data frame of ranked set samples with columns `rank` for ranks and `y` for data values.
#' @param alpha A numeric value specifying the confidence level for the interval.
#' @param mu0 A numeric value indicating the hypothesized value of the mean.
#'
#' @return
#'   \item{RSS_mean}{The RSS mean estimate.}
#'   \item{CI}{The confidence interval for the population mean.}
#'   \item{-2*Log.LR}{The empirical log likelihood ratio test statistic.}
#'   \item{p.value}{The p-value for the test.}
#'
#' @references
#'
#' S. Ahn, X. Wang, C. Moon, and J. Lim. (2024) New scheme of empirical likelihood method for ranked set sampling: Applications to two one sample problems. International Statistical Review.
#'
#' @seealso
#' \code{\link{rss.simulation}}: used for simulating Ranked Set Samples (RSS), which can serve as input.
#'
#' \code{\link{rss.sampling}}: used for sampling Ranked Set Samples (RSS) from a population data set, providing input data.
#'
#' @examples
#' ## Unbalanced RSS with a set size 3 and different sample sizes of 6, 10, and 8 for each stratum,
#' ## using imperfect ranking from a normal distribution with a mean of 0.
#' rss.data<-rss.simulation(H=3,nsamp=c(6,10,8),dist="normal", rho=0.8,delta=0)
#'
#' ## RSS empirical likelihood ratio test
#' \donttest{
#' rss.ELR.test(data=rss.data, alpha=0.05, mu0=0)
#' }
#'
#' @export
rss.ELR.test <- function(data, alpha=0.05, mu0){

  if( (alpha > 0) & (alpha < 1)){
    if(!all(c("rank", "y") %in% colnames(data))) {
      stop("The input data must contain 'rank' and 'y' variables.")
    }

    H = length(unique(data$rank))
    nsamp = table(data$rank)
    rss.mu=mean(tapply(data$y,data$rank,mean))

    chi.level = stats::qchisq(1-alpha,1)
    result=try(emplik::findUL2(step=0.005,fun=rss.EL.solve, MLE=rss.mu, level=chi.level, x=data))
    #if(class(result)=="try-error"){
    if(inherits(result, "try-error")){
      rss.el=list(Low=NA, Up=NA)
    }else{
      rss.el = emplik::findUL2(step=0.005,fun=rss.EL.solve, MLE=rss.mu, level=chi.level, x=data)
    }

    CI.low = rss.el$Low
    CI.up = rss.el$Up

    LLR = rss.EL.solve(mu0,data)$'-2LLR'
    if( LLR == "NA"){
      warning("Warning: Convergence failed due to extreme deviation, leading to infinite LRT and p-value = 0.")
      pval = 0
      LLR = "Inf"
    }else{
      pval = 1-stats::pchisq(LLR, df=1)
    }

    result <- list(RSS_mean = rss.mu, CI = c(CI.low, CI.up), "-2*log.LR" = LLR, p.value = pval)
    return(result)

  }else stop("alpha is out of bound.", call. = F)

}

#' @noRd
rss.EL.solve<-function(mu0,x){
  data=x
  H = length(unique(data$rank))
  nh = table(data$rank)
  nhsum=sum(nh)

  model <- function(x) {
    F.nu <- sum(sapply(1:H, function(h) {
      sum((data$y[data$rank == h] - mu0) / (x[h] + x[H + 1] * (data$y[data$rank == h] - mu0)))
    }))
    F.lambda <- sapply(1:H, function(h) {
      sum(1 / (x[h] + x[H + 1] * (data$y[data$rank == h] - mu0))) - 1 / H
    })
    c(F.lambda, F.nu)
  }

  init.nu<-seq(-100,100,by=5)
  i<-1
  while(i <= length(init.nu)){

    nu<- init.nu[i]
    lambda.h<-H*nh

    sol = tryCatch(rootSolve::multiroot(f=model, start=c(lambda.h,nu))$root,
                   warning=function(w){
                     if(grepl("diagonal element is zero", w$message)){
                       invokeRestart("muffleWarning")
                     }
                   }
    )
    if(is(sol)[1]=="NULL"){i<-i+1; next}
    if(sum(is.na(sol))>0){i<-i+1; next}

    lambda.h<-sol[1:H]; nu<-sol[H+1]

    lambda.h.vec<-rep(lambda.h,nh)
    prh<-1/(lambda.h.vec+nu*(data$y-mu0))

    if(sum(prh<0)==0){break}
    i<-i+1
  }
  if( sum(prh<0) !=0 ){
    new.LLR <- "NA"
  }else{
    new.LLR <- 2*sum(log(1/prh/(H*rep(nh,nh))))
  }
  return(list("-2LLR"=new.LLR))
}