From 39fa5071c9daec776fa5924072791fa1edebdd80 Mon Sep 17 00:00:00 2001 From: Daan de Jong Date: Thu, 30 Nov 2023 12:53:29 +0100 Subject: [PATCH] new default for p_select --- .Rhistory | 280 +++++++++++++++---------------- NEWS.md | 2 + R/checks_fit.R | 5 +- R/hystar_fit.R | 13 +- R/hystar_sim.R | 3 +- R/plot.R | 30 ++-- R/zzz.R | 8 +- _pkgdown.yml | 1 + man/hystar_fit.Rd | 5 +- tests/testthat/test-checks_fit.R | 4 + 10 files changed, 187 insertions(+), 164 deletions(-) diff --git a/.Rhistory b/.Rhistory index 6e63a07..9ed06d8 100644 --- a/.Rhistory +++ b/.Rhistory @@ -1,143 +1,3 @@ -library(devtools) -install() -devtools::load_all(".") -devtools::document() -?binom -?rbinom -1-250*999/1000 -1-(250*999/1000) -250*999/1000 -1-(250**(999/1000)) -1-(250^(999/1000)) -250^(999/1000) -999/1000 -250^.999 -.999*.999 -1 - .999 ^ 250 -citation(package = "hystar") -qnorm(.975, mean = 0, sd = 1) -usethis::use_pkgdown() -pkgdown::build_site() -packageVersion("hystar") -usethis::use_version() -devtools::load_all(".") -?hystar_fit -z <- z_sim(n_t = 200, n_switches = 5, start_regime = 1) -sim <- hystar_sim(z = z, r = c(-.5, .5), d = 2, phi_R0 = c(0, .6), phi_R1 = 1) -fit <- hystar_fit(sim$data) -devtools::load_all(".") -fit <- hystar_fit(sim$data) -head(sim$data) -plot(sim) -fit <- hystar_fit(sim$data) -devtools::load_all(".") -fit <- hystar_fit(sim$data) -fit$st_errors -devtools::load_all(".") -fit <- hystar_fit(sim$data) -fit$st_errors -devtools::load_all(".") -fit <- hystar_fit(sim$data) -fit$st_errors -confint(fit) -coef(fit) -usethis::use_version() -?do.call -check() -devtools::check() -devtools::load_all(".") -fit -fit$tar -?mapply -?identity -devtools::load_all(".") -y <- rnorm(100) -R <- c(rep(0, 50), rep(1, 50) -) -rv <- c(1, 1) -X <- matrix(c(rep(1, 99), y[1:99]), ncol = 2) -solve(t(X) %*% X) -X <- cbind(matrix(c(rep(1, 99), y[1:99]), ncol = 2) * (1-R), matrix(c(rep(1, 99), y[1:99]), ncol = 2) * R) -matrix(c(rep(1, 99), y[1:99]), ncol = 2) -matrix(c(rep(1, 99), y[1:99]), ncol = 2) * R -R <- matrix(c(rep(0, 50), rep(1, 50)), ncol = 2) -R -R <- matrix(c(rep(0, 50), rep(1, 50)), ncol = 2, nrow = 100) -R -X <- cbind(matrix(c(rep(1, 99), y[1:99]), ncol = 2) * (1-R), matrix(c(rep(1, 99), y[1:99]), ncol = 2) * R) -1_r -1-R -R <- matrix(c(rep(0, 49), rep(1, 50)), ncol = 2, nrow = 9) -R <- matrix(c(rep(0, 49), rep(1, 50)), ncol = 2, nrow = 99) -X <- cbind(matrix(c(rep(1, 99), y[1:99]), ncol = 2) * (1-R), matrix(c(rep(1, 99), y[1:99]), ncol = 2) * R) -X -solve(t(X) %*% X) -create_acov_mat(rnorm(100)) -create_acov_mat(1:5, 1:100) -acf(y) -y -plot(y) -acf(y) -?acf -acf(y, plot = FALSE) -compute_acov_vec(y, 10) -acf(y, plot = FALSE, type = "cov") -create_acov_mat(compute_acov_vec(y, 10)) -create_acov_mat(compute_acov_vec(y, 10), y) -create_acov_mat(compute_acov_vec(y, 3), y) -mean(y) -var(y) -var(y) -var(matrix(c(y, y), ncol = 2)) -var(matrix(c(y[2:100], y[1:99]), ncol = 2)) -devtools::load_all(".") -check() -devtools::check() -devtools::load_all(".") -devtools::load_all(".") -check() -library(devtools) -check() -?hystar_fit -z <- z_sim(n_t = 200, n_switches = 5, start_regime = 1) -sim <- hystar_sim(z = z, r = c(-.5, .5), d = 2, phi_R0 = c(0, .6), phi_R1 = 1) -fit <- hystar_fit(sim$data) -fit$st_errors -fit <- hystar_fit(sim$data, p0 = 2) -fit$st_errors -fit$eff -str(fit) -devtools::load_all(".") -compute_p_values(fit$coefficients, fit$st_errors) -devtools::load_all(".") -check() -devtools::load_all(".") -document() -check() -licence() -print -is.generic(print) -is.function(print) -class(print) -class(hystar.print) -class(print.hystar_fit) -devtools::load_all(".") -class(print.hystar_fit) -?UseMethod -plot -?hystar_fit -z <- z_sim(n_t = 200, n_switches = 5, start_regime = 1) -sim <- hystar_sim(z = z, r = c(-.5, .5), d = 2, phi_R0 = c(0, .6), phi_R1 = 1) -plot(sim) -fit <- hystar_fit(sim$data) -summary(fit) -plot(fit, regime_names = c("a", "b")) -devtools::load_all(".") -warnings() -plot(fit, xlab = "test") -devtools::load_all(".") -plot(fit, xlab = "test") -plot(fit, xlab = "test", regime_names = NULL) plot(fit_hystar, main = "Depression time series from a network model", xlab = "time", @@ -510,3 +370,143 @@ print_welcome_message() devtools::load_all(".") hystar_info() print_welcome_message() +library(devtools) +devtools::load_all(".") +print_welcome_message() +hystar_info() +devtools::load_all(".") +plot(hystar_sim(z_sim())) +?packageStartupMessage +devtools::load_all(".") +devtools::load_all(".") +hystar_info() +devtools::load_all(".") +?hystar_fit +?message +?Sys.sleep +devtools::load_all(".") +hystar_fit(hystar_sim(z_sim())$data, d = 1:10) +hystar_fit(hystar_sim(z_sim())$data, d = 1:5) +hystar_fit(hystar_sim(z_sim())$data, d = 1:5) +hystar_fit(hystar_sim(z_sim())$data, d = 1:10) +hystar_fit(hystar_sim(z_sim())$data, d = 1:10) +hystar_fit(hystar_sim(z_sim())$data, d = 1:6) +hystar_fit(hystar_sim(z_sim())$data, d = 1:7) +hystar_fit(hystar_sim(z_sim())$data, d = 1:8) +hystar_fit(hystar_sim(z_sim())$data, d = 1:9) +hystar_fit(hystar_sim(z_sim())$data, d = 1:10) +hystar_fit(hystar_sim(z_sim())$data, d = 1:11) +time_eff(1:100, 10, 2, 2) +create_x(hystar_sim(z_sim())$data$y, time_eff(1:100, 10, 2, 2), 1, 1) +hystar_fit(hystar_sim(z_sim())$data, p0 = 1:10) +hystar_fit(hystar_sim(z_sim())$data, p0 = 1:9) +hystar_fit(hystar_sim(z_sim())$data, p0 = 1:8) +hystar_fit(hystar_sim(z_sim())$data, p0 = 1:7) +hystar_fit(hystar_sim(z_sim())$data, p0 = 1:6) +hystar_fit(hystar_sim(z_sim())$data, p0 = 1:5) +hystar_fit(hystar_sim(z_sim(200))$data, p0 = 1:6) +hystar_fit(hystar_sim(z_sim(500))$data, d = 1:11) +hystar_fit(hystar_sim(z_sim(50))$data, p0 = 1:6) +?try +Inf +Inf >1000 +hystar_fit(hystar_sim(z_sim(50))$data, p0 = 1:) +hystar_fit(hystar_sim(z_sim(50))$data, p0 = 1:7) +devtools::load_all(".") +hystar_fit(hystar_sim(z_sim(50))$data, p0 = 1:7) +devtools::load_all(".") +hystar_fit(hystar_sim(z_sim(50))$data, p0 = 1:7) +(1)[ , 1, drop = TRUE] +devtools::load_all(".") +hystar_fit(hystar_sim(z_sim(50))$data, p0 = 1:7) +devtools::load_all(".") +hystar_fit(hystar_sim(z_sim(50))$data, p0 = 1:7) +?withCallingHandlers +devtools::load_all(".") +hystar_fit(hystar_sim(z_sim(50))$data, p0 = 1:7) +devtools::load_all(".") +hystar_fit(hystar_sim(z_sim(50))$data, p0 = 1:7) +devtools::load_all(".") +hystar_fit(hystar_sim(z_sim(150))$data, d = 1:6) +?message +devtools::load_all(".") +hystar_fit(hystar_sim(z_sim(150))$data, d = 1:6) +devtools::load_all(".") +hystar_fit(hystar_sim(z_sim(150))$data, d = 1:6) +devtools::load_all(".") +hystar_fit(hystar_sim(z_sim(150))$data, d = 1:6) +devtools::load_all(".") +hystar_fit(hystar_sim(z_sim(150))$data, d = 1:6) +devtools::load_all(".") +hystar_fit(hystar_sim(z_sim(150))$data, d = 1:6) +hystar_fit(hystar_sim(z_sim(150))$data, d = 1:6) +?txtProgressBar +txtProgressBar +devtools::load_all(".") +hystar_fit(hystar_sim(z_sim(150))$data, d = 1:6) +devtools::load_all(".") +devtools::load_all(".") +hystar_fit(hystar_sim(z_sim(150))$data, d = 1:6) +hystar_fit(hystar_sim(z_sim(150))$data, p0 = 1:5, p1 = 1:5) +hystar_fit(hystar_sim(z_sim(500))$data, p0 = 1:2, p1 = 1:2) +devtools::load_all(".") +hystar_fit(hystar_sim(z_sim(200))$data, p1 = 1:2) +stringi::stri_escape_unicode(🎉) +stringi::stri_escape_unicode("🎉") +document() +devtools::load_all(".") +devtools::load_all(".") +devtools::load_all(".") +devtools::load_all(".") +devtools::load_all(".") +devtools::load_all(".") +devtools::load_all(".") +devtools::load_all(".") +devtools::load_all(".") +devtools::load_all(".") +devtools::load_all(".") +devtools::load_all(".") +check( +) +devtools::load_all(".") +" +__ __ +/ /_ __ ______/ /_________ +/ _ / // (_ -/ _/ _ / __\ +/_//_/\_, /___)\__/\_,_/_/ +devtools::load_all(".") +devtools::load_all(".") +devtools::load_all(".") +hystar_fit(hystar_sim(z_sim())$data) +hystar_fit(hystar_sim(z_sim(300))$data, d = 3) +hystar_fit(hystar_sim(z_sim(300))$data, d = 1:3) +?hystar +hystar_fit(hystar_sim(z_sim(300))$data, d = 1:3)$ic +devtools::load_all(".") +hystar_fit(hystar_sim(z_sim(50))$data, d = 1:3)$ic +devtools::load_all(".") +hystar_fit(hystar_sim(z_sim(50))$data, d = 1:3)$ic +devtools::load_all(".") +hystar_fit(hystar_sim(z_sim(50))$data, d = 1:3)$ic +check() +devtools::load_all(".") +hystar_fit(hystar_sim(z_sim())$data) +devtools::load_all(".") +hystar_fit(hystar_sim(z_sim())$data) +devtools::load_all(".") +?invisible +test() +check() +document() +?hystar_info +?hystar +devtools::load_all(".") +check() +check() +hystar_fit(hystar_sim(z_sim()$data) +) +hystar_fit(hystar_sim(z_sim())$data) +devtools::load_all(".") +hystar_fit(hystar_sim(z_sim())$data) +hystar_fit(hystar_sim(z_sim())$data)$ic +usethis::use_version() diff --git a/NEWS.md b/NEWS.md index 501ac1c..9c0acbe 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,7 @@ # `hystar` 1.2.1.9000 (development version) +* Default argument for `p_select` in `hystar_fit()` now shows all four options. The new default is `p_select = c("bic", "aic", "aicc", "aiccp")`, and the first option `"bic"` is taken when the user doesn't provide the argument. + # `hystar` 1.2.1 * Added a welcome message with a logo and first directions for help. diff --git a/R/checks_fit.R b/R/checks_fit.R index d28b448..fd54c28 100644 --- a/R/checks_fit.R +++ b/R/checks_fit.R @@ -30,7 +30,7 @@ check_data <- function(data) { "You provided an object of class ", class(data)), call. = FALSE) - data <- data[, c(1, 2)] + data <- data[, c(1, 2), drop = FALSE] if (!is.numeric(data)) error_numeric(data) @@ -97,11 +97,12 @@ check_rz <- function(r, z) { } check_p_select <- function(p_select) { + if (length(p_select) > 1) p_select <- p_select[1] if (!is.character(p_select)) stop(paste0("`p_select` must be of type character."), call. = FALSE) p_select <- tolower(p_select) - choices <- c("aic", "aicc", "bic") + choices <- c("aic", "aicc", "bic", "aiccp") p_select <- tryCatch( error = function(cond) { stop(paste0("'p_select' must be one of these: ", diff --git a/R/hystar_fit.R b/R/hystar_fit.R index 725947a..24f5ad9 100644 --- a/R/hystar_fit.R +++ b/R/hystar_fit.R @@ -40,9 +40,10 @@ #' equal to `p0`. #' @param p_select The information criterion that should be minimized to select #' the orders \eqn{p_0} and \eqn{p_1}. Choices: +#' * `"bic"` (default, Bayesian Information Criterion) #' * `"aic"` (Akaike Information Criterion) #' * `"aicc"` (Corrected Akaike Information Criterion) -#' * `"bic"` (default, Bayesian Information Criterion) +#' * `"aiccp"` (Change-point Akaike Information Criterion) #' @param thin `TRUE` (default) or `FALSE`. Only relevant when `r` is a vector. #' * If `TRUE` (default), the search space for the thresholds are the #' \eqn{100a\%, 100(a+0.01)\%, \dots, 100b\%} percentiles of `z`. @@ -101,8 +102,13 @@ #' * `nobs()` #' #' @export -hystar_fit <- function(data, r = c(.1, .9), d = 0L, p0 = 1L, p1 = 1L, p_select = "bic", - thin = FALSE, tar = FALSE) { +hystar_fit <- function(data, + r = c(.1, .9), + d = 0L, + p0 = 1L, p1 = 1L, + p_select = c("bic", "aic", "aicc", "aiccp"), + thin = FALSE, + tar = FALSE) { check_data(data) if (is.vector(data)) { y <- z <- data @@ -110,6 +116,7 @@ hystar_fit <- function(data, r = c(.1, .9), d = 0L, p0 = 1L, p1 = 1L, p_select = y <- data[, 1] z <- data[, 2] } + p_select <- check_hystar_fit_input(z, d, p0, p1, p_select, r, thin, tar) eff <- time_eff(y, max(d), max(p0), max(p1)) x <- create_x(y, eff, max(p0), max(p1)) diff --git a/R/hystar_sim.R b/R/hystar_sim.R index aa79d75..99e3253 100644 --- a/R/hystar_sim.R +++ b/R/hystar_sim.R @@ -108,7 +108,8 @@ hystar_sim <- function(z, d = 0, phi_R0 = c(0, .5), phi_R1 = c(2, .5), - resvar = c(1, 1), start_regime = NULL) { + resvar = c(1, 1), + start_regime = NULL) { temp <- check_hystar_sim_input(z, r, d, phi_R0, phi_R1, resvar, start_regime) z <- temp$z diff --git a/R/plot.R b/R/plot.R index 4550efc..85a9ccf 100644 --- a/R/plot.R +++ b/R/plot.R @@ -1,22 +1,30 @@ #' @export plot.hystar_fit <- function(x, y = NULL, main = "Fitted HysTAR model", ...) { plot_hystar(hystar_object = x, main = main, ...) - return(NULL) + invisible(NULL) } #' @export plot.hystar_sim <- function(x, y = NULL, main = "Simulated HysTAR model", ...) { plot_hystar(hystar_object = x, main = main, ...) - return(NULL) + invisible(NULL) } #' @importFrom graphics par plot_hystar <- function(hystar_object, main, - regimes_name_color = c("Regime 0" = "white", - "Regime 1" = "grey85", - "Unknown" = "#AA001137"), + name_regime0 = "Regime 0", + name_regime1 = "Regime 1", + color_regime0 = "white", + color_regime1 = "grey85", ...) { + + regimes_name_color <- c(name_regime0 = name_regime0, + color_regime0 = color_regime0, + name_regime1 = name_regime1, + color_regime1 = color_regime1, + name_regime_unknown = "Unknown", + color_regime_unknown = "#AA001137") old <- par(mfrow = c(2, 1)) on.exit(par(old), add = TRUE) plot_z(hystar_object, main, regimes_name_color, ...) @@ -82,13 +90,13 @@ plot_y <- function(hystar_object, #' @importFrom graphics legend plot_legend <- function(hystar_object, regimes_name_color) { n_unknown <- sum(is.na(hystar_object$data$R)) - if (n_unknown == 0) regimes_name_color <- regimes_name_color[-3] + if (n_unknown == 0) regimes_name_color <- regimes_name_color[-c(5, 6)] starting_regime <- hystar_object$data$R[!is.na(hystar_object$data$R)][1] legend_position <- if (starting_regime == 1) "bottomleft" else "topleft" legend( x = legend_position, - legend = names(regimes_name_color), - fill = regimes_name_color, + legend = regimes_name_color[c("name_regime0", "name_regime1")], + fill = regimes_name_color[c("color_regime0", "color_regime1")], bg = "grey95", cex = .6 ) @@ -114,7 +122,7 @@ plot_background <- function(hystar_object, regimes_name_color, type) { ybottom = plot_region_borders["ybottom"], xright = plot_region_borders["xright"], ytop = plot_region_borders["ytop"], - col = regimes_name_color[[1]], + col = regimes_name_color["color_regime0"], border = NA ) # Background coloring for regime 1 @@ -123,7 +131,7 @@ plot_background <- function(hystar_object, regimes_name_color, type) { ybottom = rep(plot_region_borders["ybottom"], times = nrow(switch_points_matrix)), xright = switch_points_matrix[, 2], ytop = rep(plot_region_borders["ytop"], times = nrow(switch_points_matrix)), - col = regimes_name_color[[2]], + col = regimes_name_color["color_regime1"], border = NA ) # Background for the unknown regimes: @@ -136,7 +144,7 @@ plot_background <- function(hystar_object, regimes_name_color, type) { ybottom = plot_region_borders["ybottom"], xright = k + 1, ytop = plot_region_borders["ytop"], - col = regimes_name_color[[3]], + col = regimes_name_color["color_regime_unknown"], border = NA ) } diff --git a/R/zzz.R b/R/zzz.R index 79ad11a..4a86783 100644 --- a/R/zzz.R +++ b/R/zzz.R @@ -1,5 +1,5 @@ .onAttach <- function(libname, pkgname) { -hystar_string <- " +packageStartupMessage(" __ __ / /_ __ ______/ /_________ / _ / // (_ -/ _/ _ / __\\ @@ -9,9 +9,7 @@ hystar_string <- " Estimation and simulation of the HysTAR Model. For function help, run `?hystar_fit`, `?hystar_sim` or `?z_sim`. For more information, run `hystar_info()` (opens a URL in your browser). -" - -packageStartupMessage(hystar_string) +") } #' Get more information about the hystar package @@ -23,7 +21,7 @@ packageStartupMessage(hystar_string) #' @importFrom utils browseURL hystar_info <- function() { browseURL("https://daandejongen.github.io/hystar/") - return(NULL) + invisible(NULL) } diff --git a/_pkgdown.yml b/_pkgdown.yml index 6c62429..8b66225 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -8,3 +8,4 @@ reference: - hystar_fit - hystar_sim - z_sim + - hystar_info diff --git a/man/hystar_fit.Rd b/man/hystar_fit.Rd index 8a852b6..1b52100 100644 --- a/man/hystar_fit.Rd +++ b/man/hystar_fit.Rd @@ -10,7 +10,7 @@ hystar_fit( d = 0L, p0 = 1L, p1 = 1L, - p_select = "bic", + p_select = c("bic", "aic", "aicc", "aiccp"), thin = FALSE, tar = FALSE ) @@ -48,9 +48,10 @@ equal to \code{p0}.} \item{p_select}{The information criterion that should be minimized to select the orders \eqn{p_0} and \eqn{p_1}. Choices: \itemize{ +\item \code{"bic"} (default, Bayesian Information Criterion) \item \code{"aic"} (Akaike Information Criterion) \item \code{"aicc"} (Corrected Akaike Information Criterion) -\item \code{"bic"} (default, Bayesian Information Criterion) +\item \code{"aiccp"} (Change-point Akaike Information Criterion) }} \item{thin}{\code{TRUE} (default) or \code{FALSE}. Only relevant when \code{r} is a vector. diff --git a/tests/testthat/test-checks_fit.R b/tests/testthat/test-checks_fit.R index 0d14a7f..3b9c519 100644 --- a/tests/testthat/test-checks_fit.R +++ b/tests/testthat/test-checks_fit.R @@ -25,6 +25,10 @@ test_that("thin must be TRUE or FALSE", { expect_error(check_thin("a"), "TRUE or FALSE") }) +test_that("tar must be TRUE or FALSE", { + expect_error(check_tar(2), "TRUE or FALSE") +}) + test_that("p_select must be a valid choice", { expect_error(hystar_fit(data.frame(y = 1:10, z = 1:10), p_select = 1), "character")