Fix SF Quantities

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)

R/JointModel.R

@@ -142,9 +142,12 @@ as.character.JointModel <- function(x, ...) {
 
 #' @rdname write_stan
 #' @export
-write_stan.JointModel <- function(object, file_path) {
-    fi <- file(file_path, open = "w")
-    writeLines(as.character(object), fi)
+write_stan.JointModel <- function(object, destination, ...) {
+    if (is_connection(destination)) {
+        return(writeLines(as.character(object), con = destination))
+    }
+    fi <- file(destination, open = "w")
+    writeLines(as.character(object), con = fi)
     close(fi)
 }
 

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`
 #'

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.

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 arguments
 #'
 #' @export
-write_stan <- function(object, file_path) {
+write_stan <- function(object, destination, ...) {
     UseMethod("write_stan")
 }
 

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")
+}

_pkgdown.yml

@@ -143,6 +143,7 @@ reference:
       - as.StanModule.JointModel
       - as.StanModule.Link
       - as.StanModule.LinkComponent
+      - as.StanModule.JointModelSamples
       - as.character.Parameter
       - as.character.Prior
       - show-object

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]);

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)
+})

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)
+})

tests/testthat/test-LongitudinalRandomSlope.R

@@ -295,3 +295,25 @@ 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)
+})

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)
+})

vignettes/model_fitting.Rmd

@@ -237,7 +237,7 @@ 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"))
 ```