Merge pull request #12 from reconhub/issue_11

fix issue 11, still needs handling numer approx pb

.travis.yml

@@ -22,8 +22,8 @@ matrix:
     - os: osx
       osx_image: xcode8.3
 
-r_github_packages:
-  - jimhester/covr
+# r_github_packages:
+#   - jimhester/covr
 
 warnings_are_errors: true
 

DESCRIPTION

@@ -1,10 +1,11 @@
 Package: earlyR
 Title: Estimation of Transmissibility in the Early Stages of a Disease Outbreak
-Version: 0.0.1
+Version: 0.0.3
 Authors@R: c(person("Thibaut", "Jombart", email = "[email protected]", role = c("aut", "cre")),
     person("Anne", "Cori", email = "[email protected]", role = c("aut")),
     person("Pierre", "Nouvellet", email = "[email protected]", role = c("aut")),
-    person("Janetta", "Skarp", email = "[email protected]", role = c("aut")))
+    person("Janetta", "Skarp", email = "[email protected]", role = c("aut")),
+    person("Zhian N.", "Kamvar", email = "[email protected]", role = c("ctb")))
 Description: Implements a simple, likelihood-based estimation of the reproduction number (R0) using a branching process with a Poisson likelihood. This model requires knowledge of the serial interval distribution, and dates of symptom onsets. Infectiousness is determined by weighting R0 by the probability mass function of the serial interval on the corresponding day. It is a simplified version of the model introduced by Cori et al. (2013) <doi:10.1093/aje/kwt133>.
 Depends: R (>= 3.3.0)
 License: MIT + file LICENSE

NEWS.md

@@ -1,3 +1,9 @@
+# earlyR 0.0.3
+
+## Bug fixes
+
+- fix issue #11 where the likelihood surface was incorrectly calculated.
+  This unfortunately reverts the issue from #4
 
 # earlyR 0.0.2
 

R/internals.R

@@ -53,16 +53,12 @@
 
 
 
-## Sanitize very low log-likelihood values and put them back on their original
-## scale, summing to 1. We need to consider the max of log-like could be 0.
+## This function converts log-likelihood values back to their original
+## scales. It also aims to sanitize very low log-likelihood values which will
+## cause numerical approximation problems when converted to the original scale.
 
 loglike_to_density <- function(x) {
-  out <- x
-  x_max <- max(x, na.rm = TRUE)
-  if (x_max != 0) {
-    out <- out / abs(x_max)
-  }
-  out <- exp(out)
-  out <- out / sum(out)
+  out <- exp(x)
+  out <- out / sum(out, na.rm = TRUE)
   out
 }

tests/testthat/test_get_R.R

@@ -24,7 +24,8 @@ test_that("Estimation robust to heading zeros", {
   dat <- as.Date("2019-01-01") + sample(1:10, 50, replace = TRUE)
   i <- incidence(dat)
   i_early <- incidence(dat,
-                       first_date = as.Date("2018-12-15"))
+                       first_date = as.Date("2018-12-15"),
+                       standard = TRUE)
   ## example with a function for SI
   si <- distcrete("gamma", interval = 1L,
                   shape = 1.5,
@@ -132,3 +133,46 @@ test_that("Errors are thrown when they should", {
 
 })
 
+
+test_that("issue 11 isn't compromised", {
+
+  params <- list(shape = 1.88822148063956, scale = 4.5613778865727)
+  si     <- distcrete::distcrete("gamma", 
+                                 shape = params$shape, 
+                                 scale = params$scale, 
+                                 interval = 1L, w = 0L)
+  i_df <- data.frame(
+         dates = c("2014-04-07", "2014-04-08", "2014-04-09", "2014-04-10",
+                   "2014-04-11", "2014-04-12", "2014-04-13", "2014-04-14",
+                   "2014-04-15", "2014-04-16", "2014-04-17", "2014-04-18", "2014-04-19",
+                   "2014-04-20", "2014-04-21", "2014-04-22", "2014-04-23",
+                   "2014-04-24", "2014-04-25", "2014-04-26", "2014-04-27", "2014-04-28",
+                   "2014-04-29", "2014-04-30", "2014-05-01", "2014-05-02",
+                   "2014-05-03", "2014-05-04", "2014-05-05", "2014-05-06", "2014-05-07",
+                   "2014-05-08", "2014-05-09", "2014-05-10", "2014-05-11",
+                   "2014-05-12", "2014-05-13", "2014-05-14", "2014-05-15", "2014-05-16",
+                   "2014-05-17", "2014-05-18", "2014-05-19", "2014-05-20", "2014-05-21",
+                   "2014-05-22", "2014-05-23", "2014-05-24", "2014-05-25",
+                   "2014-05-26", "2014-05-27", "2014-05-28", "2014-05-29", "2014-05-30",
+                   "2014-05-31", "2014-06-01", "2014-06-02", "2014-06-03",
+                   "2014-06-04", "2014-06-05", "2014-06-06", "2014-06-07", "2014-06-08",
+                   "2014-06-09", "2014-06-10", "2014-06-11", "2014-06-12",
+                   "2014-06-13", "2014-06-14", "2014-06-15", "2014-06-16", "2014-06-17"),
+        counts = c(1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 2L,
+                   0L, 0L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 2L, 0L, 1L, 1L, 1L, 3L,
+                   2L, 2L, 2L, 1L, 1L, 3L, 3L, 3L, 0L, 2L, 3L, 3L, 3L, 1L, 3L, 1L,
+                   2L, 2L, 3L, 4L, 10L, 1L, 1L, 1L, 2L, 0L, 5L, 3L, 0L, 4L, 5L, 5L,
+                   1L, 5L, 4L, 3L, 1L, 1L, 2L, 5L, 5L, 4L)
+  )
+
+  i_trunc <- as.incidence(i_df[-1], dates = as.Date(i_df$dates))
+  i_trunc
+  the_R <- get_R(i_trunc, si = si)
+  # The estimate should be 1.24
+  expect_identical(round(the_R$R_ml, 2), 1.24)
+
+  # The likelihood surface should pretty much be contained between 0.24 and 2.24
+  expect_equal(sum(with(the_R, R_like[R_grid <= 2.24 & R_grid >= 0.24])), 1)
+
+})
+