Skip to content

Commit

Permalink
Rename obs_at to obs_time and obs_t to relative_obs_time
Browse files Browse the repository at this point in the history
  • Loading branch information
athowes committed Sep 20, 2024
1 parent 7ee66a8 commit 7ffdab0
Show file tree
Hide file tree
Showing 12 changed files with 55 additions and 57 deletions.
2 changes: 1 addition & 1 deletion .github/workflows/check-cmdstan.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -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
)
Expand Down
10 changes: 5 additions & 5 deletions R/latent_gamma.R
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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]

Expand All @@ -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
Expand Down
14 changes: 7 additions & 7 deletions R/latent_individual.R
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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
Expand All @@ -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,
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down
10 changes: 5 additions & 5 deletions R/latent_lognormal.R
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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]

Expand All @@ -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)
Expand Down
32 changes: 16 additions & 16 deletions R/observe.R
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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))
)
}

Expand All @@ -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)
}
6 changes: 3 additions & 3 deletions inst/make_hexsticker.R
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down
4 changes: 2 additions & 2 deletions inst/stan/latent_individual/functions.stan
Original file line number Diff line number Diff line change
Expand Up @@ -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);
}
4 changes: 1 addition & 3 deletions man/filter_obs_by_ptime.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

2 changes: 1 addition & 1 deletion man/observe_process.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

2 changes: 1 addition & 1 deletion tests/testthat/test-latent_gamma.R
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
4 changes: 2 additions & 2 deletions tests/testthat/test-latent_individual.R
Original file line number Diff line number Diff line change
Expand Up @@ -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),
Expand All @@ -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),
Expand Down
22 changes: 11 additions & 11 deletions vignettes/ebola.Rmd
Original file line number Diff line number Diff line change
Expand Up @@ -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 |>
Expand All @@ -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 |>
Expand All @@ -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 |>
Expand Down Expand Up @@ -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)
```

Expand All @@ -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)) +
Expand All @@ -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 |>
Expand All @@ -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 |>
Expand Down Expand Up @@ -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)) +
Expand All @@ -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 |>
Expand All @@ -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 |>
Expand Down

0 comments on commit 7ffdab0

Please sign in to comment.