Weibull delay support + distspec alignment (0.2.0)

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: cfrnow
 Title: Real-Time Case Fatality Ratio via a Mixture-Cure Survival Model
-Version: 0.1.0.9000
+Version: 0.2.0
 Authors@R:
     person("Sebastian", "Funk", , "sebastian.funk@lshtm.ac.uk",
            role = c("aut", "cre"))
@@ -18,7 +18,7 @@ BugReports: https://github.com/sbfnk/cfrnow/issues
 Encoding: UTF-8
 Roxygen: list(markdown = TRUE)
 Depends:
-    distspec,
+    distspec (>= 0.1.0),
     R (>= 4.1.0)
 Imports:
     brms,

NEWS.md

@@ -1,5 +1,11 @@
-# cfrnow (development version)
+# cfrnow 0.2.0
 
+* `fit_cfr()` and `simulate_linelist()` support a `Weibull()` onset-to-death (and
+  recovery) delay, alongside `LogNormal()` and `Gamma()`.
+* Delay parameterisation now uses distspec's exported `natural_params()` in place
+  of an internal helper, tracking the distspec API.
+* Added a "Stratified and partially-pooled CFR" vignette covering no-, complete-
+  and partial-pooling CFR fits and per-group `summary()` output.
 * `pp_check_cfr()` runs a posterior-predictive check on a fit: it draws replicate
   line-list outcomes from the posterior, replays the real-time truncation, and
   compares the observed death counts (plus recoveries in a two-outcome fit) and

R/cure_model.R

@@ -16,7 +16,7 @@
 #'
 #' The delay distribution's location is `mu` (as `epidist` expects); the cure
 #' probability `cfr` is an additional dpar with a logit link. Supported delay
-#' families are `lognormal()` and `Gamma()`.
+#' families are `lognormal()`, `Gamma()` and `Weibull()`.
 #'
 #' @name cfrnow-cure-model
 #' @keywords internal
