Skip to content

Commit

Permalink
function parameters format aligned
Browse files Browse the repository at this point in the history
  • Loading branch information
@eunseopkim
eunseopkim committed Feb 2, 2024
1 parent b86a791 commit 37b82f1
Show file tree
Hide file tree
Showing 4 changed files with 91 additions and 52 deletions.
15 changes: 10 additions & 5 deletions R/el_aov.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,11 +3,16 @@
#' @description
#' Fits an one-way analysis of variance model with empirical likelihood.
#'
#' @param formula A formula object. It must specify variables for response and treatment as 'response ~ treatment'.
#' @param data A data frame containing the variables in the formula.
#' @param maxit Maximum number of iterations for optimization. Defaults to 10000.
#' @param abstol Absolute convergence tolerance for optimization. Defaults to 1e-08.
#' @return A list with class \code{c("el_aov", "melt")}.
#' @param formula
#' A formula object. It must specify variables for response and treatment as 'response ~ treatment'.
#' @param data
#' A data frame containing the variables in the formula.
#' @param maxit
#' Maximum number of iterations for optimization. Defaults to 10000.
#' @param abstol
#' Absolute convergence tolerance for optimization. Defaults to 1e-08.
#' @return
#' A list with class \code{c("el_aov", "melt")}.
#' @references
#' Owen, A (1991).
#' "Empirical Likelihood for Linear Models."
Expand Down
33 changes: 22 additions & 11 deletions R/el_pairwise.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,29 +7,40 @@
#' are available: asymptotic Monte Carlo (AMC) and nonparametric bootstrap
#' (NB).
#'
#' @param formula An object of class [`formula`] (or one that can be coerced to
#' @param formula
#' An object of class [`formula`] (or one that can be coerced to
#' that class) for a symbolic description of the model to be fitted. It must
#' specify variables for response, treatment, and block as 'response ~
#' treatment | block'. Note that the use of vertical bar (|) separating
#' treatment and block.
#' @param data A data frame, list or environment (or object coercible by
#' @param data
#' A data frame, list or environment (or object coercible by
#' [as.data.frame()] to a data frame) containing the variables in `formula`.
#' @param control An optional single character that specifies the treatment for
#' @param control
#' An optional single character that specifies the treatment for
#' comparisons with control.
#' @param k A single integer for \eqn{k} in \eqn{k}-FWER. Defaults to 1.
#' @param alpha A single numeric for the overall significance level. Defaults to
#' @param k
#' A single integer for \eqn{k} in \eqn{k}-FWER. Defaults to 1.
#' @param alpha
#' A single numeric for the overall significance level. Defaults to
#' `0.05`.
#' @param method A single character for the procedure to be used; either `AMC`
#' @param method
#' A single character for the procedure to be used; either `AMC`
#' or `NB` is supported. Defaults to `AMC`.
#' @param B A single integer for the number of Monte Carlo samples for the AMC
#' @param B
#' A single integer for the number of Monte Carlo samples for the AMC
#' (number of bootstrap replicates for the NB).
#' @param nthreads A single integer for the number of threads for parallel
#' @param nthreads
#' A single integer for the number of threads for parallel
#' computation via OpenMP (if available). Defaults to `1`.
#' @param maxit A single integer for the maximum number of iterations for
#' @param maxit
#' A single integer for the maximum number of iterations for
#' constrained minimization of empirical likelihood. Defaults to `10000`.
#' @param abstol A single numeric for the the absolute convergence tolerance for
#' @param abstol
#' A single numeric for the the absolute convergence tolerance for
#' optimization. Defaults to `1e-08`.
#' @return A list with class \code{c("pairwise", "melt")}.
#' @return
#' A list with class \code{c("pairwise", "melt")}.
#' @references
#' Kim E, MacEachern SN, Peruggia M (2023).
#' "Empirical likelihood for the analysis of experimental designs."
Expand Down
86 changes: 53 additions & 33 deletions R/el_test.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,15 +27,14 @@
#' ## test for equal means
#' data("clothianidin")
#' el_test(clo ~ trt | blk, clothianidin,
#' lhs = matrix(c(1, -1, 0, 0,
#' 0, 1, -1, 0,
#' 0, 0, 1, -1), byrow = TRUE, nrow = 3))
#' lhs = matrix(c(
#' 1, -1, 0, 0,
#' 0, 1, -1, 0,
#' 0, 0, 1, -1
#' ), byrow = TRUE, nrow = 3)
#' )
#' @export
el_test <- function(formula,
data,
lhs,
rhs = NULL,
maxit = 1e+04,
el_test <- function(formula, data, lhs, rhs = NULL, maxit = 1e+04,
abstol = 1e-08) {
## check formula
f <- attributes(terms(formula))
Expand All @@ -45,14 +44,15 @@ el_test <- function(formula,
length(f$variables) != 3L,
# no other formula
typeof(f$variables[[3L]]) != "language" ||
length(f$variables[[3L]]) != 3L,
length(f$variables[[3L]]) != 3L,
# "|" operator needed
f$variables[[3L]][[1L]] != "|",
# no transformation of variables
typeof(f$variables[[3L]][[2L]]) != "symbol" ||
typeof(f$variables[[3L]][[3L]]) != "symbol",
typeof(f$variables[[3L]][[3L]]) != "symbol",
# distinct variables for treatment and block
f$variables[[3L]][[2L]] == f$variables[[3L]][[3L]])
f$variables[[3L]][[2L]] == f$variables[[3L]][[3L]]
)
) {
stop("invalied model formula. specify formula as 'response ~ treatment | block'")
}
Expand Down Expand Up @@ -87,10 +87,11 @@ el_test <- function(formula,
c <- unclass(table(mf[[3L]], mf[[2L]]))
# model matrix
x <- reshape(mf[order(mf[[2L]]), ],
idvar = names(mf)[3L],
timevar = names(mf)[2L],
v.names = names(mf)[1L],
direction = "wide")
idvar = names(mf)[3L],
timevar = names(mf)[2L],
v.names = names(mf)[1L],
direction = "wide"
)
x <- x[order(x[[names(mf)[3L]]]), ]
# replace NA with 0
x[is.na(x)] <- 0
Expand All @@ -111,7 +112,8 @@ el_test <- function(formula,

## test hypothesis
out <- ELtest(gbd$model_matrix, gbd$incidence_matrix, lhs, rhs,
threshold = nrow(lhs) * 500, maxit, abstol)
threshold = nrow(lhs) * 500, maxit, abstol
)
out$trt <- gbd$trt
out$model.matrix <- gbd$model_matrix
out$incidence.matrix <- gbd$incidence_matrix
Expand Down Expand Up @@ -140,9 +142,10 @@ confint.el_test <- function(object, parm, level = 0.95, ...) {
}
stopifnot(complete.cases(pnames))
if (!missing(level) &&
(length(level) != 1L || !is.finite(level) ||
level < 0 || level > 1))
(length(level) != 1L || !is.finite(level) ||
level < 0 || level > 1)) {
stop("'conf.level' must be a single number between 0 and 1")
}
if (level == 0) {
ci <- matrix(rep(cf, 2L), ncol = 2L)
} else if (level == 1) {
Expand All @@ -151,10 +154,12 @@ confint.el_test <- function(object, parm, level = 0.95, ...) {
} else {
cutoff <- qchisq(level, 1L)
optcfg <- object$optim$control
ci <- EL_confint2(object$data, object$optim$type, cf, cutoff,
optcfg$maxit, optcfg$abstol, optcfg$threshold)
ci <- EL_confint2(
object$data, object$optim$type, cf, cutoff,
optcfg$maxit, optcfg$abstol, optcfg$threshold
)
}
a <- (1 - level)/2
a <- (1 - level) / 2
a <- c(a, 1 - a)
pct <- paste(round(100 * a, 1L), "%")
dimnames(ci) <- list(pnames, pct)
Expand All @@ -167,32 +172,47 @@ print.el_test <- function(x, digits = getOption("digits"), ...) {
cat("Empirical Likelihood Test:", x$optim$type, "\n")
cat("\n")
out <- character()
if (!is.null(x$statistic))
out <- c(out, paste("Chisq", names(x$statistic), "=",
format(x$statistic, digits = max(1L, digits - 2L))))
if (!is.null(x$df))
if (!is.null(x$statistic)) {
out <- c(out, paste(
"Chisq", names(x$statistic), "=",
format(x$statistic, digits = max(1L, digits - 2L))
))
}
if (!is.null(x$df)) {
out <- c(out, paste("df", "=", x$df))
}
if (!is.null(x$p.value)) {
fp <- format.pval(x$p.value, digits = max(1L, digits - 3L))
out <- c(out, paste("p-value", if (startsWith(fp, "<")) fp else paste("=",
fp)))
out <- c(out, paste("p-value", if (startsWith(fp, "<")) {
fp
} else {
paste(
"=",
fp
)
}))
}
cat(strwrap(paste(out, collapse = ", ")), sep = "\n")
if (!is.null(x$alternative)) {
cat("alternative hypothesis: ")
if (!is.null(x$null.value)) {
if (length(x$null.value) == 1L) {
alt.char <- switch(x$alternative, two.sided = "not equal to",
less = "less than", greater = "greater than")
alt.char <- switch(x$alternative,
two.sided = "not equal to",
less = "less than",
greater = "greater than"
)
cat("true ", names(x$null.value), " is ", alt.char,
" ", x$null.value, "\n", sep = "")
}
else {
" ", x$null.value, "\n",
sep = ""
)
} else {
cat(x$alternative, "\nnull values:\n", sep = "")
print(x$null.value, digits = digits, ...)
}
} else {
cat(x$alternative, "\n", sep = "")
}
else cat(x$alternative, "\n", sep = "")
}
if (!is.null(x$coefficients)) {
cat("maximum EL estimates:\n")
Expand Down
9 changes: 6 additions & 3 deletions man/el_test.Rd
View file

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

0 comments on commit 37b82f1

Please sign in to comment.