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 |>