@@ -81,8 +81,8 @@ assert_epidist.epidist_cure_model <- function(data, ...) {
 }
 
 .assert_delay_family <- function(family, what = "delay") {
-  if (!family$family %in% c("lognormal", "gamma")) {
-    stop("cfrnow supports lognormal() and Gamma() ", what, "s only.",
+  if (!family$family %in% c("lognormal", "gamma", "weibull")) {
+    stop("cfrnow supports lognormal(), Gamma() and Weibull() ", what, "s only.",
       call. = FALSE
     )
   }
@@ -130,19 +130,18 @@ epidist_model_prior.epidist_cure_model <- function(data, formula, ...) NULL
 
 # The native Stan parameterisation brms uses for a family's mu-form, read out of
 # a dummy model (as epidist does), so primarycensored gets the right arguments.
+# The mu-form arguments may themselves contain parentheses (e.g. weibull's
+# `mu ./ tgamma(1 + 1 ./ shape)`), so the args run to the statement's closing
+# `);` rather than the first `)`.
 .family_stan_param <- function(family_name) {
   code <- brms::make_stancode(y ~ 1,
     data = data.frame(y = c(1, 2)),
     family = family_name
   )
-  pat <- sprintf("target \\+= %s_(lpdf|lpmf)\\(Y \\| ([^)]+)\\)", family_name)
-  hits <- grep("mu", regmatches(code, gregexpr(pat, code))[[1]],
-    fixed = TRUE, value = TRUE
-  )
-  sub(")", "", sub(
-    sprintf("target \\+= %s_(lpdf|lpmf)\\(Y \\| ", family_name),
-    "", hits[1]
-  ), fixed = TRUE)
+  pat <- sprintf("%s_l[a-z]+\\(Y \\| (.+?)\\);", family_name)
+  stmts <- regmatches(code, gregexpr(pat, code, perl = TRUE))[[1]]
+  stmt <- grep("mu", stmts, fixed = TRUE, value = TRUE)[1]
+  sub(pat, "\\1", stmt, perl = TRUE)
 }
 
 # Fill a `<<hole>>` template from inst/stan with a named list of replacements.

R/fit_cfr.R

@@ -8,8 +8,9 @@
 #'
 #' The delay's native parameters may each be a fixed number (held fixed; fixing
 #' the whole delay gives the Ghani/Nishiura estimator) or a `Normal()` prior
-#' (co-estimated). The family (`LogNormal()` or `Gamma()`) sets the delay
-#' distribution. `cfr_prior` is a `Beta()`; it matters because the CFR is weakly
+#' (co-estimated). The family (`LogNormal()`, `Gamma()` or `Weibull()`) sets the
+#' delay distribution. `cfr_prior` is a `Beta()`; it matters because the CFR is
+#' weakly
 #' identified early on (`Beta(1, 1)` is uniform, `Beta(1, 9)` favours a low CFR,
 #' `Beta(6.6, 13.4)` suits a high-fatality pathogen).
 #'
@@ -22,8 +23,8 @@
 #' @param data Output of [prepare_cfr_data()], or an `epidist_cure_model` /
 #'   data frame with `y`, `outcome`, `pwindow`, `swindow`.
 #' @param delay Onset-to-death delay as a \pkg{distspec} distribution
-#'   ([distspec::LogNormal()] or [distspec::Gamma()]) whose native parameters
-#'   are fixed numbers or `Normal()` priors.
+#'   ([distspec::LogNormal()], [distspec::Gamma()] or [distspec::Weibull()])
+#'   whose native parameters are fixed numbers or `Normal()` priors.
 #' @param cfr_prior CFR prior as a [distspec::Beta()]. Defaults to `Beta(1, 1)`.
 #' @param recovery_delay Optional onset-to-recovery delay (same form as `delay`)
 #'   for the two-outcome fit; may use a different family from `delay`.

R/simulate_linelist.R

@@ -20,7 +20,8 @@
 #' @param n Number of cases.
 #' @param cfr True case fatality ratio.
 #' @param delay Onset-to-death delay: a distspec distribution
-#'   ([distspec::LogNormal()] or [distspec::Gamma()]) with fixed parameters.
+#'   ([distspec::LogNormal()], [distspec::Gamma()] or [distspec::Weibull()])
+#'   with fixed parameters.
 #' @param recovery Optional onset-to-recovery delay (same form as `delay`); when
 #'   given, non-fatal cases get a `recovery_date`.
 #' @param onset_start First possible onset date.
@@ -74,15 +75,26 @@ simulate_linelist <- function(n = 200, cfr = 0.5, delay, recovery = NULL,
 #' @noRd
 sample_delay <- function(n, delay) {
   fam <- get_distribution(delay)
-  pars <- get_parameters(delay)[delay_native_order(fam)]
+  pars <- get_parameters(delay)[natural_params(delay)]
   if (!all(vapply(pars, is.numeric, logical(1)))) {
     stop("simulate_linelist() needs a delay with fixed parameters (numbers), ",
       "not priors.",
       call. = FALSE
     )
   }
-  switch(fam,
+  out <- switch(fam,
     lognormal = stats::rlnorm(n, pars[["meanlog"]], pars[["sdlog"]]),
-    gamma = stats::rgamma(n, shape = pars[["shape"]], rate = pars[["rate"]])
+    gamma = stats::rgamma(n, shape = pars[["shape"]], rate = pars[["rate"]]),
+    weibull = stats::rweibull(
+      n,
+      shape = pars[["shape"]], scale = pars[["scale"]]
+    )
   )
+  if (is.null(out)) {
+    stop("simulate_linelist() supports LogNormal(), Gamma() and Weibull() ",
+      "delays only.",
+      call. = FALSE
+    )
+  }
+  out
 }

R/summarise_cfr.R

@@ -11,14 +11,19 @@ naive_cfr <- function(n_deaths, n_cases) {
 #'
 #' `loc_var` is the location intercept and `scale_var` the second parameter's
 #' intercept; the second parameter is log-linked. Handles lognormal (loc is
-#' meanlog) and gamma (loc is log-mean, second parameter is shape).
+#' meanlog), gamma and weibull (loc is log-mean, second parameter is the shape).
 #' @noRd
 .delay_moments <- function(dr, loc_var, scale_var, family) {
   loc <- dr[[loc_var]]
   sc <- exp(dr[[scale_var]])
   if (family == "lognormal") {
     dmean <- exp(loc + sc^2 / 2)
     list(mean = dmean, sd = sqrt(exp(sc^2) - 1) * dmean)
+  } else if (family == "weibull") {
+    # brms weibull: mu is the mean (log link), sc is the shape k
+    dmean <- exp(loc)
+    cv <- sqrt(gamma(1 + 2 / sc) / gamma(1 + 1 / sc)^2 - 1)
+    list(mean = dmean, sd = dmean * cv)
   } else {
     dmean <- exp(loc)
     list(mean = dmean, sd = dmean / sqrt(sc))

R/translate.R

@@ -66,6 +66,37 @@
   )
 }
 
+# distspec Weibull() uses (shape, scale); brms weibull uses (mu = mean, shape),
+# both log-linked. Build the mu prior from mean = scale * gamma(1 + 1 / shape),
+# propagating shape/scale uncertainty to log(mu) by the delta method.
+.weibull_mu_prior <- function(shape_spec, scale_spec, dpar) {
+  mu_of <- function(k, lambda) lambda * gamma(1 + 1 / k)
+  if (!is.null(shape_spec$fixed) && !is.null(scale_spec$fixed)) {
+    return(.prior_row(
+      list(fixed = mu_of(shape_spec$fixed, scale_spec$fixed)), dpar, TRUE
+    ))
+  }
+  kc <- shape_spec$fixed %||% shape_spec$mean
+  lc <- scale_spec$fixed %||% scale_spec$mean
+  # d log(mu)/d shape uses digamma(1 + 1/k); sign drops on squaring
+  d_shape <- digamma(1 + 1 / kc) / kc^2
+  vlog <- sqrt(
+    ((scale_spec$sd %||% 0) / lc)^2 + (d_shape * (shape_spec$sd %||% 0))^2
+  )
+  prior_args <- if (nzchar(dpar)) {
+    list(class = "Intercept", dpar = dpar)
+  } else {
+    list(class = "Intercept")
+  }
+  do.call(
+    brms::set_prior,
+    c(
+      list(sprintf("normal(%.6f, %.6f)", log(lc) + lgamma(1 + 1 / kc), vlog)),
+      prior_args
+    )
+  )
+}
+
 # A distspec delay -> list(family, prior). `main = TRUE` targets the death delay
 # (mu is the main dpar); `main = FALSE` the recovery delay (r-prefixed dpars).
 .delay_family_prior <- function(delay, main = TRUE) {
@@ -78,25 +109,39 @@
   fam <- get_distribution(delay)
   pars <- get_parameters(delay)
   loc_dpar <- if (main) "" else "rmu"
-  if (fam == "lognormal") {
-    scale_dpar <- if (main) "sigma" else "rsigma"
-    prior <- c(
-      .prior_row(.param_spec(pars$meanlog), loc_dpar, FALSE),
-      .prior_row(.param_spec(pars$sdlog), scale_dpar, TRUE)
-    )
-    list(family = brms::lognormal(), prior = prior)
-  } else if (fam == "gamma") {
-    scale_dpar <- if (main) "shape" else "rshape"
-    sh <- .param_spec(pars$shape)
-    rate_spec <- .param_spec(pars$rate)
-    prior <- c(
-      .gamma_mu_prior(sh, rate_spec, loc_dpar),
-      .prior_row(sh, scale_dpar, TRUE)
+  scale_dpar <- if (main) "sigma" else "rsigma"
+  shape_dpar <- if (main) "shape" else "rshape"
+  out <- switch(fam,
+    lognormal = {
+      prior <- c(
+        .prior_row(.param_spec(pars$meanlog), loc_dpar, FALSE),
+        .prior_row(.param_spec(pars$sdlog), scale_dpar, TRUE)
+      )
+      list(family = brms::lognormal(), prior = prior)
+    },
+    gamma = {
+      sh <- .param_spec(pars$shape)
+      prior <- c(
+        .gamma_mu_prior(sh, .param_spec(pars$rate), loc_dpar),
+        .prior_row(sh, shape_dpar, TRUE)
+      )
+      list(family = stats::Gamma(link = "log"), prior = prior)
+    },
+    weibull = {
+      sh <- .param_spec(pars$shape)
+      prior <- c(
+        .weibull_mu_prior(sh, .param_spec(pars$scale), loc_dpar),
+        .prior_row(sh, shape_dpar, TRUE)
+      )
+      list(family = brms::weibull(), prior = prior)
+    }
+  )
+  if (is.null(out)) {
+    stop("cfrnow supports LogNormal(), Gamma() and Weibull() delays only.",
+      call. = FALSE
     )
-    list(family = stats::Gamma(link = "log"), prior = prior)
-  } else {
-    stop("cfrnow supports LogNormal() and Gamma() delays only.", call. = FALSE)
   }
+  out
 }
 
 # A Normal(m, s) on logit(cfr) whose induced mean/sd on the [0, 1] CFR scale

tests/testthat/test-families.R

@@ -2,7 +2,8 @@ test_that("simulated delays match the requested mean and sd", {
   set.seed(1)
   specs <- list(
     LogNormal(mean = 12.75, sd = 7),
-    Gamma(mean = 12.75, sd = 7)
+    Gamma(mean = 12.75, sd = 7),
+    Weibull(mean = 12.75, sd = 7)
   )
   for (delay in specs) {
     ll <- simulate_linelist(n = 20000, cfr = 1, delay = delay)
@@ -28,7 +29,7 @@ test_that("a recovery delay adds recoveries for non-fatal cases only", {
 })
 
 test_that("unsupported families and prior-parameter delays are rejected", {
-  expect_error(delay_native_order("weibull"), "unsupported")
+  expect_error(sample_delay(5, Exponential(rate = 1)), "supports") # exp not supported
   expect_error(sample_delay(5, LogNormal(
     meanlog = Normal(2, 0.1), sdlog = Normal(0.5, 0.1)
   )), "fixed")

tests/testthat/test-fit_cfr.R

@@ -60,6 +60,45 @@ test_that("a gamma delay recovers the CFR and delay moments", {
   expect_equal(s[s$quantity == "delay_mean", "q50"], 8, tolerance = 0.8)
 })
 
+test_that("a weibull delay recovers the CFR and delay moments", {
+  skip_if_no_cmdstan()
+  set.seed(3)
+  ll <- simulate_linelist(
+    n = 2000, cfr = 0.4, onset_days = 40,
+    delay = Weibull(shape = 1.5, scale = 13)
+  )
+  d <- prepare_cfr_data(ll, obs_time = max(ll$onset_date) - 2)
+  s <- summary(fit_quick(d,
+    delay = Weibull(shape = Normal(1.5, 0.4), scale = Normal(13, 3)),
+    cfr_prior = Beta(1, 1)
+  ))
+  true_mean <- 13 * gamma(1 + 1 / 1.5) # Weibull mean = scale * Gamma(1 + 1/shape)
+  expect_lt(s[s$quantity == "cfr", "q2.5"], 0.4)
+  expect_gt(s[s$quantity == "cfr", "q97.5"], 0.4)
+  expect_equal(s[s$quantity == "delay_mean", "q50"], true_mean, tolerance = 1.5)
+})
+
+test_that("a Weibull recovery delay compiles and gives sensible moments", {
+  skip_if_no_cmdstan()
+  set.seed(5)
+  ll <- simulate_linelist(
+    n = 2000, cfr = 0.4, onset_days = 40,
+    delay = Gamma(mean = 12.75, sd = 7),
+    recovery = Weibull(mean = 21, sd = 9)
+  )
+  d <- prepare_cfr_data(ll, obs_time = max(ll$onset_date) - 2)
+  fit <- fit_quick(d,
+    delay = Gamma(shape = Normal(3.3, 1), rate = Normal(0.26, 0.08)),
+    recovery_delay = Weibull(shape = Normal(2.4, 0.6), scale = Normal(23.7, 4)),
+    cfr_prior = Beta(1, 1)
+  )
+  expect_equal(fit$cfrnow$recovery_family, "weibull") # differs from gamma death
+  s <- summary(fit)
+  rm <- s[s$quantity == "recovery_mean", ] # onset-to-recovery ~21 (heavy censoring)
+  expect_gt(rm[["q50"]], 17)
+  expect_lt(rm[["q50"]], 26)
+})
+
 test_that("a fixed delay runs the Ghani/Nishiura estimator (delay held constant)", {
   skip_if_no_cmdstan()
   set.seed(4)