From 7ffdab0d62363cc92559bd144e627586d81d5f04 Mon Sep 17 00:00:00 2001 From: athowes Date: Fri, 20 Sep 2024 14:30:13 +0100 Subject: [PATCH] Rename obs_at to obs_time and obs_t to relative_obs_time --- .github/workflows/check-cmdstan.yaml | 2 +- R/latent_gamma.R | 10 +++---- R/latent_individual.R | 14 +++++----- R/latent_lognormal.R | 10 +++---- R/observe.R | 32 +++++++++++----------- inst/make_hexsticker.R | 6 ++-- inst/stan/latent_individual/functions.stan | 4 +-- man/filter_obs_by_ptime.Rd | 4 +-- man/observe_process.Rd | 2 +- tests/testthat/test-latent_gamma.R | 2 +- tests/testthat/test-latent_individual.R | 4 +-- vignettes/ebola.Rmd | 22 +++++++-------- 12 files changed, 55 insertions(+), 57 deletions(-) diff --git a/.github/workflows/check-cmdstan.yaml b/.github/workflows/check-cmdstan.yaml index cd6f35e43..6e0feec0f 100644 --- a/.github/workflows/check-cmdstan.yaml +++ b/.github/workflows/check-cmdstan.yaml @@ -61,7 +61,7 @@ jobs: run: | dummy_obs <- dplyr::tibble(case = 1L, ptime = 1, stime = 2, delay_daily = 1, delay_lwr = 1, delay_upr = 2, ptime_lwr = 1, - ptime_upr = 2, stime_lwr = 1, stime_upr = 2, obs_at = 100, + ptime_upr = 2, stime_lwr = 1, stime_upr = 2, obs_time = 100, censored = "interval", censored_obs_time = 10, ptime_daily = 1, stime_daily = 1 ) diff --git a/R/latent_gamma.R b/R/latent_gamma.R index ec3d2aab0..ab6fb9d90 100644 --- a/R/latent_gamma.R +++ b/R/latent_gamma.R @@ -11,14 +11,14 @@ posterior_predict_latent_gamma <- function(i, prep, ...) { # nolint: object_leng mu <- brms::get_dpar(prep, "mu", i = i) shape <- brms::get_dpar(prep, "shape", i = i) - obs_t <- prep$data$vreal1[i] + relative_obs_time <- prep$data$vreal1[i] pwindow <- prep$data$vreal2[i] swindow <- prep$data$vreal3[i] .predict <- function(s) { - d_censored <- obs_t + 1 + d_censored <- relative_obs_time + 1 # while loop to impose the truncation - while (d_censored > obs_t) { + while (d_censored > relative_obs_time) { p_latent <- stats::runif(1, 0, 1) * pwindow d_latent <- stats::rgamma(1, shape = shape[s], scale = mu[s] / shape[s]) s_latent <- p_latent + d_latent @@ -57,7 +57,7 @@ log_lik_latent_gamma <- function(i, prep) { mu <- brms::get_dpar(prep, "mu", i = i) shape <- brms::get_dpar(prep, "shape", i = i) y <- prep$data$Y[i] - obs_t <- prep$data$vreal1[i] + relative_obs_time <- prep$data$vreal1[i] pwindow <- prep$data$vreal2[i] swindow <- prep$data$vreal3[i] @@ -74,7 +74,7 @@ log_lik_latent_gamma <- function(i, prep) { } d <- y - pwindow + swindow - obs_time <- obs_t - pwindow + obs_time <- relative_obs_time - pwindow lpdf <- stats::dgamma(d, shape = shape, scale = mu / shape, log = TRUE) lcdf <- stats::pgamma( obs_time, shape = shape, scale = mu / shape, log.p = TRUE diff --git a/R/latent_individual.R b/R/latent_individual.R index b28312db8..9a4a2341e 100644 --- a/R/latent_individual.R +++ b/R/latent_individual.R @@ -12,7 +12,7 @@ assert_latent_individual_input <- function(data) { assert_names( names(data), must.include = c("case", "ptime_lwr", "ptime_upr", - "stime_lwr", "stime_upr", "obs_at") + "stime_lwr", "stime_upr", "obs_time") ) assert_integer(data$case, lower = 0) assert_numeric(data$ptime_lwr, lower = 0) @@ -21,7 +21,7 @@ assert_latent_individual_input <- function(data) { assert_numeric(data$stime_lwr, lower = 0) assert_numeric(data$stime_upr, lower = 0) assert_true(all(data$stime_upr - data$stime_lwr > 0)) - assert_numeric(data$obs_at, lower = 0) + assert_numeric(data$obs_time, lower = 0) } #' Prepare latent individual model @@ -46,7 +46,7 @@ as_latent_individual.data.frame <- function(data) { class(data) <- c("epidist_latent_individual", class(data)) data <- data |> mutate( - obs_t = .data$obs_at - .data$ptime_lwr, + relative_obs_time = .data$obs_time - .data$ptime_lwr, pwindow = ifelse( stime_lwr < .data$ptime_upr, stime_upr - .data$ptime_lwr, @@ -81,14 +81,14 @@ epidist_validate.epidist_latent_individual <- function(data) { assert_names( names(data), must.include = c("case", "ptime_lwr", "ptime_upr", - "stime_lwr", "stime_upr", "obs_at", - "obs_t", "pwindow", "woverlap", + "stime_lwr", "stime_upr", "obs_time", + "relative_obs_time", "pwindow", "woverlap", "swindow", "delay", "row_id") ) if (nrow(data) > 1) { assert_factor(data$row_id) } - assert_numeric(data$obs_t, lower = 0) + assert_numeric(data$relative_obs_time, lower = 0) assert_numeric(data$pwindow, lower = 0) assert_numeric(data$woverlap, lower = 0) assert_numeric(data$swindow, lower = 0) @@ -156,7 +156,7 @@ epidist_formula.epidist_latent_individual <- function(data, family, formula, formula <- brms:::validate_formula(formula, data = data, family = family) formula <- stats::update( - formula, delay | vreal(obs_t, pwindow, swindow) ~ . + formula, delay | vreal(relative_obs_time, pwindow, swindow) ~ . ) # Using this here for checking purposes diff --git a/R/latent_lognormal.R b/R/latent_lognormal.R index b098656cc..a9de72c26 100644 --- a/R/latent_lognormal.R +++ b/R/latent_lognormal.R @@ -12,14 +12,14 @@ posterior_predict_latent_lognormal <- function(i, prep, ...) { # nolint: object_ mu <- brms::get_dpar(prep, "mu", i = i) sigma <- brms::get_dpar(prep, "sigma", i = i) - obs_t <- prep$data$vreal1[i] + relative_obs_time <- prep$data$vreal1[i] pwindow <- prep$data$vreal2[i] swindow <- prep$data$vreal3[i] .predict <- function(s) { - d_censored <- obs_t + 1 + d_censored <- relative_obs_time + 1 # while loop to impose the truncation - while (d_censored > obs_t) { + while (d_censored > relative_obs_time) { p_latent <- stats::runif(1, 0, 1) * pwindow d_latent <- stats::rlnorm(1, meanlog = mu[s], sdlog = sigma[s]) s_latent <- p_latent + d_latent @@ -60,7 +60,7 @@ 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] + relative_obs_time <- prep$data$vreal1[i] pwindow <- prep$data$vreal2[i] swindow <- prep$data$vreal3[i] @@ -80,7 +80,7 @@ log_lik_latent_lognormal <- function(i, prep) { } d <- y - pwindow + swindow - obs_time <- obs_t - pwindow + obs_time <- relative_obs_time - pwindow lpdf <- stats::dlnorm(d, meanlog = mu, sdlog = sigma, log = TRUE) lcdf <- stats::plnorm(obs_time, meanlog = mu, sdlog = sigma, log.p = TRUE) return(lpdf - lcdf) diff --git a/R/observe.R b/R/observe.R index cdb5ef38d..0d321a575 100644 --- a/R/observe.R +++ b/R/observe.R @@ -12,7 +12,7 @@ #' * `delay_daily`: Given by `stime_daily - ptime_daily` #' * `delay_lwr`: Given by `delay_daily - 1` (or 0 if `delay_daily < 1`) #' * `delay_upr`: Given by `delay_daily + 1` -#' * `obs_at`: The maximum value of `stime` +#' * `obs_time`: The maximum value of `stime` #' #' @param linelist ... #' @family observe @@ -30,7 +30,7 @@ observe_process <- function(linelist) { delay_daily = .data$stime_daily - .data$ptime_daily, delay_lwr = purrr::map_dbl(.data$delay_daily, ~ max(0, . - 1)), delay_upr = .data$delay_daily + 1, - obs_at = ceiling(max(.data$stime)) + obs_time = ceiling(max(.data$stime)) ) } @@ -44,45 +44,45 @@ observe_process <- function(linelist) { filter_obs_by_obs_time <- function(linelist, obs_time) { linelist |> mutate( - obs_at = obs_time, obs_time = obs_time - .data$ptime, - censored_obs_time = .data$obs_at - .data$ptime_lwr, + censored_obs_time = .data$obs_time - .data$ptime_lwr, censored = "interval" ) |> - filter(.data$stime_upr <= .data$obs_at) + filter(.data$stime_upr <= .data$obs_time) } #' Filter observations based on the observation time of primary events #' #' @param linelist ... #' @param obs_time ... -#' @param obs_at ... +#' @param obs_time ... #' @family observe #' @autoglobal #' @export filter_obs_by_ptime <- function(linelist, obs_time, - obs_at = c("obs_secondary", "max_secondary")) { - obs_at <- match.arg(obs_at) + obs_time_type = + c("obs_secondary", "max_secondary")) { + obs_time <- match.arg(obs_time) pfilt_t <- obs_time truncated_linelist <- linelist |> mutate(censored = "interval") |> filter(.data$ptime_upr <= pfilt_t) - if (obs_at == "obs_secondary") { + if (obs_time_type == "obs_secondary") { # Update observation time to be the same as the maximum secondary time - truncated_linelist <- mutate(truncated_linelist, obs_at = .data$stime_upr) - } else if (obs_at == "max_secondary") { + truncated_linelist <- mutate(truncated_linelist, obs_time = .data$stime_upr) + } else if (obs_time_type == "max_secondary") { truncated_linelist <- truncated_linelist |> - mutate(obs_at := .data$stime_upr |> max() |> ceiling()) + mutate(obs_time := .data$stime_upr |> max() |> ceiling()) } # Make observation time as specified truncated_linelist <- truncated_linelist |> mutate( - obs_time = .data$obs_at - .data$ptime, - censored_obs_time = .data$obs_at - .data$ptime_lwr + obs_time = .data$obs_time - .data$ptime, + censored_obs_time = .data$obs_time - .data$ptime_lwr ) # Set observation time to artificial observation time if needed - if (obs_at == "obs_secondary") { - truncated_linelist <- mutate(truncated_linelist, obs_at = pfilt_t) + if (obs_time_type == "obs_secondary") { + truncated_linelist <- mutate(truncated_linelist, obs_time = pfilt_t) } return(truncated_linelist) } diff --git a/inst/make_hexsticker.R b/inst/make_hexsticker.R index 745b28c43..aaf9dadb8 100644 --- a/inst/make_hexsticker.R +++ b/inst/make_hexsticker.R @@ -26,16 +26,16 @@ truncated_obs <- obs |> combined_obs <- bind_rows( truncated_obs, - mutate(obs, obs_at = max(stime_daily)) + mutate(obs, obs_time = max(stime_daily)) ) |> - mutate(obs_at = factor(obs_at)) + mutate(obs_time = factor(obs_time)) meanlog <- secondary_dist$mu[[1]] sdlog <- secondary_dist$sigma[[1]] hex_plot <- combined_obs |> ggplot() + - aes(x = delay_daily, fill = obs_at) + + aes(x = delay_daily, fill = obs_time) + geom_histogram( aes(y = after_stat(density)), binwidth = 1, position = "dodge" diff --git a/inst/stan/latent_individual/functions.stan b/inst/stan/latent_individual/functions.stan index 43ca3de3d..6aaa1dce6 100644 --- a/inst/stan/latent_individual/functions.stan +++ b/inst/stan/latent_individual/functions.stan @@ -5,9 +5,9 @@ // are/have been replaced using regex real latent_family_lpdf(vector y, dpars_A, vector pwindow, - vector swindow, array[] real obs_t) { + vector swindow, array[] real relative_obs_t) { int n = num_elements(y); vector[n] d = y - pwindow + swindow; - vector[n] obs_time = to_vector(obs_t) - pwindow; + vector[n] obs_time = to_vector(relative_obs_t) - pwindow; return family_lpdf(d | dpars_B) - family_lcdf(obs_time | dpars_B); } diff --git a/man/filter_obs_by_ptime.Rd b/man/filter_obs_by_ptime.Rd index 20b7484a4..9aec0ac44 100644 --- a/man/filter_obs_by_ptime.Rd +++ b/man/filter_obs_by_ptime.Rd @@ -7,15 +7,13 @@ filter_obs_by_ptime( linelist, obs_time, - obs_at = c("obs_secondary", "max_secondary") + obs_time_type = c("obs_secondary", "max_secondary") ) } \arguments{ \item{linelist}{...} \item{obs_time}{...} - -\item{obs_at}{...} } \description{ Filter observations based on the observation time of primary events diff --git a/man/observe_process.Rd b/man/observe_process.Rd index 6dfaf168b..bf7a8bde7 100644 --- a/man/observe_process.Rd +++ b/man/observe_process.Rd @@ -23,7 +23,7 @@ truncation. The columns added are: \item \code{delay_daily}: Given by \code{stime_daily - ptime_daily} \item \code{delay_lwr}: Given by \code{delay_daily - 1} (or 0 if \code{delay_daily < 1}) \item \code{delay_upr}: Given by \code{delay_daily + 1} -\item \code{obs_at}: The maximum value of \code{stime} +\item \code{obs_time}: The maximum value of \code{stime} } } \seealso{ diff --git a/tests/testthat/test-latent_gamma.R b/tests/testthat/test-latent_gamma.R index d567de880..0a51dd69d 100644 --- a/tests/testthat/test-latent_gamma.R +++ b/tests/testthat/test-latent_gamma.R @@ -26,7 +26,7 @@ test_that("posterior_predict_latent_gamma can generate predictions with no censo fit_gamma <- readRDS( system.file("extdata/fit_gamma.rds", package = "epidist") ) - draws <- data.frame(obs_t = 1000, pwindow = 0, swindow = 0) |> + draws <- data.frame(relative_obs_time = 1000, pwindow = 0, swindow = 0) |> tidybayes::add_predicted_draws(fit_gamma, ndraws = 100) expect_equal(draws$.draw, 1:100) pred <- draws$.prediction diff --git a/tests/testthat/test-latent_individual.R b/tests/testthat/test-latent_individual.R index 910f73c80..59e7991e0 100644 --- a/tests/testthat/test-latent_individual.R +++ b/tests/testthat/test-latent_individual.R @@ -83,7 +83,7 @@ test_that("epidist_formula.epidist_latent_individual with default settings produ expect_s3_class(form, "brmsformula") expect_equal( as_string_formula(form$formula), - "delay | vreal(obs_t, pwindow, swindow) ~ 1" + "delay | vreal(relative_obs_time, pwindow, swindow) ~ 1" ) expect_equal( as_string_formula(form$pforms$sigma), @@ -101,7 +101,7 @@ test_that("epidist_formula.epidist_latent_individual with custom formulas produc expect_s3_class(form_sex, "brmsformula") expect_equal( as_string_formula(form_sex$formula), - "delay | vreal(obs_t, pwindow, swindow) ~ sex" + "delay | vreal(relative_obs_time, pwindow, swindow) ~ sex" ) expect_equal( as_string_formula(form_sex$pforms$sigma), diff --git a/vignettes/ebola.Rmd b/vignettes/ebola.Rmd index d391bbf75..2b09a1744 100644 --- a/vignettes/ebola.Rmd +++ b/vignettes/ebola.Rmd @@ -300,7 +300,7 @@ Figure \@ref(fig:epred)B illustrates the higher mean of men as compared with wom epred_draws <- obs_prep |> as.data.frame() |> data_grid(NA) |> - mutate(obs_t = NA, pwindow = NA, swindow = NA) |> + mutate(relative_obs_time = NA, pwindow = NA, swindow = NA) |> add_epred_draws(fit, dpar = TRUE) epred_base_figure <- epred_draws |> @@ -312,7 +312,7 @@ epred_base_figure <- epred_draws |> epred_draws_sex <- obs_prep |> as.data.frame() |> data_grid(sex) |> - mutate(obs_t = NA, pwindow = NA, swindow = NA) |> + mutate(relative_obs_time = NA, pwindow = NA, swindow = NA) |> add_epred_draws(fit_sex, dpar = TRUE) epred_sex_figure <- epred_draws_sex |> @@ -324,7 +324,7 @@ epred_sex_figure <- epred_draws_sex |> epred_draws_sex_district <- obs_prep |> as.data.frame() |> data_grid(sex, district) |> - mutate(obs_t = NA, pwindow = NA, swindow = NA) |> + mutate(relative_obs_time = NA, pwindow = NA, swindow = NA) |> add_epred_draws(fit_sex_district, dpar = TRUE) epred_sex_district_figure <- epred_draws_sex_district |> @@ -355,7 +355,7 @@ For example, for the `mu` parameter in the sex-district stratified model (Figure linpred_draws_sex_district <- obs_prep |> as.data.frame() |> data_grid(sex, district) |> - mutate(obs_t = NA, pwindow = NA, swindow = NA) |> + mutate(relative_obs_time = NA, pwindow = NA, swindow = NA) |> add_linpred_draws(fit_sex_district, dpar = TRUE) ``` @@ -379,13 +379,13 @@ In this section, we demonstrate how to produce either a discrete probability mas ### Discrete probability mass function To generate a discrete probability mass function (PMF) we predict the delay distribution that would be observed with daily censoring and no right truncation. -To do this, we set each of `pwindow` and `swindow` to 1 for daily censoring, and `obs_t` to 1000 for no censoring. +To do this, we set each of `pwindow` and `swindow` to 1 for daily censoring, and `relative_obs_time` to 1000 for no censoring. Figure \@ref(fig:pmf) shows the result, where the few delays greater than 30 are omitted from the figure. ```{r} draws_pmf <- obs_prep |> as.data.frame() |> - mutate(obs_t = 1000, pwindow = 1, swindow = 1) |> + mutate(relative_obs_time = 1000, pwindow = 1, swindow = 1) |> add_predicted_draws(fit, ndraws = 1000) pmf_base_figure <- ggplot(draws_pmf, aes(x = .prediction)) + @@ -397,7 +397,7 @@ pmf_base_figure <- ggplot(draws_pmf, aes(x = .prediction)) + draws_sex_pmf <- obs_prep |> as.data.frame() |> data_grid(sex) |> - mutate(obs_t = 1000, pwindow = 1, swindow = 1) |> + mutate(relative_obs_time = 1000, pwindow = 1, swindow = 1) |> add_predicted_draws(fit_sex, ndraws = 1000) pmf_sex_figure <- draws_sex_pmf |> @@ -411,7 +411,7 @@ pmf_sex_figure <- draws_sex_pmf |> draws_sex_district_pmf <- obs_prep |> as.data.frame() |> data_grid(sex, district) |> - mutate(obs_t = 1000, pwindow = 1, swindow = 1) |> + mutate(relative_obs_time = 1000, pwindow = 1, swindow = 1) |> add_predicted_draws(fit_sex_district, ndraws = 1000) pmf_sex_district_figure <- draws_sex_district_pmf |> @@ -449,7 +449,7 @@ That is to produce continuous delay times (Figure \@ref(fig:pdf)): ```{r} draws_pdf <- obs_prep |> as.data.frame() |> - mutate(obs_t = 1000, pwindow = 0, swindow = 0) |> + mutate(relative_obs_time = 1000, pwindow = 0, swindow = 0) |> add_predicted_draws(fit, ndraws = 1000) pdf_base_figure <- ggplot(draws_pdf, aes(x = .prediction)) + @@ -461,7 +461,7 @@ pdf_base_figure <- ggplot(draws_pdf, aes(x = .prediction)) + draws_sex_pdf <- obs_prep |> as.data.frame() |> data_grid(sex) |> - mutate(obs_t = 1000, pwindow = 0, swindow = 0) |> + mutate(relative_obs_time = 1000, pwindow = 0, swindow = 0) |> add_predicted_draws(fit_sex, ndraws = 1000) pdf_sex_figure <- draws_sex_pdf |> @@ -475,7 +475,7 @@ pdf_sex_figure <- draws_sex_pdf |> draws_sex_district_pdf <- obs_prep |> as.data.frame() |> data_grid(sex, district) |> - mutate(obs_t = 1000, pwindow = 0, swindow = 0) |> + mutate(relative_obs_time = 1000, pwindow = 0, swindow = 0) |> add_predicted_draws(fit_sex_district, ndraws = 1000) pdf_sex_district_figure <- draws_sex_district_pdf |>