From 8e776d09bc7b8b7f509b7ec58059717c8cf70a40 Mon Sep 17 00:00:00 2001 From: Sebastian Funk Date: Tue, 10 Dec 2024 21:29:29 +0000 Subject: [PATCH 1/6] make D real --- inst/stan/pcens_model.stan | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/inst/stan/pcens_model.stan b/inst/stan/pcens_model.stan index c900041..8471007 100644 --- a/inst/stan/pcens_model.stan +++ b/inst/stan/pcens_model.stan @@ -6,7 +6,7 @@ functions { real partial_sum(array[] int dummy, int start, int end, array[] int d, array[] int d_upper, array[] int n, - array[] int pwindow, array[] int D, + array[] int pwindow, data array[] real D, int dist_id, array[] real params, int primary_id, array[] real primary_params) { real partial_target = 0; @@ -27,7 +27,7 @@ data { array[N] int d_upper; // observed delays upper bound array[N] int n; // number of occurrences for each delay array[N] int pwindow; // primary censoring window - array[N] int D; // maximum delay + array[N] real D; // maximum delay int dist_id; // distribution identifier int primary_id; // primary distribution identifier int n_params; // number of distribution parameters From 52ea1ddd743862886ac699e18cc288144b853802 Mon Sep 17 00:00:00 2001 From: Sebastian Funk Date: Wed, 11 Dec 2024 10:40:49 +0000 Subject: [PATCH 2/6] add regression test --- tests/testthat/test-pcd_cmdstan_model.R | 69 +++++++++++++++++++++++++ 1 file changed, 69 insertions(+) diff --git a/tests/testthat/test-pcd_cmdstan_model.R b/tests/testthat/test-pcd_cmdstan_model.R index 24aa0d7..2a60cff 100644 --- a/tests/testthat/test-pcd_cmdstan_model.R +++ b/tests/testthat/test-pcd_cmdstan_model.R @@ -342,3 +342,72 @@ test_that( ) } ) + +test_that("pcd_cmdstan_model recovers true values with no bound on D", { + # Simulate data + set.seed(123) + n <- 2000 + true_meanlog <- 1.5 + true_sdlog <- 0.5 + + simulated_delays <- rprimarycensored( + n = n, + rdist = rlnorm, + meanlog = true_meanlog, + sdlog = true_sdlog, + pwindow = 1, + D = Inf + ) + + simulated_data <- data.frame( + delay = simulated_delays, + delay_upper = simulated_delays + 1, + pwindow = 1, + relative_obs_time = Inf + ) + + delay_counts <- simulated_data |> + dplyr::summarise( + n = dplyr::n(), + .by = c(pwindow, relative_obs_time, delay, delay_upper) + ) + + # Prepare data for Stan + stan_data <- pcd_as_stan_data( + delay_counts, + dist_id = 1, # Lognormal + primary_id = 1, # Uniform + param_bounds = list(lower = c(-Inf, 0), upper = c(Inf, Inf)), + primary_param_bounds = list(lower = numeric(0), upper = numeric(0)), + priors = list(location = c(0, 1), scale = c(1, 1)), + primary_priors = list(location = numeric(0), scale = numeric(0)) + ) + + # Fit model + model <- suppressMessages(suppressWarnings(pcd_cmdstan_model())) + fit <- suppressMessages(suppressWarnings(model$sample( + data = stan_data, + seed = 123, + chains = 4, + parallel_chains = 4, + refresh = 0, + show_messages = FALSE, + iter_warmup = 500 + ))) + + # Extract posterior + posterior <- fit$draws(c("params[1]", "params[2]"), format = "df") + + # Check mean estimates + expect_equal(mean(posterior$`params[1]`), true_meanlog, tolerance = 0.05) + expect_equal(mean(posterior$`params[2]`), true_sdlog, tolerance = 0.05) + + # Check credible intervals + ci_meanlog <- quantile(posterior$`params[1]`, c(0.05, 0.95)) + ci_sdlog <- quantile(posterior$`params[2]`, c(0.05, 0.95)) + + expect_gt(true_meanlog, ci_meanlog[1]) + expect_lt(true_meanlog, ci_meanlog[2]) + expect_gt(true_sdlog, ci_sdlog[1]) + expect_lt(true_sdlog, ci_sdlog[2]) +}) From 3187f422e81a1353427fa904485d5ab52cb0f0eb Mon Sep 17 00:00:00 2001 From: Sebastian Funk Date: Wed, 11 Dec 2024 10:42:58 +0000 Subject: [PATCH 3/6] add news item --- NEWS.md | 1 + 1 file changed, 1 insertion(+) diff --git a/NEWS.md b/NEWS.md index 03dfe1c..46a0b5c 100644 --- a/NEWS.md +++ b/NEWS.md @@ -11,6 +11,7 @@ Development release. ## Bug fixes - Added a missing `@family` tag to the `pcens` functions. This omission resulted in the Weibull analytical solution not being visible in the package documentation. +- Changed `D` to be of type real in order to support infinite `relative_obs_time`. # primarycensored 1.0.0 From 6e18515ac6e46521549e8affece8436aec7768fd Mon Sep 17 00:00:00 2001 From: Sam Abbott Date: Wed, 11 Dec 2024 13:35:58 +0000 Subject: [PATCH 4/6] Update NEWS.md --- NEWS.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index b51023f..f1dadb0 100644 --- a/NEWS.md +++ b/NEWS.md @@ -12,7 +12,7 @@ Development release. - Added a missing `@family` tag to the `pcens` functions. This omission resulted in the Weibull analytical solution not being visible in the package documentation. - Changed a call to `size()` to use `num_elements()` instead as an underlying type conversion was causing issues on some platforms. -- Changed `D` to be of type real in order to support infinite `relative_obs_time`. +- Changed `D` to be of type real in `pcens_model.stan` in order to support infinite `relative_obs_time`. # primarycensored 1.0.0 From dd7894e7d092870cde1a3133828ef0c1980ff430 Mon Sep 17 00:00:00 2001 From: Sam Abbott Date: Wed, 11 Dec 2024 13:37:10 +0000 Subject: [PATCH 5/6] Update tests/testthat/test-pcd_cmdstan_model.R --- tests/testthat/test-pcd_cmdstan_model.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/testthat/test-pcd_cmdstan_model.R b/tests/testthat/test-pcd_cmdstan_model.R index 2a60cff..636bbe1 100644 --- a/tests/testthat/test-pcd_cmdstan_model.R +++ b/tests/testthat/test-pcd_cmdstan_model.R @@ -388,7 +388,7 @@ test_that("pcd_cmdstan_model recovers true values with no bound on D", { fit <- suppressMessages(suppressWarnings(model$sample( data = stan_data, seed = 123, - chains = 4, + chains = 2, parallel_chains = 4, refresh = 0, show_messages = FALSE, From c59bad8a1e6a08825d01bfd803b4f945762c0fc0 Mon Sep 17 00:00:00 2001 From: Sam Abbott Date: Wed, 11 Dec 2024 13:37:17 +0000 Subject: [PATCH 6/6] Update tests/testthat/test-pcd_cmdstan_model.R --- tests/testthat/test-pcd_cmdstan_model.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/testthat/test-pcd_cmdstan_model.R b/tests/testthat/test-pcd_cmdstan_model.R index 636bbe1..147c31f 100644 --- a/tests/testthat/test-pcd_cmdstan_model.R +++ b/tests/testthat/test-pcd_cmdstan_model.R @@ -389,7 +389,7 @@ test_that("pcd_cmdstan_model recovers true values with no bound on D", { data = stan_data, seed = 123, chains = 2, - parallel_chains = 4, + parallel_chains = 2 refresh = 0, show_messages = FALSE, iter_warmup = 500