Merge pull request #7 from ropensci/dev

ael wip

.github/workflows/test-coverage.yaml

@@ -2,7 +2,7 @@
 # Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help
 on:
   push:
-    branches: [main, master]
+    branches: [main, master, dev]
   pull_request:
     branches: [main, master]
 

DESCRIPTION

@@ -1,7 +1,7 @@
 Type: Package
 Package: melt
 Title: Multiple Empirical Likelihood Tests
-Version: 1.11.2
+Version: 1.11.3
 Authors@R: c(
     person("Eunseop", "Kim", , "markean@pm.me", role = c("aut", "cph", "cre")),
     person("Steven", "MacEachern", role = c("ctb", "ths")),

NAMESPACE

@@ -29,6 +29,7 @@ exportMethods(sigTests)
 exportMethods(summary)
 exportMethods(weights)
 importFrom(Rcpp,sourceCpp)
+importFrom(checkmate,assert_choice)
 importFrom(checkmate,assert_class)
 importFrom(checkmate,assert_int)
 importFrom(checkmate,assert_logical)

NEWS.md

@@ -1,3 +1,11 @@
+# melt 1.11.3
+## NEW FEATURES
+* New `"ael"` option has been added in the `calibrate` argument of `elt()` for adjusted empirical likelihood calibration.
+
+## MINOR IMPROVEMENTS
+* The package vignette has been updated.
+
+
 # melt 1.11.2
 ## MINOR IMPROVEMENTS
 * The package vignette has been updated.

R/AllClasses.R

@@ -24,6 +24,8 @@
 #' @slot nthreads A single integer for the number of threads for parallel
 #'   computation via OpenMP (if available).
 #' @slot seed A single integer for the seed for random number generation.
+#' @slot an A single numeric representing the scaling factor for adjusted
+#'   empirical likelihood calibration.
 #' @slot b A single integer for the number of bootstrap replicates.
 #' @slot m A single integer for the number of Monte Carlo samples.
 #' @aliases ControlEL
@@ -33,7 +35,7 @@ setClass("ControlEL",
   slots = c(
     maxit = "integer", maxit_l = "integer", tol = "numeric", tol_l = "numeric",
     step = "ANY", th = "ANY", verbose = "logical", keep_data = "logical",
-    nthreads = "integer", seed = "ANY", b = "integer", m = "integer"
+    nthreads = "integer", seed = "ANY", an = "ANY", b = "integer", m = "integer"
   )
 )
 

R/AllGenerics.R

@@ -283,7 +283,7 @@ setGeneric("elmt", function(object,
 
 #' Empirical likelihood test
 #'
-#' Tests a linear hypothesis.
+#' Tests a linear hypothesis with various calibration options.
 #'
 #' @param object An object that inherits from \linkS4class{EL}.
 #' @param rhs A numeric vector or a column matrix for the right-hand side of
@@ -295,9 +295,9 @@ setGeneric("elmt", function(object,
 #'   of parameters. Or a character vector with a symbolic description of the
 #'   hypothesis is allowed. Defaults to `NULL`. See ‘Details’.
 #' @param alpha A single numeric for the significance level. Defaults to `0.05`.
-#' @param calibrate A single character for the calibration method. It is
-#'   case-insensitive and must be one of `"chisq"`, `"boot"`, or `"f"`.
-#'   Defaults to `"chisq"`. See ‘Details’.
+#' @param calibrate A single character representing the calibration method. It
+#'   is case-insensitive and must be one of `"ael"`, `"boot"`, `"chisq"`, or
+#'   `"f"`. Defaults to `"chisq"`. See ‘Details’.
 #' @param control An object of class \linkS4class{ControlEL} constructed by
 #'   [el_control()]. Defaults to `NULL` and inherits the `control` slot in
 #'   `object`.
@@ -318,12 +318,13 @@ setGeneric("elmt", function(object,
 #'   1. If `lhs` is `NULL`, \eqn{L} is set to the identity matrix and the
 #'   problem reduces to evaluating at \eqn{r} as \eqn{l(r)}.
 #'
-#'   `calibrate` specifies the calibration method used. Three methods are
-#'   available: `"chisq"` (chi-square calibration), `"boot"` (bootstrap
-#'   calibration), and `"f"` (\eqn{F} calibration). `"boot"` is applicable only
-#'   when `lhs` is `NULL`. The `nthreads`, `seed`, and `B` slots in `control`
-#'   apply to the bootstrap procedure. `"f"` is applicable only to the mean
-#'   parameter when `lhs` is `NULL`.
+#'   `calibrate` specifies the calibration method used. Four methods are
+#'   available: `"ael"` (adjusted empirical likelihood calibration), `"boot"`
+#'   (bootstrap calibration), `"chisq"` (chi-square calibration), and `"f"`
+#'   (\eqn{F} calibration). When `lhs` is not `NULL`, only `"chisq"` is
+#'   available. `"f"` is applicable only to the mean parameter. The `an` slot in
+#'   `control` applies specifically to `"ael"`, while the `nthreads`, `seed`,
+#'   and `B` slots apply to `"boot"`.
 #' @return An object of class of \linkS4class{ELT}. If `lhs` is non-`NULL`, the
 #'   `optim` slot corresponds to that of \linkS4class{CEL}. Otherwise, it
 #'   corresponds to that of \linkS4class{EL}.
@@ -333,6 +334,11 @@ setGeneric("elmt", function(object,
 #'   \emph{Statistical Methods & Applications}, **19**(4), 463--476.
 #'   \doi{10.1007/s10260-010-0137-9}.
 #' @references
+#'   Chen J, Variyath AM, Abraham B (2008).
+#'   ``Adjusted Empirical Likelihood and Its Properties.''
+#'   \emph{Journal of Computational and Graphical Statistics}, **17**(2),
+#'   426--443.
+#' @references
 #'   Kim E, MacEachern SN, Peruggia M (2024).
 #'   ``melt: Multiple Empirical Likelihood Tests in R.''
 #'   \emph{Journal of Statistical Software}, **108**(5), 1--33.
@@ -345,9 +351,15 @@ setGeneric("elmt", function(object,
 #' @seealso \linkS4class{EL}, \linkS4class{ELT}, [elmt()], [el_control()]
 #' @usage NULL
 #' @examples
-#' ## F calibration for the mean
+#' ## Adjusted empirical likelihood calibration
 #' data("precip")
 #' fit <- el_mean(precip, 32)
+#' elt(fit, rhs = 100, calibrate = "ael")
+#'
+#' ## Bootstrap calibration
+#' elt(fit, rhs = 32, calibrate = "boot")
+#'
+#' ## F calibration
 #' elt(fit, rhs = 32, calibrate = "f")
 #'
 #' ## Test of no treatment effect

R/calibrate.R

@@ -14,7 +14,7 @@
 #' @noRd
 calibrate <- function(calibrate, alpha, statistic, p, par, object, control) {
   switch(calibrate,
-    "chisq" = {
+    "ael" = {
       c(
         cv = qchisq(1 - alpha, df = p),
         pval = pchisq(statistic, df = p, lower.tail = FALSE)
@@ -32,6 +32,12 @@ calibrate <- function(calibrate, alpha, statistic, p, par, object, control) {
         control@maxit_l, control@tol_l, control@th, getWeights(object)
       )
     },
+    "chisq" = {
+      c(
+        cv = qchisq(1 - alpha, df = p),
+        pval = pchisq(statistic, df = p, lower.tail = FALSE)
+      )
+    },
     "f" = {
       n <- nrow(getData(object))
       c(

R/el_control.R

@@ -41,6 +41,9 @@
 #'   non-overlapping but reproducible sequence of random numbers. The
 #'   Xoshiro256+ pseudo-random number generator is used internally to work with
 #'   OpenMP.
+#' @param an A single numeric representing the scaling factor for adjusted
+#'   empirical likelihood calibration. It only applies to [elt()] when
+#'   `calibrate` is set to `"ael"`. Defaults to `NULL`.
 #' @param b A single integer for the number of bootstrap replicates. It only
 #'   applies to [elt()] when `calibrate` is set to `"boot"`. Defaults to
 #'   `10000`.
@@ -61,6 +64,7 @@ el_control <- function(maxit = 200L,
                        keep_data = TRUE,
                        nthreads,
                        seed = NULL,
+                       an = NULL,
                        b = 10000L,
                        m = 1e+06L) {
   maxit <- assert_int(maxit, lower = 1L, coerce = TRUE)
@@ -92,11 +96,14 @@ el_control <- function(maxit = 200L,
   if (isFALSE(is.null(seed))) {
     seed <- assert_int(seed, coerce = TRUE)
   }
+  if (isFALSE(is.null(an))) {
+    an <- assert_number(an, lower = .Machine$double.eps, finite = TRUE)
+  }
   b <- assert_int(b, lower = 1L, coerce = TRUE)
   m <- assert_int(m, lower = 1L, coerce = TRUE)
   new("ControlEL",
     maxit = maxit, maxit_l = maxit_l, tol = tol, tol_l = tol_l, step = step,
     th = th, verbose = verbose, keep_data = keep_data, nthreads = nthreads,
-    seed = seed, b = b, m = m
+    seed = seed, an = an, b = b, m = m
   )
 }

R/elt-methods.R

@@ -18,15 +18,16 @@ setMethod("elt", "EL", function(object,
   onames <- names(logProb(object))
   pnames <- names(getOptim(object)$par)
   h <- validate_hypothesis(rhs, lhs, getNumPar(object), pnames)
-  alpha <- assert_number(alpha, lower = 0, upper = 1, finite = TRUE)
-  calibrate <- validate_calibrate(calibrate)
+  assert_number(alpha, lower = 0, upper = 1, finite = TRUE)
+  calibrate <- assert_choice(tolower(calibrate), c("ael", "boot", "chisq", "f"))
   method <- getMethodEL(object)
   maxit <- control@maxit
   maxit_l <- control@maxit_l
   tol <- control@tol
   tol_l <- control@tol_l
   step <- control@step
   th <- control@th
+  an <- control@an
   w <- getWeights(object)
   if (is.null(lhs)) {
     stopifnot(
@@ -35,6 +36,17 @@ setMethod("elt", "EL", function(object,
     )
     par <- h$r
     out <- compute_EL(method, par, getData(object), maxit_l, tol_l, th, w)
+    if (isTRUE(calibrate == "ael")) {
+      if (is.null(an)) {
+        n <- nobs(object)
+        an <- log(n) / 2
+      }
+      g <- out$g
+      pg <- -an * colMeans(g)
+      g <- rbind(g, pg)
+      out <- compute_generic_EL(g, maxit_l, tol_l, th, w)
+      out$optim <- c(par = list(par), out$optim)
+    }
     optim <- validate_optim(out$optim)
     names(optim$par) <- pnames
     optim$cstr <- FALSE
@@ -50,9 +62,11 @@ setMethod("elt", "EL", function(object,
   }
   # Proceed with chi-square calibration for non-NULL `lhs`
   stopifnot(
-    "Bootstrap calibration is applicable only when `lhs` is NULL." =
+    "AEL calibration is applicable only when `lhs` is `NULL`." =
+      (calibrate != "ael"),
+    "Bootstrap calibration is applicable only when `lhs` is `NULL`." =
       (calibrate != "boot"),
-    "F calibration is applicable only when `lhs` is NULL." =
+    "F calibration is applicable only when `lhs` is `NULL`." =
       (calibrate != "f")
   )
   out <- test_hypothesis(
@@ -95,15 +109,16 @@ setMethod("elt", "QGLM", function(object,
   nm <- names(getOptim(object)$par)
   pnames <- nm[-getNumPar(object)]
   h <- validate_hypothesis(rhs, lhs, getNumPar(object) - 1L, pnames)
-  alpha <- assert_number(alpha, lower = 0, upper = 1, finite = TRUE)
-  calibrate <- validate_calibrate(calibrate)
+  assert_number(alpha, lower = 0, upper = 1, finite = TRUE)
+  assert_choice(tolower(calibrate), c("ael", "boot", "chisq", "f"))
   method <- getMethodEL(object)
   maxit <- control@maxit
   maxit_l <- control@maxit_l
   tol <- control@tol
   tol_l <- control@tol_l
   step <- control@step
   th <- control@th
+  an <- control@an
   w <- getWeights(object)
   if (is.null(lhs)) {
     stopifnot(
@@ -115,6 +130,17 @@ setMethod("elt", "QGLM", function(object,
       par <- c(h$r, object@dispersion)
     }
     out <- compute_EL(method, par, getData(object), maxit_l, tol_l, th, w)
+    if (isTRUE(calibrate == "ael")) {
+      if (is.null(an)) {
+        n <- nobs(object)
+        an <- log(n) / 2
+      }
+      g <- out$g
+      pg <- -an * colMeans(g)
+      g <- rbind(g, pg)
+      out <- compute_generic_EL(g, maxit_l, tol_l, th, w)
+      out$optim <- c(par = list(par), out$optim)
+    }
     optim <- validate_optim(out$optim)
     names(optim$par) <- nm
     optim$cstr <- TRUE
@@ -124,12 +150,14 @@ setMethod("elt", "QGLM", function(object,
     return(new("ELT",
       optim = optim, logp = setNames(out$logp, onames), logl = out$logl,
       loglr = out$loglr, statistic = out$statistic, df = length(par),
-      pval = unname(cal["pval"]), cv = unname(cal["cv"]), rhs = h$r, lhs = h$l,
+      pval = unname(cal["pval"]), cv = unname(cal["cv"]), rhs = par, lhs = h$l,
       alpha = alpha, calibrate = calibrate, control = control
     ))
   }
   stopifnot(
-    "Bootstrap calibration is applicable only when `lhs` is NULL." =
+    "AEL calibration is applicable only when `lhs` is `NULL`." =
+      (calibrate != "ael"),
+    "Bootstrap calibration is applicable only when `lhs` is `NULL`." =
       (calibrate != "boot"),
     "F calibration is applicable only to the mean." = (calibrate != "f")
   )
@@ -172,11 +200,12 @@ setMethod("elt", "SD", function(object,
   onames <- names(logProb(object))
   pnames <- names(getOptim(object)$par)
   h <- validate_hypothesis(rhs, lhs, getNumPar(object), pnames)
-  alpha <- assert_number(alpha, lower = 0, upper = 1, finite = TRUE)
-  calibrate <- validate_calibrate(calibrate)
+  assert_number(alpha, lower = 0, upper = 1, finite = TRUE)
+  assert_choice(tolower(calibrate), c("ael", "boot", "chisq", "f"))
   maxit_l <- control@maxit_l
   tol_l <- control@tol_l
   th <- control@th
+  an <- control@an
   w <- getWeights(object)
   if (is.null(lhs)) {
     par <- h$r
@@ -188,6 +217,17 @@ setMethod("elt", "SD", function(object,
         (par >= .Machine$double.eps)
     )
     out <- compute_EL("sd", par, getData(object), maxit_l, tol_l, th, w)
+    if (isTRUE(calibrate == "ael")) {
+      if (is.null(an)) {
+        n <- nobs(object)
+        an <- log(n) / 2
+      }
+      g <- out$g
+      pg <- -an * colMeans(g)
+      g <- rbind(g, pg)
+      out <- compute_generic_EL(g, maxit_l, tol_l, th, w)
+      out$optim <- c(par = list(par), out$optim)
+    }
     optim <- validate_optim(out$optim)
     names(optim$par) <- pnames
     optim$cstr <- FALSE
@@ -200,6 +240,8 @@ setMethod("elt", "SD", function(object,
     ))
   }
   stopifnot(
+    "AEL calibration is applicable only when `lhs` is `NULL`." =
+      (calibrate != "ael"),
     "Bootstrap calibration is applicable only when `lhs` is `NULL`." =
       (calibrate != "boot"),
     "F calibration is applicable only when `lhs` is `NULL`." =
@@ -233,8 +275,8 @@ setMethod("elt", "missing", function(object,
                                      alpha = 0.05,
                                      calibrate = "chisq",
                                      control = NULL) {
-  alpha <- assert_number(alpha, lower = 0, upper = 1, finite = TRUE)
-  calibrate <- validate_calibrate(calibrate)
+  assert_number(alpha, lower = 0, upper = 1, finite = TRUE)
+  assert_choice(tolower(calibrate), c("ael", "boot", "chisq", "f"))
   if (is.null(control)) {
     control <- el_control()
   } else {

R/melt-package.R

@@ -1,4 +1,5 @@
 ## usethis namespace: start
+#' @importFrom checkmate assert_choice
 #' @importFrom checkmate assert_class
 #' @importFrom checkmate assert_int
 #' @importFrom checkmate assert_logical

R/print-methods.R

@@ -57,8 +57,9 @@ setMethod("print", "ELT", function(x,
     x@rhs, x@lhs, names(getOptim(x)$par)[seq_len(ncol(x@lhs))], digits
   ), sep = "\n")
   method <- switch(x@calibrate,
-    "chisq" = "Chi-square",
+    "ael" = "Adjusted EL",
     "boot" = "Bootstrap",
+    "chisq" = "Chi-square",
     "f" = "F"
   )
   cat(paste0(
@@ -209,8 +210,9 @@ setMethod(
       x@rhs, x@lhs, names(getOptim(x)$par)[seq_len(ncol(x@lhs))], digits
     ), sep = "\n")
     method <- switch(x@calibrate,
-      "chisq" = "Chi-square",
+      "ael" = "Adjusted EL",
       "boot" = "Bootstrap",
+      "chisq" = "Chi-square",
       "f" = "F"
     )
     cat(paste0(
@@ -265,7 +267,7 @@ setMethod(
       cat("\nLagrange multipliers:\n")
       print.default(getOptim(x)$lambda, digits = digits, ...)
       cat("\nMaximum EL estimates:\n")
-      print.default(coef(x)[,1L], digits = digits, ...)
+      print.default(coef(x)[, 1L], digits = digits, ...)
       cat(paste(
         "\nlogL:", format.default(logL(x), digits = digits, ...),
         ", logLR:", format.default(logLR(x), digits = digits, ...)
@@ -337,7 +339,7 @@ setMethod(
       cat("\nLagrange multipliers:\n")
       print.default(getOptim(x)$lambda, digits = digits, ...)
       cat("\nMaximum EL estimates:\n")
-      print.default(coef(x)[,1L], digits = digits, ...)
+      print.default(coef(x)[, 1L], digits = digits, ...)
       cat(paste(
         "\nlogL:", format.default(logL(x), digits = digits, ...),
         ", logLR:", format.default(logLR(x), digits = digits, ...)