diff --git a/.github/workflows/test-coverage.yaml b/.github/workflows/test-coverage.yaml index 960234cd..258c0087 100644 --- a/.github/workflows/test-coverage.yaml +++ b/.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] diff --git a/DESCRIPTION b/DESCRIPTION index f023615a..1811680b 100644 --- a/DESCRIPTION +++ b/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")), diff --git a/NAMESPACE b/NAMESPACE index 12da1df7..38016fd2 100644 --- a/NAMESPACE +++ b/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) diff --git a/NEWS.md b/NEWS.md index b5f6efff..1c2aafe1 100644 --- a/NEWS.md +++ b/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. diff --git a/R/AllClasses.R b/R/AllClasses.R index 9bcc5311..09ee2bcf 100644 --- a/R/AllClasses.R +++ b/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" ) ) diff --git a/R/AllGenerics.R b/R/AllGenerics.R index e405196b..bc251ac3 100644 --- a/R/AllGenerics.R +++ b/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 diff --git a/R/calibrate.R b/R/calibrate.R index e8f844cf..73339873 100644 --- a/R/calibrate.R +++ b/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( diff --git a/R/el_control.R b/R/el_control.R index 98ae6214..5fa4b32a 100644 --- a/R/el_control.R +++ b/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 ) } diff --git a/R/elt-methods.R b/R/elt-methods.R index 5115ced4..d4605237 100644 --- a/R/elt-methods.R +++ b/R/elt-methods.R @@ -18,8 +18,8 @@ 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 @@ -27,6 +27,7 @@ setMethod("elt", "EL", function(object, 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,8 +109,8 @@ 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 @@ -104,6 +118,7 @@ setMethod("elt", "QGLM", function(object, 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 { diff --git a/R/melt-package.R b/R/melt-package.R index fec3dcfe..fa4cbfd6 100644 --- a/R/melt-package.R +++ b/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 diff --git a/R/print-methods.R b/R/print-methods.R index 915a0fd2..513a4cd6 100644 --- a/R/print-methods.R +++ b/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, ...) diff --git a/R/validate.R b/R/validate.R index c94d74a6..fc681ee5 100644 --- a/R/validate.R +++ b/R/validate.R @@ -1,23 +1,3 @@ -#' Validate `calibrate` -#' -#' Validate `calibrate` in [elt()]. -#' -#' @param calibrate A single character. -#' @return A single character. -#' @noRd -validate_calibrate <- function(calibrate) { - assert_string(calibrate) - table <- c("chisq", "boot", "f") - calibrate <- table[pmatch(tolower(calibrate), table = table)] - if (isTRUE(is.na(calibrate))) { - stop(gettextf( - "`calibrate` must be one of %s, %s, or %s.", - dQuote("chisq"), dQuote("boot"), dQuote("f") - ), domain = NA) - } - calibrate -} - #' Validate `cv` #' #' Validate `cv` in [confint()] and [confreg()]. diff --git a/README.Rmd b/README.Rmd index 5b853ee4..b89610c1 100644 --- a/README.Rmd +++ b/README.Rmd @@ -20,7 +20,7 @@ knitr::opts_chunk$set( [![Project Status: Active - The project has reached a stable, usable state and is being actively developed.](https://www.repostatus.org/badges/latest/active.svg)](https://www.repostatus.org/#active) [![R-CMD-check](https://github.com/ropensci/melt/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/ropensci/melt/actions/workflows/R-CMD-check.yaml) [![pkgcheck](https://github.com/ropensci/melt/workflows/pkgcheck/badge.svg)](https://github.com/ropensci/melt/actions?query=workflow%3Apkgcheck) -[![Codecov test coverage](https://codecov.io/gh/ropensci/melt/branch/master/graph/badge.svg)](https://app.codecov.io/gh/ropensci/melt?branch=master) +[![Codecov test coverage](https://codecov.io/gh/ropensci/melt/branch/main/graph/badge.svg)](https://app.codecov.io/gh/ropensci/melt?branch=main) [![CRAN status](https://www.r-pkg.org/badges/version/melt)](https://CRAN.R-project.org/package=melt) [![runiverse](https://ropensci.r-universe.dev/badges/melt)](http://ropensci.r-universe.dev/ui/#package:melt) [![ropensci review](https://badges.ropensci.org/550_status.svg)](https://github.com/ropensci/software-review/issues/550) @@ -67,7 +67,7 @@ melt provides an intuitive API for performing the most common data analysis task * `el_glm()` fits a generalized linear model with empirical likelihood. * `confint()` computes confidence intervals for model parameters. * `confreg()` computes confidence region for model parameters. -* `elt()` tests a linear hypothesis. +* `elt()` tests a hypothesis with various calibration options. * `elmt()` performs multiple testing simultaneously. @@ -78,7 +78,19 @@ set.seed(971112) ## Test for the mean data("precip") -el_mean(precip, par = 30) +(fit <- el_mean(precip, par = 30)) + + +## Adjusted empirical likelihood calibration +elt(fit, rhs = 30, calibrate = "ael") + + +## Bootstrap calibration +elt(fit, rhs = 30, calibrate = "boot") + + +## F calibration +elt(fit, rhs = 30, calibrate = "f") ## Linear model diff --git a/README.md b/README.md index 858d170f..76cb294b 100644 --- a/README.md +++ b/README.md @@ -11,7 +11,7 @@ developed.](https://www.repostatus.org/badges/latest/active.svg)](https://www.re [![R-CMD-check](https://github.com/ropensci/melt/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/ropensci/melt/actions/workflows/R-CMD-check.yaml) [![pkgcheck](https://github.com/ropensci/melt/workflows/pkgcheck/badge.svg)](https://github.com/ropensci/melt/actions?query=workflow%3Apkgcheck) [![Codecov test -coverage](https://codecov.io/gh/ropensci/melt/branch/master/graph/badge.svg)](https://app.codecov.io/gh/ropensci/melt?branch=master) +coverage](https://codecov.io/gh/ropensci/melt/branch/main/graph/badge.svg)](https://app.codecov.io/gh/ropensci/melt?branch=main) [![CRAN status](https://www.r-pkg.org/badges/version/melt)](https://CRAN.R-project.org/package=melt) [![runiverse](https://ropensci.r-universe.dev/badges/melt)](http://ropensci.r-universe.dev/ui/#package:melt) @@ -69,7 +69,7 @@ analysis tasks: - `el_glm()` fits a generalized linear model with empirical likelihood. - `confint()` computes confidence intervals for model parameters. - `confreg()` computes confidence region for model parameters. -- `elt()` tests a linear hypothesis. +- `elt()` tests a hypothesis with various calibration options. - `elmt()` performs multiple testing simultaneously. ## Usage @@ -80,7 +80,7 @@ set.seed(971112) ## Test for the mean data("precip") -el_mean(precip, par = 30) +(fit <- el_mean(precip, par = 30)) #> #> Empirical Likelihood #> @@ -93,6 +93,51 @@ el_mean(precip, par = 30) #> EL evaluation: converged +## Adjusted empirical likelihood calibration +elt(fit, rhs = 30, calibrate = "ael") +#> +#> Empirical Likelihood Test +#> +#> Hypothesis: +#> par = 30 +#> +#> Significance level: 0.05, Calibration: Adjusted EL +#> +#> Statistic: 7.744, Critical value: 3.841 +#> p-value: 0.005389 +#> EL evaluation: converged + + +## Bootstrap calibration +elt(fit, rhs = 30, calibrate = "boot") +#> +#> Empirical Likelihood Test +#> +#> Hypothesis: +#> par = 30 +#> +#> Significance level: 0.05, Calibration: Bootstrap +#> +#> Statistic: 8.285, Critical value: 3.84 +#> p-value: 0.0041 +#> EL evaluation: converged + + +## F calibration +elt(fit, rhs = 30, calibrate = "f") +#> +#> Empirical Likelihood Test +#> +#> Hypothesis: +#> par = 30 +#> +#> Significance level: 0.05, Calibration: F +#> +#> Statistic: 8.285, Critical value: 3.98 +#> p-value: 0.005318 +#> EL evaluation: converged + + ## Linear model data("mtcars") fit_lm <- el_lm(mpg ~ disp + hp + wt + qsec, data = mtcars) @@ -226,13 +271,13 @@ summary(fit_glm) #> #> Coefficients: #> Estimate Chisq Pr(>Chisq) -#> (Intercept) -0.10977 0.090 0.763757 +#> (Intercept) -0.10977 0.090 0.764 #> log(mass) 0.24750 425.859 < 2e-16 *** -#> fruit 0.04654 11.584 0.000665 *** +#> fruit 0.04654 29.024 7.15e-08 *** #> foliage -19.40632 65.181 6.83e-16 *** #> varGZ -0.25760 17.308 3.18e-05 *** -#> trtSpray 0.06724 0.860 0.353820 -#> trtFurrow -0.03634 0.217 0.641379 +#> trtSpray 0.06724 0.860 0.354 +#> trtFurrow -0.03634 0.217 0.641 #> trtSeed 0.34790 19.271 1.13e-05 *** #> --- #> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 diff --git a/codemeta.json b/codemeta.json index ec6e74ac..416a70f0 100644 --- a/codemeta.json +++ b/codemeta.json @@ -8,7 +8,7 @@ "codeRepository": "https://github.com/ropensci/melt", "issueTracker": "https://github.com/ropensci/melt/issues", "license": "https://spdx.org/licenses/GPL-2.0", - "version": "1.11.2", + "version": "1.11.3", "programmingLanguage": { "@type": "ComputerLanguage", "name": "R", @@ -291,7 +291,7 @@ }, "SystemRequirements": null }, - "fileSize": "300699.916KB", + "fileSize": "367508.966KB", "citation": [ { "@type": "ScholarlyArticle", @@ -334,7 +334,7 @@ } ], "releaseNotes": "https://github.com/ropensci/melt/blob/master/NEWS.md", - "contIntegration": ["https://github.com/ropensci/melt/actions/workflows/R-CMD-check.yaml", "https://github.com/ropensci/melt/actions?query=workflow%3Apkgcheck", "https://app.codecov.io/gh/ropensci/melt?branch=master"], + "contIntegration": ["https://github.com/ropensci/melt/actions/workflows/R-CMD-check.yaml", "https://github.com/ropensci/melt/actions?query=workflow%3Apkgcheck", "https://app.codecov.io/gh/ropensci/melt?branch=main"], "developmentStatus": "https://www.repostatus.org/#active", "review": { "@type": "Review", diff --git a/inst/WORDLIST b/inst/WORDLIST index c67b68eb..25142869 100644 --- a/inst/WORDLIST +++ b/inst/WORDLIST @@ -39,7 +39,6 @@ LM Lazar MELE Mroz -NA's Neonicotinoid Obregon OpenMP @@ -49,7 +48,6 @@ Pederson Poveda QGLM Qin -Qu Raphson Rcpp RcppEigen @@ -59,35 +57,28 @@ R’ SES SPM SX -Signif Sinica Statistica Thiamethoxam UScrime +Variyath Wilks Xoshiro Zhang Zhao Zhu agricolae -byrow -calpha carData -cd -cex chisq clothianidin coef coercible colMeans -conf confint confreg conv -cr critVal cstr -ctrl cv df doi @@ -111,7 +102,6 @@ glht glm gmm hc -hcyes iteratively jss lfp @@ -123,55 +113,39 @@ logProb logl loglr logp -lwd -lwg -lwr maxit misspecification momentfit multcomp -ncol neonicotinoid nonconvex npoints -nrow nthreads optim optimality pVal parallelization -parm -pch piecewiseSEM pkgVersion pkgcheck priori pval -quasipoisson rOpenSci ridgeline ropensci runiverse -scb sd stars’ straightneck studentization -tck th thiamethoxam tidyverse tol translocates trt -trtFurrow -trtSeed -trtSpray -upr -varGZ vcov wc -wcyes xs ys µg diff --git a/man/ControlEL-class.Rd b/man/ControlEL-class.Rd index b9d45903..5e7f8813 100644 --- a/man/ControlEL-class.Rd +++ b/man/ControlEL-class.Rd @@ -43,6 +43,9 @@ computation via OpenMP (if available).} \item{\code{seed}}{A single integer for the seed for random number generation.} +\item{\code{an}}{A single numeric representing the scaling factor for adjusted +empirical likelihood calibration.} + \item{\code{b}}{A single integer for the number of bootstrap replicates.} \item{\code{m}}{A single integer for the number of Monte Carlo samples.} diff --git a/man/el_control.Rd b/man/el_control.Rd index b9ca6712..c076619f 100644 --- a/man/el_control.Rd +++ b/man/el_control.Rd @@ -15,6 +15,7 @@ el_control( keep_data = TRUE, nthreads, seed = NULL, + an = NULL, b = 10000L, m = 1000000L ) @@ -69,6 +70,10 @@ non-overlapping but reproducible sequence of random numbers. The Xoshiro256+ pseudo-random number generator is used internally to work with OpenMP.} +\item{an}{A single numeric representing the scaling factor for adjusted +empirical likelihood calibration. It only applies to \code{\link[=elt]{elt()}} when +\code{calibrate} is set to \code{"ael"}. Defaults to \code{NULL}.} + \item{b}{A single integer for the number of bootstrap replicates. It only applies to \code{\link[=elt]{elt()}} when \code{calibrate} is set to \code{"boot"}. Defaults to \code{10000}.} diff --git a/man/elt.Rd b/man/elt.Rd index e7516df5..55685c5c 100644 --- a/man/elt.Rd +++ b/man/elt.Rd @@ -32,9 +32,9 @@ hypothesis is allowed. Defaults to \code{NULL}. See ‘Details’.} \item{alpha}{A single numeric for the significance level. Defaults to \code{0.05}.} -\item{calibrate}{A single character for the calibration method. It is -case-insensitive and must be one of \code{"chisq"}, \code{"boot"}, or \code{"f"}. -Defaults to \code{"chisq"}. See ‘Details’.} +\item{calibrate}{A single character representing the calibration method. It +is case-insensitive and must be one of \code{"ael"}, \code{"boot"}, \code{"chisq"}, or +\code{"f"}. Defaults to \code{"chisq"}. See ‘Details’.} \item{control}{An object of class \linkS4class{ControlEL} constructed by \code{\link[=el_control]{el_control()}}. Defaults to \code{NULL} and inherits the \code{control} slot in @@ -46,7 +46,7 @@ An object of class of \linkS4class{ELT}. If \code{lhs} is non-\code{NULL}, the corresponds to that of \linkS4class{EL}. } \description{ -Tests a linear hypothesis. +Tests a linear hypothesis with various calibration options. } \details{ \code{\link[=elt]{elt()}} performs the constrained minimization of \eqn{l(\theta)} @@ -68,17 +68,24 @@ is performed with the right-hand side \eqn{r} and the left-hand side problem reduces to evaluating at \eqn{r} as \eqn{l(r)}. } -\code{calibrate} specifies the calibration method used. Three methods are -available: \code{"chisq"} (chi-square calibration), \code{"boot"} (bootstrap -calibration), and \code{"f"} (\eqn{F} calibration). \code{"boot"} is applicable only -when \code{lhs} is \code{NULL}. The \code{nthreads}, \code{seed}, and \code{B} slots in \code{control} -apply to the bootstrap procedure. \code{"f"} is applicable only to the mean -parameter when \code{lhs} is \code{NULL}. +\code{calibrate} specifies the calibration method used. Four methods are +available: \code{"ael"} (adjusted empirical likelihood calibration), \code{"boot"} +(bootstrap calibration), \code{"chisq"} (chi-square calibration), and \code{"f"} +(\eqn{F} calibration). When \code{lhs} is not \code{NULL}, only \code{"chisq"} is +available. \code{"f"} is applicable only to the mean parameter. The \code{an} slot in +\code{control} applies specifically to \code{"ael"}, while the \code{nthreads}, \code{seed}, +and \code{B} slots apply to \code{"boot"}. } \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 @@ -104,6 +111,11 @@ Adimari G, Guolo A (2010). \emph{Statistical Methods & Applications}, \strong{19}(4), 463--476. \doi{10.1007/s10260-010-0137-9}. +Chen J, Variyath AM, Abraham B (2008). +``Adjusted Empirical Likelihood and Its Properties.'' +\emph{Journal of Computational and Graphical Statistics}, \strong{17}(2), +426--443. + Kim E, MacEachern SN, Peruggia M (2024). ``melt: Multiple Empirical Likelihood Tests in R.'' \emph{Journal of Statistical Software}, \strong{108}(5), 1--33. diff --git a/src/EL.cpp b/src/EL.cpp index 4ca055f0..4ad540e8 100644 --- a/src/EL.cpp +++ b/src/EL.cpp @@ -18,17 +18,18 @@ EL::EL(const std::string method, const double tol_l, const double th, const Eigen::Ref &wt) - : par{par0}, - l{Eigen::VectorXd::Zero(par0.size())}, - mele_fn{set_mele_fn(method)}, - w{wt}, - maxit_l{maxit_l}, + : maxit_l{maxit_l}, tol_l{tol_l}, th{th}, n{static_cast(x.rows())}, - g_fn{EL::set_g_fn(method)} + g_fn{EL::set_g_fn(method)}, + par{par0}, + l{Eigen::VectorXd::Zero(par0.size())}, + mele_fn{set_mele_fn(method)}, + w{wt}, + g{g_fn(x, par)} { - set_el(g_fn(x, par), wt); + set_el(g, wt); } EL::EL(const Eigen::Ref &g, @@ -36,15 +37,16 @@ EL::EL(const Eigen::Ref &g, const double tol_l, const double th, const Eigen::Ref &wt) - : par{}, - l{Eigen::VectorXd::Zero(g.cols())}, - mele_fn{}, - w{wt}, - maxit_l{maxit_l}, + : maxit_l{maxit_l}, tol_l{tol_l}, th{th}, n{static_cast(g.rows())}, - g_fn{} + g_fn{}, + par{}, + l{Eigen::VectorXd::Zero(g.cols())}, + mele_fn{}, + w{wt}, + g{g} { set_el(g, wt); } diff --git a/src/EL.h b/src/EL.h index e961ed3e..edfe05d9 100644 --- a/src/EL.h +++ b/src/EL.h @@ -7,6 +7,18 @@ class EL { +private: + // members + const int maxit_l; // maximum number of iterations + const double tol_l; // relative convergence tolerance + const double th; // threshold value for negative log-likelihood ratio + const int n; // sample size + // estimating function + const std::function &, + const Eigen::Ref &)> + g_fn; + public: // members const Eigen::VectorXd par; // parameter value specified @@ -18,6 +30,7 @@ class EL int iter{0}; // iterations performed in optimization bool conv{false}; // convergence status const Eigen::ArrayXd w; // weights + const Eigen::MatrixXd g; // constructors EL(const std::string method, @@ -50,17 +63,7 @@ class EL // log-likelihood double loglik() const; -private: - // members - const int maxit_l; // maximum number of iterations - const double tol_l; // relative convergence tolerance - const double th; // threshold value for negative log-likelihood ratio - const int n; // sample size - // estimating function - const std::function &, - const Eigen::Ref &)> - g_fn; + }; class CEL diff --git a/src/compute_EL.cpp b/src/compute_EL.cpp index 57d9cdb0..2bca2718 100644 --- a/src/compute_EL.cpp +++ b/src/compute_EL.cpp @@ -21,6 +21,7 @@ Rcpp::List compute_EL(const std::string method, Rcpp::Named("lambda") = el.l, Rcpp::Named("iterations") = el.iter, Rcpp::Named("convergence") = el.conv), + Rcpp::Named("g") = el.g, Rcpp::Named("logp") = el.logp(x), Rcpp::Named("logl") = el.loglik(), Rcpp::Named("loglr") = -el.nllr, diff --git a/tests/testthat/test-el_control.R b/tests/testthat/test-el_control.R index 667d0f65..29061ef3 100644 --- a/tests/testthat/test-el_control.R +++ b/tests/testthat/test-el_control.R @@ -16,4 +16,5 @@ test_that("Invalid `control`.", { expect_error(el_control(nthreads = c(10, 20))) expect_error(el_control(nthreads = 0)) expect_warning(el_control(nthreads = .Machine$integer.max)) + expect_error(el_control(an = 0)) }) diff --git a/tests/testthat/test-elt.R b/tests/testthat/test-elt.R index d693007d..a899378a 100644 --- a/tests/testthat/test-elt.R +++ b/tests/testthat/test-elt.R @@ -72,6 +72,9 @@ test_that("When elt == evaluation.", { rhs <- c(0, 0) wfit4 <- elt(wfit3, rhs = rhs, lhs = lhs) expect_equal(getOptim(wfit3)$lambda, getOptim(wfit4)$lambda) + expect_error(elt(fit, lhs = c("par"), rhs = 10, calibrate = "ael")) + out <- elt(fit, rhs = 10, calibrate = "ael") + expect_s4_class(out, "ELT") }) test_that("`conv()` method and calibration.", { @@ -81,8 +84,8 @@ test_that("`conv()` method and calibration.", { test_that("Probabilities add up to 1.", { fit <- el_mean(precip, par = 60) - elt <- elt(fit, rhs = 65) - expect_equal(sum(exp(logProb(elt))), 1, tolerance = 1e-07) + out <- elt(fit, rhs = 65) + expect_equal(sum(exp(logProb(out))), 1, tolerance = 1e-07) }) test_that("Vector `lhs`.", { @@ -130,11 +133,14 @@ test_that("`SD` class.", { expect_error(elt(fit, rhs = -1)) expect_error(elt(fit, rhs = rhs, lhs = lhs, calibrate = "boot")) expect_error(elt(fit, rhs = rhs, lhs = lhs, calibrate = "f")) + expect_error(elt(fit, rhs = -1, lhs = 1, calibrate = "ael")) expect_error(elt(fit, rhs = -1, lhs = 1, calibrate = "f")) out <- elt(fit, rhs = 1) out2 <- elt(fit, rhs = 1, lhs = 2) + out3 <- elt(fit, rhs = 1, calibrate = "ael") expect_s4_class(out, "ELT") expect_s4_class(out2, "ELT") + expect_s4_class(out3, "ELT") }) test_that("`QGLM` class.", { @@ -142,7 +148,9 @@ test_that("`QGLM` class.", { family = quasipoisson("log"), data = mtcars ) out <- elt(fit, rhs = coef(fit)) - expect_s4_class(out, "ELT") out2 <- elt(fit, lhs = c("mpg + cyl")) + out3 <- elt(fit, rhs = c(11, 0, 7), calibrate = "ael") + expect_s4_class(out, "ELT") expect_s4_class(out2, "ELT") + expect_s4_class(out3, "ELT") }) diff --git a/vignettes/article.Rnw b/vignettes/article.Rnw index c127cd66..251eddbb 100644 --- a/vignettes/article.Rnw +++ b/vignettes/article.Rnw @@ -769,16 +769,16 @@ specify the problem in Equation~\ref{eq:constraint}. % \begin{itemize} \item \fct{elt}: Tests a linear hypothesis with EL. It returns an object of - class \class{ELT} that contains the test statistic, the critical value, and - the level of the test. -Several calibration options discussed in Section~\ref{sec:2.2} are available, - and the \(p\)~value is computed by the calibration method chosen. +class \class{ELT} that contains the test statistic, the critical value, and the +level of the test. Various calibration options are available, including some of +those discussed in Section~\ref{sec:2.2} and the adjusted empirical likelihood +calibration method proposed by \citet{chen2008adjusted}. The \(p\)~value is +computed based on the chosen calibration method. \item \fct{elmt}: Tests multiple linear hypotheses simultaneously with EL. Each test can be considered as one instance of \fct{elt}, where the marginal - test statistic has an asymptotic chi-square distribution under the - corresponding null hypothesis. -It returns an object of class \class{ELMT} with slots similar to those in - \class{ELT}. +test statistic has an asymptotic chi-square distribution under the corresponding +null hypothesis. It returns an object of class \class{ELMT} with slots similar +to those in \class{ELT}. \end{itemize} An \class{ELT} object also has the \code{optim} slot, which does not necessarily correspond to the EL optimization. The user can supply an arbitrary parameter @@ -832,7 +832,7 @@ The default values are shown below. \begin{Code} el_control(maxit = 200L, maxit_l = 25L, tol = 1e-06, tol_l = 1e-06, step = NULL, th = NULL, verbose = FALSE, keep_data = TRUE, nthreads, - seed = NULL, b = 10000L, m = 1000000L) + seed = NULL, an = NULL, b = 10000L, m = 1000000L) \end{Code} Specifically, \code{nthreads} specifies the number of threads for parallel computation via \pkg{OpenMP} (if available). By default, it is set to half the diff --git a/vignettes/references.bib b/vignettes/references.bib index dc079e98..e3abd63d 100644 --- a/vignettes/references.bib +++ b/vignettes/references.bib @@ -59,6 +59,17 @@ @Article{chaudhuri2017hamiltonian year = {2017}, } +@Article{chen2008adjusted, + title = {Adjusted Empirical Likelihood and its Properties}, + author = {Jiahua Chen and Asokan Mulayath Variyath and Bovas Abraham}, + journal = {Journal of Computational and Graphical Statistics}, + volume = {17}, + number = {2}, + pages = {426--443}, + year = {2008}, + doi = {10.1198/106186008X321068}, +} + @Article{chen2009effects, author = {Song Xi Chen and Liang Peng and Ying-Li Qin}, title = {Effects of Data Dimension on Empirical Likelihood}, @@ -723,7 +734,7 @@ @Manual{melt title = {\pkg{melt}: Multiple Empirical Likelihood Tests}, author = {Eunseop Kim}, year = {2024}, - note = {\proglang{R} package version 1.11.1}, + note = {\proglang{R} package version 1.11.3}, url = {https://CRAN.R-project.org/package=melt}, }