diff --git a/DESCRIPTION b/DESCRIPTION index 869c25c5..5d9c7f20 100755 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -57,7 +57,8 @@ Suggests: flexsurv, purrr, vdiffr, - prodlim + prodlim, + loo Config/testthat/edition: 3 Collate: 'DataSubject.R' diff --git a/NAMESPACE b/NAMESPACE index 5d4f2475..d214f054 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -200,6 +200,7 @@ export(brierScore) export(compileStanModel) export(enableLink) export(generateQuantities) +export(getParameters) export(getPredictionNames) export(initialValues) export(linkDSLD) @@ -211,6 +212,7 @@ export(merge) export(prior_beta) export(prior_cauchy) export(prior_gamma) +export(prior_init_only) export(prior_invgamma) export(prior_logistic) export(prior_loglogistic) @@ -267,6 +269,7 @@ importFrom(ggplot2.utils,geom_km) importFrom(glue,as_glue) importFrom(stats,.checkMFClasses) importFrom(stats,acf) +importFrom(stats,as.formula) importFrom(stats,delete.response) importFrom(stats,model.frame) importFrom(stats,model.matrix) diff --git a/R/Prior.R b/R/Prior.R index 9708844a..c08ee946 100755 --- a/R/Prior.R +++ b/R/Prior.R @@ -337,7 +337,7 @@ prior_beta <- function(a, b) { ) } -#' Only Initial Values Specification +#' Initial Values Specification #' #' @param dist (`Prior`)\cr a prior Distribution #' @family Prior @@ -346,6 +346,7 @@ prior_beta <- function(a, b) { #' This is primarily used for hierarchical parameters whose distributions #' are fixed within the model and cannot be altered by the user. #' +#' @export prior_init_only <- function(dist) { Prior( parameters = list(), diff --git a/R/StanModel.R b/R/StanModel.R index c1cb7f6e..01fcdab7 100644 --- a/R/StanModel.R +++ b/R/StanModel.R @@ -69,7 +69,7 @@ as.list.StanModel <- function(x, ...) { #' @rdname getParameters #' @export -getParameters.StanModel <- function(object) object@parameters +getParameters.StanModel <- function(object, ...) object@parameters #' @export diff --git a/R/defaults.R b/R/defaults.R index dda61bdf..37c20fd8 100755 --- a/R/defaults.R +++ b/R/defaults.R @@ -11,7 +11,7 @@ NULL #' @rdname getParameters #' @export -getParameters.default <- function(object) { +getParameters.default <- function(object, ...) { if (missing(object) || is.null(object)) { return(NULL) } diff --git a/R/generics.R b/R/generics.R index 9a8874b8..18a22092 100755 --- a/R/generics.R +++ b/R/generics.R @@ -102,9 +102,10 @@ as.StanModule <- function(object, ...) { #' from a model. #' #' @param object where to obtain the parameters from. +#' @param ... additional options. #' -#' @keywords internal -getParameters <- function(object) { +#' @export +getParameters <- function(object, ...) { UseMethod("getParameters") } @@ -409,6 +410,9 @@ resolvePromise.default <- function(object, ...) { enableLink <- function(object, ...) { UseMethod("enableLink") } +enableLink.default <- function(object, ...) { + object +} @@ -439,6 +443,7 @@ as_formula <- function(x, ...) { UseMethod("as_formula") } +#' @importFrom stats as.formula #' @export as_formula.default <- function(x, ...) { as.formula(x, ...) diff --git a/_pkgdown.yml b/_pkgdown.yml index 08fd1358..7b63965e 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -124,6 +124,7 @@ reference: - StanModel - StanModule - write_stan + - getParameters - title: Convenience Functions contents: diff --git a/man/getParameters.Rd b/man/getParameters.Rd index b96cd9f3..12f616aa 100644 --- a/man/getParameters.Rd +++ b/man/getParameters.Rd @@ -9,18 +9,20 @@ \alias{getParameters.default} \title{\code{getParameters}} \usage{ -getParameters(object) +getParameters(object, ...) -\method{getParameters}{StanModel}(object) +\method{getParameters}{StanModel}(object, ...) \method{getParameters}{LinkComponent}(object, ...) \method{getParameters}{Link}(object, ...) -\method{getParameters}{default}(object) +\method{getParameters}{default}(object, ...) } \arguments{ \item{object}{where to obtain the parameters from.} + +\item{...}{additional options.} } \description{ Extract any modelling parameters as a \code{\link{ParameterList}} object @@ -34,4 +36,3 @@ Other LinkComponent: \code{\link{initialValues}()} } \concept{LinkComponent} -\keyword{internal} diff --git a/man/prior_init_only.Rd b/man/prior_init_only.Rd index a7a27ffe..21a8bf13 100644 --- a/man/prior_init_only.Rd +++ b/man/prior_init_only.Rd @@ -2,7 +2,7 @@ % Please edit documentation in R/Prior.R \name{prior_init_only} \alias{prior_init_only} -\title{Only Initial Values Specification} +\title{Initial Values Specification} \usage{ prior_init_only(dist) } diff --git a/vignettes/custom-model-dsld.stan b/vignettes/custom-model-dsld.stan new file mode 100644 index 00000000..03ea2403 --- /dev/null +++ b/vignettes/custom-model-dsld.stan @@ -0,0 +1,12 @@ +functions { + // Provide definition for the dsld link function + matrix link_dsld_contrib(matrix time, matrix link_function_inputs) { + int nrow = rows(time); + int ncol = cols(time); + // broadcast input vectors to match the size of the time matrix + matrix[nrow, ncol] baseline = rep_matrix(link_function_inputs[,1], ncol); + matrix[nrow, ncol] shrinkage = rep_matrix(link_function_inputs[,2], ncol); + matrix[nrow, ncol] growth = rep_matrix(link_function_inputs[,3], ncol); + return growth - baseline .* shrinkage .* exp(- shrinkage .* time); + } +} diff --git a/vignettes/custom-model-enable-link.stan b/vignettes/custom-model-enable-link.stan new file mode 100644 index 00000000..82b3413a --- /dev/null +++ b/vignettes/custom-model-enable-link.stan @@ -0,0 +1,7 @@ +transformed parameters { + // Define matrix required for link functions + matrix[n_subjects, 3] link_function_inputs; + link_function_inputs[,1] = baseline_idv; + link_function_inputs[,2] = shrinkage_idv; + link_function_inputs[,3] = growth_idv; +} diff --git a/vignettes/custom-model.Rmd b/vignettes/custom-model.Rmd new file mode 100644 index 00000000..3a602689 --- /dev/null +++ b/vignettes/custom-model.Rmd @@ -0,0 +1,354 @@ +--- +title: "Fitting a Custom Longitudinal Model" +package: jmpost +output: + rmarkdown::html_vignette: + toc: true +link-citations: true +vignette: > + %\VignetteIndexEntry{Fitting a Custom Longitudinal Model} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +editor_options: + chunk_output_type: console +--- + + +```{r setup, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>" +) +``` + +## Introduction + +This vignette shows a complete example for how to fit a custom longitudinal model. Note that full +details for the various different interfaces can be found in the "Extending jmpost" vignette. +This example implements the Wang, Sung et al. 2009 mixed exponential decay and linear growth model +along with an exponential survival model. +In particular the following model will be implemented: + +**Longitudinal Model**: +$$ +\begin{align} +y_{i}(t) &\sim N\left(\mu_i(t),\ \sigma^2\right) \\ \\ +\mu_i(t) &= b_i e^{-s_it} + g_i t \\ \\ +b_i &\sim \text{LogNormal}(\mu_b,\ \sigma_b) \\ +s_i &\sim \text{LogNormal}(\mu_s,\ \sigma_s) \\ +g_i &\sim \text{LogNormal}(\mu_g,\ \sigma_g) \\ +\end{align} +$$ + +Where: +* $i$ is the subject index +* $y_{i}(t)$ is the observed tumour size measurement for subject $i$ at time $t$ +* $\mu_i(t)$ is the expected tumour size measurement for subject $i$ at time $t$ +* $b_i$ is the subject baseline tumour size measurement +* $s_i$ is the subject kinetics shrinkage parameter +* $g_i$ is the subject kinetics tumour growth parameter +* $\mu_{\theta}$ is the population mean for parameter $\theta$ +* $\omega_{\theta}$ is the population standard deviation for parameter $\theta$. + +**Survival Model**: +$$ +\log(h_i(t)) = \log(\lambda_0) + X_i \beta + G(t \mid b_i, s_i, g_i) +$$ + +Where: + +- $\lambda_0$ is the baseline hazard rate. This is because for this +example we are using an exponential survival model e.g. $h_0(t) = \lambda_0$ +- $t$ is the event time +- $G(.)$ is a link function that maps the subjects tumour growth parameters to a contribution to +the log-hazard function +- $X_i$ is the subjects covariate design matrix +- $\beta$ is the corresponding coefficients vector to scale the design matrix covariates +contribution to the log-hazard function + +For this example we will just consider the derivative of the growth function as the link function, +e.g. +$$ +G(t \mid b_i, s_i, g_i) = -s_i b_i e^{-s_i t} + g_i +$$ + + +To keep the example simple, a number of features that have been implemented in the package's +internal models will be skipped; you may wish to consider adding these if implementing +this model in a real project. +In particular the following have been omitted from this example: + +- Handling for censored observations (e.g. observations that are below the limit of quantification) +- Separate populations per study / arm +- Non-centred parameterisation for the hierarchical parameters (this parameterisation leads to better +performance if you have small numbers of observations per each subject). +- Handling negative observation time (e.g. observations that are taken before the start of the study) + +For reference the following libraries will be used during this example: + +```{R} +library(jmpost) +library(ggplot2) +library(dplyr) +library(loo) +``` + +## Generating Simulated Data + +In order to be confident that our model is working correctly we will first generate some simulated +data. This will allow us to compare the true parameter values with the estimated parameter values. +This can be done using the `SimJointData` constructor function as follows: + +```{R} +# Define our simulation parameters + object +SimWang <- setClass( + "SimWang", + contains = "SimLongitudinal", + slots = c( + times = "numeric", + mu_b = "numeric", + mu_s = "numeric", + mu_g = "numeric", + omega_b = "numeric", + omega_s = "numeric", + omega_g = "numeric", + sigma = "numeric", + link_dsld = "numeric" + ) +) + +# Method to generate individual subjects parameters from the hierarchical distributions +sampleSubjects.SimWang <- function(object, subjects_df) { + nsub <- nrow(subjects_df) + subjects_df$b <- stats::rlnorm(nsub, log(object@mu_b), object@omega_b) + subjects_df$s <- stats::rlnorm(nsub, log(object@mu_s), object@omega_s) + subjects_df$g <- stats::rlnorm(nsub, log(object@mu_g), object@omega_g) + subjects_df +} + +# Method to generate observations for each indiviudal subject +sampleObservations.SimWang <- function(object, times_df) { + nobs <- nrow(times_df) + calc_mu <- function(time, b, s, g) b * exp(-s * time) + g * time + calc_dsld <- function(time, b, s, g) -s * b * exp(-s * time) + g + + times_df$mu_sld <- calc_mu(times_df$time, times_df$b, times_df$s, times_df$g) + times_df$dsld <- calc_dsld(times_df$time, times_df$b, times_df$s, times_df$g) + times_df$sld <- stats::rnorm(nobs, times_df$mu_sld, object@sigma) + times_df$log_haz_link <- object@link_dsld * times_df$dsld + times_df +} + + +# Generate simulated data +set.seed(1622) +joint_data_sim <- SimJointData( + design = list(SimGroup(80, "Arm-A", "Study-X")), + survival = SimSurvivalExponential( + lambda = (1 / 400) * 365, + time_max = 4, + time_step = 1 / 365, + lambda_censor = 1 / 9000, + beta_cat = c("A" = 0, "B" = -0.1, "C" = 0.5), + beta_cont = 0.3 + ), + longitudinal = SimWang( + times = c(1, 50, 100, 200, 300, 400, 600, + 800, 1000, 1300, 1600) / 365, + mu_b = 60, + mu_s = 2, + mu_g = 10, + omega_b = 0.3, + omega_s = 0.3, + omega_g = 0.3, + sigma = 1.5, + link_dsld = 0.2 + ), + .silent = TRUE +) + +dat_lm <- joint_data_sim@longitudinal +dat_os <- joint_data_sim@survival + + +# Select 6 random subjects to plot +dat_lm_plot <- dat_lm |> + filter(pt %in% sample(dat_os$pt, 6)) + +ggplot(dat_lm_plot, aes(x = time, y = sld, group = pt, color = pt)) + + geom_line() + + geom_point() + + labs(x = "Time (years)", y = "Tumour Size", col = "Subject") + + theme_bw() + + theme(legend.position = "bottom") +``` + + +## Defining the Longitudinal Model + +The longitudinal model can be implemented by extending the `LongitudinalModel` class. This can be +done as follows: + +```{R} +WangModel <- setClass( + "WangModel", + contains = "LongitudinalModel" +) + +longmodel <- WangModel( + LongitudinalModel( + name = "Wang", + stan = StanModule("custom-model.stan"), + parameters = ParameterList( + Parameter(name = "mu_baseline", prior = prior_lognormal(log(60), 1), size = 1), + Parameter(name = "mu_shrinkage", prior = prior_lognormal(log(2), 1), size = 1), + Parameter(name = "mu_growth", prior = prior_lognormal(log(10), 1), size = 1), + Parameter(name = "sigma_baseline", prior = prior_lognormal(0.3, 1), size = 1), + Parameter(name = "sigma_shrinkage", prior = prior_lognormal(0.3, 1), size = 1), + Parameter(name = "sigma_growth", prior = prior_lognormal(0.3, 1), size = 1), + Parameter(name = "sigma", prior = prior_lognormal(1.5, 1), size = 1), + # The following is only required if we want jmpost to generate + # initial values automatically for us + Parameter( + name = "baseline_idv", + prior = prior_init_only(prior_lognormal(log(60), 1)), + size = "n_subjects" + ), + Parameter( + name = "shrinkage_idv", + prior = prior_init_only(prior_lognormal(log(2), 1)), + size = "n_subjects" + ), + Parameter( + name = "growth_idv", + prior = prior_init_only(prior_lognormal(log(10), 1)), + size = "n_subjects" + ) + ) + ) +) +``` + +Please note that the `parameters` argument is used to specify the priors for the model and +that the `name` argument for the `Parameter`'s objects must match the name of the parameter used +within the corresponding Stan code. + +The `StanModule` object contains all of the stan code used to implement the model. For this +particular model the Stan code specified in the `custom-model.stan` file is as follows: +```{R, results='asis', echo=FALSE} +x <- readLines("./custom-model.stan") +cat(c("\```stan", x, "\```"), sep = "\n") +``` + +## Defining the Link Function + +As stated in the introduction, the link function for this model is going to be the derivative +of the growth function. This can be implemented in using the `jmpost` framework as follows: + +```{R} +enableLink.WangModel <- function(object, ...) { + object@stan <- merge( + object@stan, + StanModule("custom-model-enable-link.stan") + ) + object +} + + +link <- LinkComponent( + stan = StanModule("custom-model-dsld.stan"), + prior = prior_normal(0, 1), + key = "link_dsld" +) +``` + +Where the Stan code for the `custom-model-enable-link.stan` file is as follows: +```{R, results='asis', echo=FALSE} +x <- readLines("./custom-model-enable-link.stan") +cat(c("\```stan", x, "\```"), sep = "\n") +``` + +And the Stan code for the `custom-model-dsld.stan` file is as follows: +```{R, results='asis', echo=FALSE} +x <- readLines("./custom-model-dsld.stan") +cat(c("\```stan", x, "\```"), sep = "\n") +``` + +Note that as only one link function has been defined the `enableLink` method is not strictly +necessary and the Stan code it contains could have been implemented directly in the `stan` slot +of the `LinkComponent` object. However, if you have multiple link functions the `enableLink` +method is required to avoid duplicating the implementation of the `link_function_inputs` matrix. + + +## Sampling the Joint model + +With our longitudinal model defined we can now specify and sample from the full joint model. +```{R} +model <- JointModel( + longitudinal = longmodel, + survival = SurvivalExponential(lambda = prior_gamma(1, 1)), + link = link +) + +joint_data <- DataJoint( + subject = DataSubject( + data = dat_os, + subject = "pt", + arm = "arm", + study = "study" + ), + survival = DataSurvival( + data = dat_os, + formula = Surv(time, event) ~ cov_cat + cov_cont + ), + longitudinal = DataLongitudinal( + data = dat_lm, + formula = sld ~ time + ) +) + +model_samples <- sampleStanModel( + model, + data = joint_data, + iter_warmup = 800, + iter_sampling = 1000, + chains = 1, + refresh = 0, + parallel_chains = 1 +) + +vars <- c( + "mu_baseline", "mu_shrinkage", "mu_growth", "sigma", + "link_dsld", "sm_exp_lambda" +) +as.CmdStanMCMC(model_samples)$summary(vars) +``` + + +## Generating Quantities of Interest + +If you look at the Stan code for the above model you will notice that several of the +requirement elements for the generated quantities have been defined thus enabling us to +generate quantities both at a population and individual subject level. + + +This can be done at the patient level via: + +```{R} +selected_patients <- head(dat_os$pt, 4) +long_quantities_idv <- LongitudinalQuantities( + model_samples, + grid = GridFixed(subjects = selected_patients) +) +autoplot(long_quantities_idv) +``` + +Or at the population level via: + +```{R} +long_quantities_pop <- LongitudinalQuantities( + model_samples, + grid = GridPopulation() +) +autoplot(long_quantities_pop) +``` diff --git a/vignettes/custom-model.stan b/vignettes/custom-model.stan new file mode 100644 index 00000000..a71b89e4 --- /dev/null +++ b/vignettes/custom-model.stan @@ -0,0 +1,82 @@ + +functions { + // Expected tumour size value + vector sld(vector tumour_time, vector baseline, vector shrinkage, vector growth) { + vector[rows(tumour_time)] tumour_value; + tumour_value = baseline .* exp(- shrinkage .* tumour_time) + + growth .* tumour_time; + return tumour_value; + } + + // Define required function for enabling generated quantities + vector lm_predict_individual_patient(vector time, matrix long_gq_parameters) { + return sld( + time, + long_gq_parameters[,1], // baseline + long_gq_parameters[,2], // shrinkage + long_gq_parameters[,3] // growth + ); + } +} + +parameters{ + // Declare individual patient parameters + vector[n_subjects] baseline_idv; + vector[n_subjects] shrinkage_idv; + vector[n_subjects] growth_idv; + + // Declare population level parameters + real mu_baseline; + real mu_shrinkage; + real mu_growth; + real sigma_baseline; + real sigma_shrinkage; + real sigma_growth; + + // Declare standard deviation for the overall model error + real sigma; +} + +transformed parameters{ + + // Calculated the fitted Tumour values + vector[n_tumour_all] Ypred = sld( + tumour_time, + baseline_idv[subject_tumour_index], + shrinkage_idv[subject_tumour_index], + growth_idv[subject_tumour_index] + ); + + // Calculate per observation log-likelihood for {loo} integration + // These values are automatically added to the target for you + Ypred_log_lik = vect_normal_log_dens( + tumour_value, + Ypred, + rep_vector(sigma, n_tumour_all) // broadcast sigma to the length of Ypred + ); +} + +model { + // Define the heirarchical relationship between the individual + // and population level parameters + baseline_idv ~ lognormal(log(mu_baseline), sigma_baseline); + shrinkage_idv ~ lognormal(log(mu_shrinkage), sigma_shrinkage); + growth_idv ~ lognormal(log(mu_growth), sigma_growth); +} + + +generated quantities { + // Enable individual subject predictions / quantities e.g. + // `GridFixed()` / `GridObservation()` / `GridGrouped()` / `GridEven + matrix[n_subjects, 3] long_gq_parameters; + long_gq_parameters[, 1] = baseline_idv; + long_gq_parameters[, 2] = shrinkage_idv; + long_gq_parameters[, 3] = growth_idv; + + // Enable Population level predictions / quantities by taking the median of the + // hierarchical distribution e.g. `GridPopulation()` + matrix[gq_n_quant, 3] long_gq_pop_parameters; + long_gq_pop_parameters[, 1] = rep_vector(mu_baseline, gq_n_quant); + long_gq_pop_parameters[, 2] = rep_vector(mu_shrinkage, gq_n_quant); + long_gq_pop_parameters[, 3] = rep_vector(mu_growth, gq_n_quant); +} diff --git a/vignettes/extending-jmpost.Rmd b/vignettes/extending-jmpost.Rmd index a71c623f..c4c51c94 100644 --- a/vignettes/extending-jmpost.Rmd +++ b/vignettes/extending-jmpost.Rmd @@ -41,7 +41,8 @@ distribution only. Survival distributions are implemented as S4 classes that inherit from `SurvivalModel`. The two main components of the `SurvivalModel` object are the `stan` and `parameters` slots. The `parameters` slot is a `ParametersList()` object which is used to specify the prior distributions -of each parameter used by the survival distribution as well as for the design matrix coefficients. +of each parameter used by the survival distribution as well as for the design matrix coefficients +(please see the "Prior Specification" section below). The `stan` slot is a `StanModule()` object which needs to define the following: @@ -103,8 +104,142 @@ SurvivalWeibullPH <- function( ``` +## Custom Longitudinal Models +Similar to the survival model the longitudinal models are implemented as S4 classes that inherit +from the `LongitudinalModel` class. The two main components of the `LongitudinalModel` object are +the `stan` and `parameters` slots which specify the underlying Stan code and prior distributions +for the models parameters respectively (please see the "Prior Specification" section below). + +Unlike the survival distributions, the longitudinal models are a lot more flexible and have less +constraints on how they are implemented. That is there aren't any specific variables or functions +that you need to define. + +That being said there are a several optional features of `jmpost` that do require the use of specific +interfaces if you want to enable them for your model. + +### 1) `loo` integration + +If you want to use the `loo` package to calculate the leave-one-out cross-validation +then you need to populate the `Ypred_log_lik` vector. This vector should contain the log-likelihood +contribution for each individual tumour observation. This vector +is automatically 0-initialised, thus all your code needs to do is populate it. +```stan +transformed parameters { + Ypred_log_lik = vect_normal_log_dens( + tumour_value, + expected_tumour_value, + rep_vector(lm_rs_sigma, n_tumour_obs) + ); +} +``` + +Where: +- `tumour_value`, `n_tumour_obs` are predefined data objects (see the "Longitudinal Data Objects" +section below) +- `expected_tumour_value` is the expected value of the tumour assessment for each observation +- `lm_rs_sigma` is the standard deviation of the tumour assessment +- `vect_normal_log_dens` is a vectorised version of the normal log-likelihood function provided by +`jmpost` (this is opposed to Stan's inbuilt `normal_lpdf` function which returns a single value +of the sum all of the log-likelihoods) + +### 2) Individual Subject Generated Quantity Integration + +In order to calculate the individual subject generated quantities (via +`GridFixed()` / `GridGrouped()` / etc) you need to define a Stan function with the signature: +```stan +vector lm_predict_individual_patient(vector time, matrix long_gq_parameters) +``` + +Where: + +- `time` is a vector of timepoints for which to calculate the generated quantity +- `long_gq_parameters` is a matrix with one row per subject and one column per parameter + +Likewise, the `long_gq_parameters` object also needs to be defined for the model in the generated +quantities block as a matrix with 1 row per subject and 1 column per parameter. This structure is +to allow for models with subject specific parameters, in particular random effects models. +If your model has the same parameters for all subjects, for example a fixed effects model, +then the value should be repeated for each subject. Subject's values should be in the same order +as their factor levels in the `DataJoint` object. The following is an example from the +`LongitudinalRandomSlope` model: +```stan +generated quantities { + matrix[n_subjects, 2] long_gq_parameters; + long_gq_parameters[, 1] = lm_rs_ind_intercept; + long_gq_parameters[, 2] = lm_rs_ind_rnd_slope; +} +``` +Where: + +- `lm_rs_ind_intercept` and `lm_rs_ind_rnd_slope` are the individual subject's random intercept and +random slope parameters respectively. + +Note that the `long_gq_parameters` matrix should be structured +as your `lm_predict_individual_patient()` function would expect it to be for the +`long_gq_parameters` argument. + +### 3) Population Generated Quantity Integration + +A common use case is to calculate the quantities based on the "population" level parameters which is +supported in `jmpost` via the `GridPopulation()` function. What +this means in practice though is often model and parameterisation specific. +For example some models would +take the median of the distribution whilst others might take the mean or set the random effects +offset to 0. As such if you wish for your model to be compatible with the `GridPopulation()` +then you need to declare and populate the `long_gq_pop_parameters` object with the following +signature: +```stan +matrix[gq_n_quant, 2] long_gq_pop_parameters; +``` + +Note that the number of rows is `gq_n_quant`. This number will be set to the unique number of +combinations of the arm and study factors in the `DataJoint` object. To support populating this +object two additional variables are provided for you namely `gq_long_pop_study_index` and +`gq_long_pop_arm_index` which are vectors that contain the corresponding index of the study and arm +variables for each row in the `long_gq_pop_parameters` matrix. The following is an example +from the `LongitudinalRandomSlope` model: +```stan +generated quantities { + matrix[gq_n_quant, 2] long_gq_pop_parameters; + long_gq_pop_parameters[, 1] = to_vector(lm_rs_intercept[gq_long_pop_study_index]); + long_gq_pop_parameters[, 2] = to_vector(lm_rs_slope_mu[gq_long_pop_arm_index]); +} +``` + +Where: + +- `lm_rs_intercept` and `lm_rs_slope_mu` are the model specific group level +intercept and slope parameters respectively. + +Note that the `long_gq_pop_parameters` matrix should be structured +as your `lm_predict_individual_patient()` function would expect it to be for the +`long_gq_parameters` argument. + +## Prior Specification + +When writing your own custom longitudinal or survival model it is important to understand +how the prior definitions are specified. By default `jmpost` will insert the Stan statements +for any prior distributions based on the `parameters` slot of the model object. +The importance of this is that it means you should not define the prior distributions in the +Stan code itself. Note that this does not apply to hierarchical parameters who must have their +distributions specified in the Stan code. For example in the `LongitudinalRandomSlope` model +their is a different random slope for each treatment arm which is specified in the Stan code as: + +```stan +model { + lm_rs_ind_rnd_slope ~ normal( + lm_rs_slope_mu[subject_arm_index], + lm_rs_slope_sigma + ); +} +``` + +There is however no prior specified in the Stan code for `lm_rs_slope_mu` or `lm_rs_slope_sigma` +as these are handled by the `parameters` slot of the model object as mentioned above. +The main reason for using this approach is that `jmpost` implements the priors in such a way +that users can change them without having to re-compile the Stan model. ## Custom Link Functions