[WIP] errors-in-variables model

DESCRIPTION

@@ -15,7 +15,7 @@ Authors@R: c(
 License: file LICENSE
 Encoding: UTF-8
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.2.0
+RoxygenNote: 7.3.2
 Imports: cmdstanr, posterior, PKPDsim, dplyr, purrr, magrittr, rlang, tidyr, 
   ggplot2, crayon
 Suggests: 

R/PKPDsim_to_stan_data.R

@@ -1,6 +1,7 @@
 #' Prepare data object for use in get_mcmc_posterior()
 #'
 #' @inheritParams new_stan_data
+#' 
 #' @keywords internal
 PKPDsim_to_stan_data <- function(
   regimen, 
@@ -12,6 +13,7 @@ PKPDsim_to_stan_data <- function(
   ruv,
   dose_cmt = 1,
   ltbs = FALSE,
+  eiv_sd = NULL,
   verbose = FALSE
 ) {
   ## Convert regimen, covariates, tdm data
@@ -42,6 +44,17 @@ PKPDsim_to_stan_data <- function(
     c("TIME", "EVID", "AMT", "CMT", "SS", "II", "ADDL", "RATE")
   )
   names(out)[which(names(out) %in% lowercase)] <- tolower(lowercase)
+
+  if(!is.null(eiv_sd)) {
+    if(length(eiv_sd) != 2) {
+      stop("Argument `eiv_sd` needs to be a numeric vector of length 2.")
+    }
+    out$time_sd <- dplyr::if_else(
+      out$evid == 0, 
+      eiv_sd[1], 
+      eiv_sd[2]
+    )
+  }
 
   ## Population parameters and IIV
   for(key in names(parameters)) {

R/get_parameter_tables.R

@@ -1,20 +1,23 @@
 #' Get parameter tables from posterior and prior
 #' 
-#' @param post posterior object from `get_mcmc_posterior`
-#' @param long return tables in long format? Default is `FALSE`
+#' @inheritParams plot_params
 #' 
 #' @export
+#' 
 get_parameter_tables <- function(
   post,
+  params = NULL,
   long = FALSE
 ) {
   if(is.null(post$settings) || is.null(post$draws_df)) {
     stop("Provided posterior object does not contain expected info.")
   }
-  params <- setdiff(
-    names(post$settings$init),
-    post$data$fixed # remove fixed parameters
-  )
+  if(is.null(params)) {
+    params <- setdiff(
+      names(post$settings$init),
+      post$data$fixed # remove fixed parameters
+    )
+  }
   prior_params <- paste0("prior_", params)
   suppressWarnings({
     par_table <- post$draws_df %>% 

R/new_stan_data.R

@@ -17,6 +17,8 @@
 #' log-scale (which then becomes an approximate proportional error).
 #' @param ltbs use log-transform-both-sides approach for observations? Default 
 #' is `FALSE`.
+#' @param eiv_sd vector of length 2 with sd for observation and 
+#' dose / other event types respectively (evid = 0, 1+).
 #' @param verbose verbosity
 #' 
 #' @return Named list suitable for passing on to Torsten.
@@ -66,6 +68,7 @@ new_stan_data <- function(
   ruv,
   dose_cmt = 1,
   ltbs = FALSE,
+  eiv_sd = NULL,
   verbose = FALSE
 ) {
 
@@ -79,6 +82,7 @@ new_stan_data <- function(
     ruv = ruv,
     dose_cmt = dose_cmt,
     ltbs = ltbs,
+    eiv_sd = eiv_sd,
     verbose = verbose
   )
 

R/new_stan_model.R

@@ -21,6 +21,8 @@
 #' @param custom_ipred list of custom definitions for individual predictions for 
 #' each observation type. Sometimes this is needed when more flexibility is 
 #' needed than the `scale` parameter provides.
+#' @param errors_in_variables Either `TRUE` or `FALSE`, should EiV code be 
+#' implemented in the model to probabilistically interpret the time variable?
 #' @param verbose verbosity
 #' @inheritParams new_stan_data
 #' 
@@ -65,6 +67,7 @@ new_stan_model <- function(
   obs_cmt = NULL,
   scale = NULL,
   custom_ipred = NULL,
+  errors_in_variables = FALSE,
   verbose = FALSE
 ) {
   solver <- match.arg(solver)
@@ -118,7 +121,8 @@ new_stan_model <- function(
     obs_types = obs_types,
     obs_cmt = obs_cmt,
     scale = scale,
-    custom_ipred
+    custom_ipred,
+    errors_in_variables = errors_in_variables
   )
 
   model_code <- generate_stan_code(

R/parse_model_definitions.R

@@ -14,7 +14,8 @@ parse_model_definitions <- function(
   obs_types,
   obs_cmt,
   scale,
-  custom_ipred
+  custom_ipred,
+  errors_in_variables
 ) {
 
   def <- list()
@@ -55,7 +56,15 @@ parse_model_definitions <- function(
     "real<lower = 0> ruv_prop_", obs_types, ";", cr, "  ",
     "real<lower = 0> ruv_add_", obs_types, ";"
   )
-
+
+  if(errors_in_variables) {
+    def[["parameters_idv"]] <- paste0(
+      paste0(
+        "array[n_t] real x;"), # will use uniform distribution by defaults
+        sep = "\n"
+      )
+  }
+
   ## Input data, hardcoded
   def[["input_data"]] <- c(
     "int<lower = 1> cmt[n_t];",
@@ -67,7 +76,13 @@ parse_model_definitions <- function(
     "real rate[n_t];",
     "real ii[n_t];"
   )
-
+  if(errors_in_variables) {
+    def[["input_data"]] <- c(
+      def[["input_data"]],
+      "real time_sd[n_t];"
+    )
+  }
+
   ## Parse covariates
   if(!is.null(covariate_definitions)) {
     def[["covariate_data"]] <- paste0(
@@ -104,13 +119,14 @@ parse_model_definitions <- function(
   def[["parameter_definitions"]] <- paste0("real<lower=0> ", names(parameters_sampled), ";")
 
   ## Transformed parameters block:
+  time <- ifelse(errors_in_variables, "x", "time")
   solver_args <- switch(
     solver,
-    "pmx_solve_onecpt" = "time, amt, rate, ii, evid, cmt, addl, ss, theta",
-    "pmx_solve_twocpt" = "time, amt, rate, ii, evid, cmt, addl, ss, theta",
-    "pmx_solve_rk45" = "ode, n_cmt, time, amt, rate, ii, evid, cmt, addl, ss, theta, 1e-5, 1e-8, 1e5",
-    "pmx_solve_adams" = "ode, n_cmt, time, amt, rate, ii, evid, cmt, addl, ss, theta, 1e-5, 1e-8, 1e5",
-    "pmx_solve_bdf" = "ode, n_cmt, time, amt, rate, ii, evid, cmt, addl, ss, theta, 1e-5, 1e-8, 1e5"
+    "pmx_solve_onecpt" = paste0(time, ", amt, rate, ii, evid, cmt, addl, ss, theta"),
+    "pmx_solve_twocpt" = paste0(time, ", amt, rate, ii, evid, cmt, addl, ss, theta"),
+    "pmx_solve_rk45" = paste0("ode, n_cmt, ", time, ", amt, rate, ii, evid, cmt, addl, ss, theta, 1e-5, 1e-8, 1e5"),
+    "pmx_solve_adams" = paste0("ode, n_cmt, ", time, ", amt, rate, ii, evid, cmt, addl, ss, theta, 1e-5, 1e-8, 1e5"),
+    "pmx_solve_bdf" = paste0("ode, n_cmt, ", time,", amt, rate, ii, evid, cmt, addl, ss, theta, 1e-5, 1e-8, 1e5")
   )
   def[["solver_call"]] <- paste0("A = ", solver, "(", solver_args, ");")
   if(!is.null(custom_ipred)) { # custom definition of individual predictions
@@ -195,6 +211,12 @@ parse_model_definitions <- function(
     names(parameters_sampled), "), omega_", 
     names(parameters_sampled), ");"
   )
+  if(errors_in_variables) {
+    def[["likelihood_idv"]] <- paste0(
+      "time ~ normal(x, time_sd);",
+      sep = "\n"
+    )
+  }
 
   def[["likelihood_observed_data"]] <- paste0(
     "if(ltbs_", obs_types, ") {", cr, "  ",

R/plot_params.R

@@ -3,15 +3,22 @@
 #' for the `CL` parameter, it will plot them as well.
 #' 
 #' @param post posterior object from `get_mcmc_posterior()`
+#' @param params vector of parameter names to extract from table with 
+#' posterior draws. If `NULL`, will take default parameters from `init` object.
 #' 
 #' @return a ggplot2 object with posterior distributions (histograms) faceted 
 #' by parameter.
 #'  
 #' @export
-plot_params <- function(post) {
+#' 
+plot_params <- function(
+  post,
+  params = NULL
+) {
 
   par_table_long <- get_parameter_tables(
-    post, 
+    post,
+    params = params,
     long = TRUE
   )
   pars <- par_table_long$posterior %>%

inst/models/template.stan

@@ -31,6 +31,8 @@ transformed data{
 parameters{
 
   // parameter_definitions
+
+  // parameters_idv
 
 }
 
@@ -56,6 +58,8 @@ model{
 
   // likelihood_observed_data
 
+  // likelihood_idv
+
 }
 
 generated quantities{