Improved treatment of doubling / halving time when growth rate estimate spans 0

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)
 ============================

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

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)
+})

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)