Issue 240: Add pointwise log likelihood for latent lognormal model (#259

) * Template for log_lik function * Add wrong but working version of the log_lik * Lint things * Add wN section, link to issue, better documentation * Lint and document * Correction to the overlap case and improve unit test
epinowcast · Sep 2, 2024 · 352e69b · 352e69b
1 parent 993f1ec
commit 352e69b
Show file tree

Hide file tree

Showing 15 changed files with 94 additions and 0 deletions.
diff --git a/NAMESPACE b/NAMESPACE
@@ -42,6 +42,7 @@ export(filter_obs_by_ptime)
 export(is_latent_individual)
 export(linelist_to_cases)
 export(linelist_to_counts)
+export(log_lik_latent_lognormal)
 export(make_relative_to_truth)
 export(observe_process)
 export(pad_zero)

diff --git a/R/latent_lognormal.R b/R/latent_lognormal.R
@@ -47,3 +47,42 @@ posterior_epred_latent_lognormal <- function(prep) { # nolint: object_length_lin
   sigma <- brms::get_dpar(prep, "sigma")
   exp(mu + sigma^2 / 2)
 }
+
+#' Calculate the pointwise log likelihood
+#'
+#' See [brms::log_lik()].
+#'
+#' @param i The index of the observation to calculate the log likelihood of
+#' @param prep The result of a call to [brms::prepare_predictions()]
+#' @family postprocess
+#' @autoglobal
+#' @export
+log_lik_latent_lognormal <- function(i, prep) {
+  mu <- brms::get_dpar(prep, "mu", i = i)
+  sigma <- brms::get_dpar(prep, "sigma", i = i)
+  y <- prep$data$Y[i]
+  obs_t <- prep$data$vreal1[i]
+  pwindow_width <- prep$data$vreal2[i]
+  swindow_width <- prep$data$vreal3[i]
+
+  # Generates values of the swindow_raw and pwindow_raw, but really these should
+  # be extracted from prep or the fitted raws somehow. See:
+  # https://github.com/epinowcast/epidist/issues/267
+  swindow_raw <- runif(prep$ndraws)
+  pwindow_raw <- runif(prep$ndraws)
+
+  swindow <- swindow_raw * swindow_width
+
+  # For no overlap calculate as usual, for overlap ensure pwindow < swindow
+  if (i %in% prep$data$noverlap) {
+    pwindow <- pwindow_raw * pwindow_width
+  } else {
+    pwindow <- pwindow_raw * swindow
+  }
+
+  d <- y - pwindow + swindow
+  obs_time <- obs_t - pwindow
+  lpdf <- dlnorm(d, meanlog = mu, sdlog = sigma, log = TRUE)
+  lcdf <- plnorm(obs_time, meanlog = mu, sdlog = sigma, log.p = TRUE)
+  return(lpdf - lcdf)
+}
diff --git a/man/add_mean_sd.Rd b/man/add_mean_sd.Rd
diff --git a/man/add_mean_sd.default.Rd b/man/add_mean_sd.default.Rd
diff --git a/man/add_mean_sd.gamma_samples.Rd b/man/add_mean_sd.gamma_samples.Rd
diff --git a/man/add_mean_sd.lognormal_samples.Rd b/man/add_mean_sd.lognormal_samples.Rd
diff --git a/man/draws_to_long.Rd b/man/draws_to_long.Rd
diff --git a/man/log_lik_latent_lognormal.Rd b/man/log_lik_latent_lognormal.Rd
diff --git a/man/make_relative_to_truth.Rd b/man/make_relative_to_truth.Rd
diff --git a/man/posterior_epred_latent_lognormal.Rd b/man/posterior_epred_latent_lognormal.Rd
diff --git a/man/posterior_predict_latent_lognormal.Rd b/man/posterior_predict_latent_lognormal.Rd
diff --git a/man/predict_delay_parameters.Rd b/man/predict_delay_parameters.Rd
diff --git a/man/summarise_draws.Rd b/man/summarise_draws.Rd
diff --git a/man/summarise_variable.Rd b/man/summarise_variable.Rd
diff --git a/tests/testthat/test-unit-predict.R b/tests/testthat/test-unit-predict.R
@@ -50,3 +50,15 @@ test_that("posterior_epred_latent_lognormal creates a array of non-negative numb
   expect_equal(ncol(epred), length(prep$data$Y))
   expect_gte(min(epred), 0)
 })
+
+test_that("log_lik_latent_lognormal produces a vector with length ndraws of finite non-NA numbers", { # nolint: line_length_linter.
+  fit <- readRDS(
+    system.file("extdata/fit.rds", package = "epidist")
+  )
+  prep <- brms::prepare_predictions(fit)
+  i <- 1
+  log_lik <- log_lik_latent_lognormal(i, prep)
+  expect_equal(length(log_lik), prep$ndraws)
+  expect_true(all(!is.na(log_lik)))
+  expect_true(all(is.finite(log_lik)))
+})