Feature: exponential delays

NEWS.md

@@ -12,6 +12,8 @@ This is a major release and contains breaking changes. If needing the old interf
 
 ## Model
 
+* Added support for parametric exponential delay distributions (note that this is comparable to an intercept only non-parametric hazard model). See [#84](https://github.com/epiforecasts/epinowcast/pull/84) by [@seabbs](https://github.com/seabbs).
+
 ## Internals
 
 * Array declarations in the stan model have been updated. To maintain compatibility with [expose_stan_fns()] (which itself depends on `rstan`) additional functionality has been added to parse stan code in this function. See [#74](https://github.com/epiforecasts/epinowcast/issues/74) by [@sbfnk](https://github.com/sbfnk) and [@seabbs](https://github.com/seabbs).

R/model-tools.R

@@ -47,7 +47,7 @@ enw_formula_as_data_list <- function(data, reference_effects, report_effects) {
 #'
 #' @param distribution Character string indicating the type of distribution to
 #' use for reference date effects. The default is to use a lognormal but other
-#' options available include: gamma distributed ("gamma").
+#' options available include the, exponential and gamma distributions. 
 #'
 #' @param nowcast Logical, defaults to `TRUE`. Should a nowcast be made using
 #' posterior predictions of the unobserved future reported notifications.
@@ -76,10 +76,13 @@ enw_opts_as_data_list <- function(distribution = "lognormal", nowcast = TRUE,
     nowcast <- TRUE
   }
   # check distribution type is supported and change to numeric
-  distribution <- match.arg(distribution, c("lognormal", "gamma"))
+  distribution <- match.arg(
+    distribution, c("exponential", "lognormal", "gamma")
+  )
   distribution <- data.table::fcase(
-    distribution %in% "lognormal", 0,
-    distribution %in% "gamma", 1
+    distribution %in% "exponential", 0,
+    distribution %in% "lognormal", 1,
+    distribution %in% "gamma", 2
   )
 
   data <- list(

R/model.R

@@ -179,8 +179,17 @@ enw_as_data_list <- function(pobs,
 enw_inits <- function(data) {
   init_fn <- function() {
     init <- list(
-      logmean_int = rnorm(1, data$logmean_int_p[1], data$logmean_int_p[2] / 10),
-      logsd_int = abs(rnorm(1, data$logsd_int_p[1], data$logsd_int_p[2] / 10)),
+      logmean_int = rnorm(1, data$logmean_int_p[1], data$logmean_int_p[2] / 10)
+    )
+    if (data$dist > 0) {
+      init$logsd_int <- abs(
+        rnorm(1, data$logsd_int_p[1], data$logsd_int_p[2] / 10)
+      )
+    }else{
+      init$logsd_int <- numeric(0)
+    }
+
+    init <- c(init, list(
       leobs_init = array(purrr::map_dbl(
         data$latest_obs[1, ] + 1,
         ~ rnorm(1, log(.), 1)
@@ -193,7 +202,7 @@ enw_inits <- function(data) {
         dim = c(data$t - 1, data$g)
       ),
       sqrt_phi = abs(rnorm(1, data$sqrt_phi_p[1], data$sqrt_phi_p[2] / 10))
-    )
+    ))
     init$logmean <- rep(init$logmean_int, data$npmfs)
     init$logsd <- rep(init$logsd_int, data$npmfs)
     init$phi <- 1 / sqrt(init$sqrt_phi)

inst/stan/chunks/debug.stan

@@ -12,10 +12,14 @@
         print(pmfs[, i]);
         print("Logmean and Logsd intercept");
         print(logmean_int);
-        print(logsd_int);
+        if (dist) {
+          print(logsd_int);
+        }
         print("Logmean and Logsd for pmf");
         print(logmean[i]);
-        print(logsd[i]);
+        if (dist) {
+          print(logsd[i]);
+        }
         print("Unique report day hazards");
         print(srdlh);
         print("Overdispersion");

inst/stan/epinowcast.stan

@@ -68,13 +68,13 @@ parameters {
   array[g] real leobs_init; // First time point for expected observations
   array[g] real<lower=0> eobs_lsd; // standard deviation of rw for primary obs
   array[g] vector[t - 1] leobs_resids; // unscaled rw for primary obs
-  real<lower=-10, upper=logdmax> logmean_int; // logmean intercept
-  real<lower=1e-3, upper=2*dmax> logsd_int; // logsd intercept
+  array[1] real<lower=-10, upper=logdmax> logmean_int; // logmean intercept
+  array[dist ? 1 : 0]real<lower=1e-3, upper=2*dmax> logsd_int; // logsd intercept
   vector[neffs] logmean_eff; // unscaled modifiers to log mean
-  vector[neffs] logsd_eff; // unscaled modifiers to log sd
+  vector[dist ? neffs : 0] logsd_eff; // unscaled modifiers to log sd
   vector[nrd_effs] rd_eff; // unscaled modifiers to report date hazard
   vector<lower=0>[neff_sds] logmean_sd; // pooled modifiers to logmean
-  vector<lower=0>[neff_sds] logsd_sd; // pooled modifiers to logsd
+  vector<lower=0>[dist ? neff_sds : 0] logsd_sd; // pooled modifiers to logsd
   vector<lower=0>[nrd_eff_sds] rd_eff_sd; // pooled modifiers to report date
   real<lower=0, upper=1e4> sqrt_phi; // Overall dispersion by group
 }
@@ -90,11 +90,13 @@ transformed parameters{
   // calculate log mean and sd parameters for each dataset from design matrices
   profile("transformed_delay_reference_date_total") {
   profile("transformed_delay_reference_date_effects") {
-  logmean = combine_effects(logmean_int, logmean_eff, d_fixed, logmean_sd,
+  logmean = combine_effects(logmean_int[1], logmean_eff, d_fixed, logmean_sd,
                             d_random);
-  logsd = combine_effects(log(logsd_int), logsd_eff, d_fixed, logsd_sd,
-                          d_random); 
-  logsd = exp(logsd);
+  if (dist) {
+    logsd = combine_effects(log(logsd_int[1]), logsd_eff, d_fixed, logsd_sd,
+                            d_random); 
+    logsd = exp(logsd);
+  }
   }
   // calculate pmfs
   profile("transformed_delay_reference_date_pmfs") {
@@ -145,14 +147,20 @@ model {
   }
   // priors for the intercept of the log normal truncation distribution
   logmean_int ~ normal(logmean_int_p[1], logmean_int_p[2]);
-  logsd_int ~ normal(logsd_int_p[1], logsd_int_p[2]);
+  if (dist) {
+    logsd_int ~ normal(logsd_int_p[1], logsd_int_p[2]);
+  }
   // priors and scaling for date of reference effects
   if (neffs) {
     logmean_eff ~ std_normal();
-    logsd_eff ~ std_normal();
+    if (dist) {
+      logsd_eff ~ std_normal();
+    }
     if (neff_sds) {
       logmean_sd ~ zero_truncated_normal(logmean_sd_p[1], logmean_sd_p[2]);
-      logsd_sd ~ zero_truncated_normal(logsd_sd_p[1], logsd_sd_p[2]);
+      if (dist) {
+        logsd_sd ~ zero_truncated_normal(logsd_sd_p[1], logsd_sd_p[2]);
+      }
     }
   }
   // priors and scaling for date of report effects

inst/stan/functions/discretised_reporting_prob.stan

@@ -4,10 +4,15 @@ vector discretised_reporting_prob(real mu, real sigma, int n, int dist) {
   vector[n] pmf; 
   vector[n] upper_cdf;
   if (dist == 0) {
+    real emu = exp(-mu);
     for (i in 1:n) {
-      upper_cdf[i] = lognormal_cdf(i | mu, sigma);
+      upper_cdf[i] = exponential_cdf(i | emu);
     }
   } else if (dist == 1) {
+    for (i in 1:n) {
+      upper_cdf[i] = lognormal_cdf(i | mu, sigma);
+    }
+  } else if (dist == 2) {
     real emu = exp(mu);
     for (i in 1:n) {
       upper_cdf[i] = gamma_cdf(i | emu, sigma);