diff --git a/NEWS.md b/NEWS.md index 5b274c1..5474c64 100644 --- a/NEWS.md +++ b/NEWS.md @@ -11,12 +11,17 @@ incidence 1.5.3 (2018-12-07) (See https://github.com/reconhub/incidence/issues/88) * `fit()` now returns correct coefficients when dates is POSIXt by converting to Date. (See https://github.com/reconhub/incidence/issues/91) +* `fit()` now returns Inf for the upper bound of the doubling / halving time if the + confidence interval on the growth rate spans 0. + (See https://github.com/reconhub/incidence/issues/28) ### MISC * `demo("incidence-demo" package = "incidence")` has been updated to show use of custom colors. -* `incidence()` no longer accepts characters as input for dates, first_date, or last_date argeuments +* `incidence()` no longer accepts characters as input for dates, first_date, or last_date arguments +* `fit()` now returns an error if growth rate estimates for groups have different signs. + (See https://github.com/reconhub/incidence/issues/28) incidence 1.5.2 (2018-11-30) ============================ diff --git a/R/extract_info.R b/R/extract_info.R index d4e73e1..9ebdb9c 100644 --- a/R/extract_info.R +++ b/R/extract_info.R @@ -14,6 +14,11 @@ extract_info <- function(reg, origin, level){ r <- stats::coef(reg)[to.keep] use.groups <- length(r) > 1 if (use.groups) { + if (all(r >= 0) | all(r <= 0)) { + # continue + } else { + stop("Growth rates of groups have different signs; fit groups separately.") + } names(r) <- reg$xlevels[[1]] # names = levels if groups } else { names(r) <- NULL # no names otherwise @@ -47,10 +52,12 @@ extract_info <- function(reg, origin, level){ o.names <- colnames(info$doubling.conf) info$doubling.conf <- info$doubling.conf[, rev(seq_along(o.names)), drop = FALSE] + info$doubling.conf[info$doubling.conf < 0] <- Inf colnames(info$doubling.conf) <- o.names } else { info$halving <- log(0.5) / r info$halving.conf <- log(0.5) / r.conf + info$halving.conf[info$halving.conf < 0] <- Inf } ## We need to store the date corresponding to 'day 0', as this will be used diff --git a/tests/testthat/rds/fit.i.rds b/tests/testthat/rds/fit.i.rds index d080762..ec64fca 100644 Binary files a/tests/testthat/rds/fit.i.rds and b/tests/testthat/rds/fit.i.rds differ diff --git a/tests/testthat/rds/fit.i.sex.rds b/tests/testthat/rds/fit.i.sex.rds index 3b6bc70..d38e334 100644 Binary files a/tests/testthat/rds/fit.i.sex.rds and b/tests/testthat/rds/fit.i.sex.rds differ diff --git a/tests/testthat/rds/print.fit.i.rds b/tests/testthat/rds/print.fit.i.rds index f09743a..21e3a08 100644 Binary files a/tests/testthat/rds/print.fit.i.rds and b/tests/testthat/rds/print.fit.i.rds differ diff --git a/tests/testthat/rds/print.fit.sex.rds b/tests/testthat/rds/print.fit.sex.rds index c72313e..3533b04 100644 Binary files a/tests/testthat/rds/print.fit.sex.rds and b/tests/testthat/rds/print.fit.sex.rds differ diff --git a/tests/testthat/test-fit.R b/tests/testthat/test-fit.R index 824f1de..579e765 100644 --- a/tests/testthat/test-fit.R +++ b/tests/testthat/test-fit.R @@ -5,7 +5,7 @@ test_that("fit", { set.seed(1) dat <- c(sample(1:50, 200, replace = TRUE, prob = 1 + exp(1:50 * 0.1)), - sample(51:100, 200, replace = TRUE, prob = rev(1 + exp(1:50 * 0.1)))) + sample(51:100, 200, replace = TRUE, prob = 1 + exp(1:50 * 0.11))) sex <- sample(c("female", "male"), 400, replace = TRUE) i <- incidence(dat, 5L) @@ -14,10 +14,10 @@ test_that("fit", { expect_warning(fit.i.sex <- fit(i.sex), "3 dates with incidence of 0 ignored for fitting") - expect_equal_to_reference(fit.i, file = "rds/fit.i.rds") - expect_equal_to_reference(fit.i.sex, file = "rds/fit.i.sex.rds") - expect_equal_to_reference(capture.output(fit.i), file = "rds/print.fit.i.rds") - expect_equal_to_reference(capture.output(fit.i.sex), file = "rds/print.fit.sex.rds") + expect_known_value(fit.i, file = "rds/fit.i.rds", update = FALSE) + expect_known_value(fit.i.sex, file = "rds/fit.i.sex.rds", update = FALSE) + expect_known_value(capture.output(fit.i), file = "rds/print.fit.i.rds", update = FALSE) + expect_known_value(capture.output(fit.i.sex), file = "rds/print.fit.sex.rds", update = FALSE) ## errors x <- incidence(c(1, 0, 0, 0), interval = 7) @@ -36,10 +36,10 @@ test_that("fit_optim_split", { i.sex <- incidence(dat, 5L, groups = sex) i.fit <- fit_optim_split(i, plot = FALSE) i.fit.sex <- fit_optim_split(i.sex, plot = FALSE) - expect_equal_to_reference(i.fit, - file = "rds/o.fit.i.rds") - expect_equal_to_reference(i.fit.sex, - file = "rds/o.fit.i.sex.rds") + expect_known_value(i.fit, + file = "rds/o.fit.i.rds", update = FALSE) + expect_known_value(i.fit.sex, + file = "rds/o.fit.i.sex.rds", update = FALSE) expect_is(i.fit$df, "data.frame") expect_is(i.fit$fit, "incidence_fit_list") @@ -73,3 +73,33 @@ test_that("fitting results are the same for incidence fits on Dates and POSIXct" expect_equal(fit(iP), fit(iD)) }) + +test_that("doubling / halving time makes sense when CI of r crosses 0", { + + set.seed(20181213) + days <- 1:14 + # estimate of r is positive + dat_cases_1 <- rnbinom(14, .5, .1) + dat_dates_1 <- rep(as.Date(Sys.Date() + days), dat_cases_1) + + i1 <- incidence(dat_dates_1) + f1 <- suppressWarnings(fit(i1)) + + expect_true(any(is.infinite(f1$info$doubling.conf))) + + # estimate of r is negative + dat_cases_2 <- rnbinom(14, .5, .1) + dat_dates_2 <- rep(as.Date(Sys.Date() + days), dat_cases_2) + + i2 <- incidence(dat_dates_2) + f2 <- suppressWarnings(fit(i2)) + + expect_true(any(is.infinite(f2$info$halving.conf))) + + # groups have different signs for r + grp <- rep(c("grp1", "grp2"), c(length(dat_dates_1), length(dat_dates_2))) + i.grp <- incidence(c(dat_dates_1, dat_dates_2), groups = grp) + + msg <- "Growth rates of groups have different signs; fit groups separately." + expect_error(suppressWarnings(fit(i.grp)), msg) +}) diff --git a/tests/testthat/test-plot.R b/tests/testthat/test-plot.R index 2451eda..ddb0afc 100644 --- a/tests/testthat/test-plot.R +++ b/tests/testthat/test-plot.R @@ -3,7 +3,7 @@ context("Test plotting") test_that("plot for incidence object", { skip_on_cran() - set.seed(1) + set.seed(12) dat <- sample(1:50, 200, replace = TRUE, prob = 1 + exp(1:50 * 0.1)) dat2 <- as.Date("2016-01-02") + dat dat3 <- as.POSIXct(dat2)