From ab2b589d4ae933c6f13b408bb31ffefcdff8d690 Mon Sep 17 00:00:00 2001 From: gowerc Date: Tue, 25 Jun 2024 16:23:29 +0100 Subject: [PATCH 1/8] Initital attempt --- NAMESPACE | 1 + R/JointModel.R | 7 ++-- R/JointModelSamples.R | 34 ++++++++++++++----- R/generics.R | 5 +-- R/utilities.R | 5 +++ inst/stan/lm-stein-fojo/quantities.stan | 6 ++-- man/as.StanModule.JointModelSamples.Rd | 24 +++++++++++++ man/write_stan.Rd | 8 +++-- tests/testthat/test-LongitudinalGSF.R | 23 +++++++++++++ tests/testthat/test-LongitudinalRandomSlope.R | 23 +++++++++++++ tests/testthat/test-LongitudinalSteinFojo.R | 25 +++++++++++++- vignettes/model_fitting.Rmd | 3 +- 12 files changed, 143 insertions(+), 21 deletions(-) create mode 100644 man/as.StanModule.JointModelSamples.Rd diff --git a/NAMESPACE b/NAMESPACE index 0f2bebee..856533ea 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -18,6 +18,7 @@ S3method(as.QuantityGenerator,GridObserved) S3method(as.QuantityGenerator,GridPopulation) S3method(as.QuantityGenerator,GridPrediction) S3method(as.StanModule,JointModel) +S3method(as.StanModule,JointModelSamples) S3method(as.StanModule,Link) S3method(as.StanModule,LinkComponent) S3method(as.StanModule,Parameter) diff --git a/R/JointModel.R b/R/JointModel.R index 183f130b..aa91d0ff 100755 --- a/R/JointModel.R +++ b/R/JointModel.R @@ -142,8 +142,11 @@ as.character.JointModel <- function(x, ...) { #' @rdname write_stan #' @export -write_stan.JointModel <- function(object, file_path) { - fi <- file(file_path, open = "w") +write_stan.JointModel <- function(object, destination) { + if (is_connection(destination)) { + return(writeLines(as.character(object), fi)) + } + fi <- file(destination, open = "w") writeLines(as.character(object), fi) close(fi) } diff --git a/R/JointModelSamples.R b/R/JointModelSamples.R index 76a1f5b4..c0df8311 100644 --- a/R/JointModelSamples.R +++ b/R/JointModelSamples.R @@ -38,7 +38,30 @@ generateQuantities.JointModelSamples <- function(object, generator, type, ...) { append(as_stan_list(object@model@parameters)) |> append(as_stan_list(generator, data = object@data, model = object@model)) + stanobj <- as.StanModule(object, generator = generator, type = type) + model <- compileStanModel(stanobj) + + devnull <- utils::capture.output( + results <- model$generate_quantities( + data = data, + fitted_params = object@results + ) + ) + return(results) +} + + +#' `JointModelSamples` -> `StanModule` +#' +#' Converts a `JointModelSamples` object into a `StanModule` object ensuring +#' that the resulting `StanModule` object is able to generate post sampling +#' quantities. +#' +#' @inheritParams generateQuantities +#' @export +as.StanModule.JointModelSamples <- function(object, generator, type, ...) { assert_that( + is(generator, "QuantityGenerator"), length(type) == 1, type %in% c("survival", "longitudinal") ) @@ -60,17 +83,10 @@ generateQuantities.JointModelSamples <- function(object, generator, type, ...) { quant_stanobj ) ) + stanobj +} - model <- compileStanModel(stanobj) - devnull <- utils::capture.output( - results <- model$generate_quantities( - data = data, - fitted_params = object@results - ) - ) - return(results) -} #' `JointModelSamples` -> Printable `Character` #' diff --git a/R/generics.R b/R/generics.R index ec0e45ee..90d305ee 100755 --- a/R/generics.R +++ b/R/generics.R @@ -44,10 +44,11 @@ NULL #' Write the Stan code for a Stan module. #' #' @param object the module. -#' @param file_path (`string`)\cr output file. +#' @param destination (`character` or `connection`)\cr Where to write stan code to. +#' @param ... Additional arguements #' #' @export -write_stan <- function(object, file_path) { +write_stan <- function(object, destination, ...) { UseMethod("write_stan") } diff --git a/R/utilities.R b/R/utilities.R index 03958256..87b9de8a 100644 --- a/R/utilities.R +++ b/R/utilities.R @@ -329,3 +329,8 @@ decompose_subjects <- function(subjects, all_subjects) { is_cmdstanr_available <- function() { requireNamespace("cmdstanr", quietly = TRUE) } + + +is_connection <- function(obj) { + inherits(obj, "connection") +} diff --git a/inst/stan/lm-stein-fojo/quantities.stan b/inst/stan/lm-stein-fojo/quantities.stan index 92c07fdb..7352540c 100644 --- a/inst/stan/lm-stein-fojo/quantities.stan +++ b/inst/stan/lm-stein-fojo/quantities.stan @@ -15,9 +15,9 @@ generated quantities { // Source - lm-stein-fojo/quantities.stan // matrix[n_subjects, 3] long_gq_parameters; - long_gq_parameters[, 1] = lm_gsf_psi_bsld; - long_gq_parameters[, 2] = lm_gsf_psi_ks; - long_gq_parameters[, 3] = lm_gsf_psi_kg; + long_gq_parameters[, 1] = lm_sf_psi_bsld; + long_gq_parameters[, 2] = lm_sf_psi_ks; + long_gq_parameters[, 3] = lm_sf_psi_kg; matrix[gq_n_quant, 3] long_gq_pop_parameters; long_gq_pop_parameters[, 1] = exp(lm_sf_mu_bsld[gq_long_pop_study_index]); diff --git a/man/as.StanModule.JointModelSamples.Rd b/man/as.StanModule.JointModelSamples.Rd new file mode 100644 index 00000000..fe70d01b --- /dev/null +++ b/man/as.StanModule.JointModelSamples.Rd @@ -0,0 +1,24 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/JointModelSamples.R +\name{as.StanModule.JointModelSamples} +\alias{as.StanModule.JointModelSamples} +\title{\code{JointModelSamples} -> \code{StanModule}} +\usage{ +\method{as.StanModule}{JointModelSamples}(object, generator, type, ...) +} +\arguments{ +\item{object}{object to obtain generated quantities from} + +\item{generator}{(\code{QuantityGenerator})\cr object that specifies which subjects and time points +to calculate the quantities at} + +\item{type}{(\code{character})\cr type of quantities to be generated, must be either "survival" or +"longitudinal".} + +\item{...}{additional options.} +} +\description{ +Converts a \code{JointModelSamples} object into a \code{StanModule} object ensuring +that the resulting \code{StanModule} object is able to generate post sampling +quantities. +} diff --git a/man/write_stan.Rd b/man/write_stan.Rd index 5b64623b..7cc055e8 100644 --- a/man/write_stan.Rd +++ b/man/write_stan.Rd @@ -5,14 +5,16 @@ \alias{write_stan.JointModel} \title{\code{write_stan}} \usage{ -write_stan(object, file_path) +write_stan(object, destination, ...) -\method{write_stan}{JointModel}(object, file_path) +\method{write_stan}{JointModel}(object, destination) } \arguments{ \item{object}{the module.} -\item{file_path}{(\code{string})\cr output file.} +\item{destination}{(\code{character} or \code{connection})\cr Where to write stan code to.} + +\item{...}{Additional arguements} } \description{ Write the Stan code for a Stan module. diff --git a/tests/testthat/test-LongitudinalGSF.R b/tests/testthat/test-LongitudinalGSF.R index 571f8a48..fb6c56bf 100644 --- a/tests/testthat/test-LongitudinalGSF.R +++ b/tests/testthat/test-LongitudinalGSF.R @@ -197,3 +197,26 @@ test_that("Can recover known distributional parameters from a full GSF joint mod expect_true(all(dat$q99 >= true_values)) expect_true(all(dat$ess_bulk > 100)) }) + + + +test_that("Quantity models pass the parser", { + mock_samples <- .JointModelSamples( + model = JointModel(longitudinal = LongitudinalGSF(centred = FALSE)), + data = structure(1, class = "DataJoint"), + results = structure(1, class = "CmdStanMCMC") + ) + stanmod <- as.StanModule( + mock_samples, + generator = QuantityGeneratorPopulation(1, "A", "B"), + type = "longitudinal" + ) + expect_stan_syntax(stanmod) + + stanmod <- as.StanModule( + mock_samples, + generator = QuantityGeneratorSubject(1, "A"), + type = "longitudinal" + ) + expect_stan_syntax(stanmod) +}) diff --git a/tests/testthat/test-LongitudinalRandomSlope.R b/tests/testthat/test-LongitudinalRandomSlope.R index ff60e437..ca151791 100644 --- a/tests/testthat/test-LongitudinalRandomSlope.R +++ b/tests/testthat/test-LongitudinalRandomSlope.R @@ -295,3 +295,26 @@ test_that("Random Slope Model left-censoring works as expected", { ) expect_gt(lmer_cor, 0.99) }) + + +test_that("Quantity models pass the parser", { + mock_samples <- .JointModelSamples( + model = JointModel(longitudinal = LongitudinalRandomSlope()), + data = structure(1, class = "DataJoint"), + results = structure(1, class = "CmdStanMCMC") + ) + stanmod <- as.StanModule( + mock_samples, + generator = QuantityGeneratorPopulation(1, "A", "B"), + type = "longitudinal" + ) + expect_stan_syntax(stanmod) + + stanmod <- as.StanModule( + mock_samples, + generator = QuantityGeneratorSubject(1, "A"), + type = "longitudinal" + ) + expect_stan_syntax(stanmod) +}) + diff --git a/tests/testthat/test-LongitudinalSteinFojo.R b/tests/testthat/test-LongitudinalSteinFojo.R index 9af17236..03d42b09 100644 --- a/tests/testthat/test-LongitudinalSteinFojo.R +++ b/tests/testthat/test-LongitudinalSteinFojo.R @@ -1,5 +1,4 @@ - test_that("LongitudinalSteinFojo works as expected with default arguments", { result <- expect_silent(LongitudinalSteinFojo()) expect_s4_class(result, "LongitudinalSteinFojo") @@ -90,6 +89,7 @@ test_that("Non-Centralised parameterisation compiles without issues", { + test_that("Can recover known distributional parameters from a SF joint model", { skip_if_not(is_full_test()) @@ -388,3 +388,26 @@ test_that("Can recover known distributional parameters from a SF joint model wit expect_true(all(dat$q99 >= true_values)) expect_true(all(dat$ess_bulk > 100)) }) + + + +test_that("Quantity models pass the parser", { + mock_samples <- .JointModelSamples( + model = JointModel(longitudinal = LongitudinalSteinFojo(centred = TRUE)), + data = structure(1, class = "DataJoint"), + results = structure(1, class = "CmdStanMCMC") + ) + stanmod <- as.StanModule( + mock_samples, + generator = QuantityGeneratorPopulation(1, "A", "B"), + type = "longitudinal" + ) + expect_stan_syntax(stanmod) + + stanmod <- as.StanModule( + mock_samples, + generator = QuantityGeneratorSubject(1, "A"), + type = "longitudinal" + ) + expect_stan_syntax(stanmod) +}) diff --git a/vignettes/model_fitting.Rmd b/vignettes/model_fitting.Rmd index 9c12088e..de50f7dd 100644 --- a/vignettes/model_fitting.Rmd +++ b/vignettes/model_fitting.Rmd @@ -237,9 +237,10 @@ It is always possible to read out the Stan code that is contained in the ```{r debug_stan} tmp <- tempfile() -write_stan(simple_model, file_path = tmp) +write_stan(simple_model, destination = tmp) first_part <- head(readLines(tmp), 10) cat(paste(first_part, collapse = "\n")) +close(tmp) ``` ## Sampling Parameters From 0ca508eeb51f4476f151eec50727b9c8912e128b Mon Sep 17 00:00:00 2001 From: gowerc Date: Tue, 25 Jun 2024 16:34:23 +0100 Subject: [PATCH 2/8] fix cicd failures --- R/generics.R | 2 +- _pkgdown.yml | 1 + man/write_stan.Rd | 2 +- 3 files changed, 3 insertions(+), 2 deletions(-) diff --git a/R/generics.R b/R/generics.R index 90d305ee..cb5082ed 100755 --- a/R/generics.R +++ b/R/generics.R @@ -45,7 +45,7 @@ NULL #' #' @param object the module. #' @param destination (`character` or `connection`)\cr Where to write stan code to. -#' @param ... Additional arguements +#' @param ... Additional arguments #' #' @export write_stan <- function(object, destination, ...) { diff --git a/_pkgdown.yml b/_pkgdown.yml index 868e3a63..9af2868b 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -141,6 +141,7 @@ reference: - as.StanModule.JointModel - as.StanModule.Link - as.StanModule.LinkComponent + - as.StanModule.JointModelSamples - as.character.Parameter - as.character.Prior - show-object diff --git a/man/write_stan.Rd b/man/write_stan.Rd index 7cc055e8..4080f495 100644 --- a/man/write_stan.Rd +++ b/man/write_stan.Rd @@ -14,7 +14,7 @@ write_stan(object, destination, ...) \item{destination}{(\code{character} or \code{connection})\cr Where to write stan code to.} -\item{...}{Additional arguements} +\item{...}{Additional arguments} } \description{ Write the Stan code for a Stan module. From d4509d7163c7cbe4978c48f7dbbeda3088de8841 Mon Sep 17 00:00:00 2001 From: gowerc Date: Tue, 25 Jun 2024 16:35:06 +0100 Subject: [PATCH 3/8] pedantic lintr --- tests/testthat/test-LongitudinalRandomSlope.R | 1 - 1 file changed, 1 deletion(-) diff --git a/tests/testthat/test-LongitudinalRandomSlope.R b/tests/testthat/test-LongitudinalRandomSlope.R index ca151791..4afbb6ff 100644 --- a/tests/testthat/test-LongitudinalRandomSlope.R +++ b/tests/testthat/test-LongitudinalRandomSlope.R @@ -317,4 +317,3 @@ test_that("Quantity models pass the parser", { ) expect_stan_syntax(stanmod) }) - From 851987ad6728570ae2b2b97f760adfa9004f43d4 Mon Sep 17 00:00:00 2001 From: gowerc Date: Wed, 26 Jun 2024 10:42:05 +0100 Subject: [PATCH 4/8] Fix broken close command --- vignettes/model_fitting.Rmd | 1 - 1 file changed, 1 deletion(-) diff --git a/vignettes/model_fitting.Rmd b/vignettes/model_fitting.Rmd index de50f7dd..c1d8f5df 100644 --- a/vignettes/model_fitting.Rmd +++ b/vignettes/model_fitting.Rmd @@ -240,7 +240,6 @@ tmp <- tempfile() write_stan(simple_model, destination = tmp) first_part <- head(readLines(tmp), 10) cat(paste(first_part, collapse = "\n")) -close(tmp) ``` ## Sampling Parameters From cc1d765d269df7285b1c88e6a820f1b25fffb51d Mon Sep 17 00:00:00 2001 From: Craig Gower-Page Date: Wed, 26 Jun 2024 10:44:04 +0100 Subject: [PATCH 5/8] Update R/JointModel.R Co-authored-by: Isaac Gravestock <83659704+gravesti@users.noreply.github.com> --- R/JointModel.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/JointModel.R b/R/JointModel.R index aa91d0ff..3d2b5e01 100755 --- a/R/JointModel.R +++ b/R/JointModel.R @@ -144,7 +144,7 @@ as.character.JointModel <- function(x, ...) { #' @export write_stan.JointModel <- function(object, destination) { if (is_connection(destination)) { - return(writeLines(as.character(object), fi)) + return(writeLines(as.character(object), con = destination)) } fi <- file(destination, open = "w") writeLines(as.character(object), fi) From 53c6dbae7cc2906f2cac3dd50c2f2d0d133fdad4 Mon Sep 17 00:00:00 2001 From: gowerc Date: Wed, 26 Jun 2024 10:44:44 +0100 Subject: [PATCH 6/8] use explicit arguments --- R/JointModel.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/JointModel.R b/R/JointModel.R index 3d2b5e01..c553b25c 100755 --- a/R/JointModel.R +++ b/R/JointModel.R @@ -147,7 +147,7 @@ write_stan.JointModel <- function(object, destination) { return(writeLines(as.character(object), con = destination)) } fi <- file(destination, open = "w") - writeLines(as.character(object), fi) + writeLines(as.character(object), con = fi) close(fi) } From e95396efd3c967c3dde6c1ea693ebda5fb849a94 Mon Sep 17 00:00:00 2001 From: gowerc Date: Wed, 26 Jun 2024 10:46:33 +0100 Subject: [PATCH 7/8] add test for claretBruno --- tests/testthat/test-LongitudinalClaretBruno.R | 22 +++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/tests/testthat/test-LongitudinalClaretBruno.R b/tests/testthat/test-LongitudinalClaretBruno.R index 4e059c8b..de2b35b8 100644 --- a/tests/testthat/test-LongitudinalClaretBruno.R +++ b/tests/testthat/test-LongitudinalClaretBruno.R @@ -273,3 +273,25 @@ test_that("Can recover known distributional parameters from a SF joint model", { expect_true(all(dat$q99 >= true_values)) expect_true(all(dat$ess_bulk > 100)) }) + + +test_that("Quantity models pass the parser", { + mock_samples <- .JointModelSamples( + model = JointModel(longitudinal = LongitudinalClaretBruno()), + data = structure(1, class = "DataJoint"), + results = structure(1, class = "CmdStanMCMC") + ) + stanmod <- as.StanModule( + mock_samples, + generator = QuantityGeneratorPopulation(1, "A", "B"), + type = "longitudinal" + ) + expect_stan_syntax(stanmod) + + stanmod <- as.StanModule( + mock_samples, + generator = QuantityGeneratorSubject(1, "A"), + type = "longitudinal" + ) + expect_stan_syntax(stanmod) +}) From 627b514ec361f2a09bbc43009946ac8da7b0c1fb Mon Sep 17 00:00:00 2001 From: gowerc Date: Wed, 26 Jun 2024 11:14:58 +0100 Subject: [PATCH 8/8] Fix R CMD check warnings --- R/JointModel.R | 2 +- R/SimLongitudinalClaretBruno.R | 2 +- man/clbr_sld.Rd | 4 ++-- man/write_stan.Rd | 2 +- 4 files changed, 5 insertions(+), 5 deletions(-) diff --git a/R/JointModel.R b/R/JointModel.R index c553b25c..19aa1115 100755 --- a/R/JointModel.R +++ b/R/JointModel.R @@ -142,7 +142,7 @@ as.character.JointModel <- function(x, ...) { #' @rdname write_stan #' @export -write_stan.JointModel <- function(object, destination) { +write_stan.JointModel <- function(object, destination, ...) { if (is_connection(destination)) { return(writeLines(as.character(object), con = destination)) } diff --git a/R/SimLongitudinalClaretBruno.R b/R/SimLongitudinalClaretBruno.R index 5ca97a1d..7c949c34 100644 --- a/R/SimLongitudinalClaretBruno.R +++ b/R/SimLongitudinalClaretBruno.R @@ -175,7 +175,7 @@ sampleSubjects.SimLongitudinalClaretBruno <- function(object, subjects_df) { #' Claret-Bruno Functionals #' -#' @param time (`numeric`)\cr time grid. +#' @param t (`numeric`)\cr time grid. #' @param b (`number`)\cr baseline sld. #' @param g (`number`)\cr growth rate. #' @param c (`number`)\cr resistance rate. diff --git a/man/clbr_sld.Rd b/man/clbr_sld.Rd index e194cf30..fea332cd 100644 --- a/man/clbr_sld.Rd +++ b/man/clbr_sld.Rd @@ -13,6 +13,8 @@ clbr_ttg(t, b, g, c, p) clbr_dsld(t, b, g, c, p) } \arguments{ +\item{t}{(\code{numeric})\cr time grid.} + \item{b}{(\code{number})\cr baseline sld.} \item{g}{(\code{number})\cr growth rate.} @@ -20,8 +22,6 @@ clbr_dsld(t, b, g, c, p) \item{c}{(\code{number})\cr resistance rate.} \item{p}{(\code{number})\cr growth inhibition.} - -\item{time}{(\code{numeric})\cr time grid.} } \value{ The function results. diff --git a/man/write_stan.Rd b/man/write_stan.Rd index 4080f495..86d2933d 100644 --- a/man/write_stan.Rd +++ b/man/write_stan.Rd @@ -7,7 +7,7 @@ \usage{ write_stan(object, destination, ...) -\method{write_stan}{JointModel}(object, destination) +\method{write_stan}{JointModel}(object, destination, ...) } \arguments{ \item{object}{the module.}