Added TRANSVIR Case Study 2, changed some internals and added sero check functions

Author: dchodge

NAMESPACE

@@ -3,7 +3,11 @@
 export(addAbkineticsModel)
 export(addCopModel)
 export(addObservationalModel)
-export(add_par_df)
+export(addPrior)
+export(check_exposures_times)
+export(check_input_data)
+export(check_sero_no_single_entries)
+export(check_sero_timings)
 export(createSeroJumpModel)
 export(data_known_exposures_ex1)
 export(data_titre_ex)

R/input_checks.R

@@ -26,7 +26,7 @@ check_time <- function(data_titre) {
 }
 
 
-addExposurePrior_checkempirical <- function(exp_prior, data_t) {
+addExposurePrior_checkempirical <- function(exp_prior, T) {
     columns <- c("time", "prob")
     existing_columns <- names(exp_prior)
 
@@ -39,91 +39,90 @@ addExposurePrior_checkempirical <- function(exp_prior, data_t) {
         stop(paste("Error: Column(s)", paste(names(column_types[!column_types %in% c("numeric", "integer")]), collapse = ", "), "are not numeric or integers in the exposure prior dataframe"))
     }
 
-    if (data_t$T != nrow(exp_prior) ) {
+    if (T != nrow(exp_prior) ) {
         warning("Warning: The number of rows in the exposure prior dataframe does not match the number of time points in the study\n")
     }
 }
 
 
-check_inputs <- function(data_sero, data_known, modeldefinition) {
+check_inputs <- function(data_sero, data_known, biomarkers, exposureTypes, exposureFitted, observationalModel, abkineticsModel, exposurePrior, exposurePriorType) {
     # CHECK inputs of modeldefinition are present
-    if(is.null(modeldefinition$biomarkers)) {
-        stop("Please define the `biomarkers` variable in `modeldefinition`")
-    }
-    if(is.null(modeldefinition$exposureTypes)) {
-        stop("Please define the `exposureTypes` variable in `modeldefinition`")
-    }
-    if(is.null(modeldefinition$exposureFitted)) {
-        stop("Please define the `exposureFitted` variable in `modeldefinition`")
-    }
-    if(is.null(modeldefinition$observationalModel)) {
-        stop("Please define the `observationalModel` structure in `modeldefinition`")
-    }
-    if(is.null(modeldefinition$abkineticsModel)) {
-        stop("Please define the `abkineticsModel` structure in `modeldefinition`")
-    }
+   # if(is.null(modeldefinition$biomarkers)) {
+  #      stop("Please define the `biomarkers`")
+  #  }
+  #  if(is.null(modeldefinition$exposureTypes)) {
+  #      stop("Please define the `exposureTypes`")
+  #  }
+  #  if(is.null(modeldefinition$exposureFitted)) {
+  #      stop("Please define the `exposureFitted`")
+  #  }
+  #  if(is.null(modeldefinition$observationalModel)) {
+  #      stop("Please define the `observationalModel``")
+   # }
+   ## if(is.null(modeldefinition$abkineticsModel)) {
+   #     stop("Please define the `abkineticsModel`")
+   # }
 
 
     # CHECK BIOMARKERS ARE WELL DEFINED
     # Check columns of data_sero match model definition 
     data_sero_name <- data_sero %>% names
-    biomarkers_md <- modeldefinition$biomarkers
-    biomarkers_obs <- modeldefinition$observationalModel$model %>% map(~.x$biomarker) %>% unlist
-    biomarkers_abkin <- modeldefinition$abkineticsModel$model %>% map(~.x$biomarker) %>% unlist %>% unique
+    biomarkers_md <- biomarkers
+    biomarkers_obs <- observationalModel$model %>% map(~.x$biomarker) %>% unlist
+    biomarkers_abkin <- abkineticsModel$model %>% map(~.x$biomarker) %>% unlist %>% unique
 
     for(b in biomarkers_md) {
         if(!b %in% data_sero_name) {
-            stop("Biomarker, ", b, ", in `modeldefinition$biomarkers` is not a column of serological data; `",
+            stop("Biomarker, ", b, ", in `biomarkers` is not a column of serological data; `",
                 paste(data_sero_name, collapse = ", "), "`")
         }
     }
      if(!identical(biomarkers_md, biomarkers_obs) ) {
-        stop("Biomarkers in observationalModel (", paste(biomarkers_obs, collapse = ", "), ") do not match biomarkers in `modeldefinition$biomarkers` (", paste(biomarkers_md, collapse = ", "), ")")
+        stop("Biomarkers in observationalModel (", paste(biomarkers_obs, collapse = ", "), ") do not match biomarkers in `biomarkers` (", paste(biomarkers_md, collapse = ", "), ")")
     }
     if(!identical(biomarkers_md, biomarkers_abkin) ) {
-        stop("Biomarkers in abkineticsModel (", paste(biomarkers_abkin, collapse = ", "), ") do not match biomarkers in `modeldefinition$biomarkers` (", paste(biomarkers_md, collapse = ", "), ")")
+        stop("Biomarkers in abkineticsModel (", paste(biomarkers_abkin, collapse = ", "), ") do not match biomarkers in `biomarkers` (", paste(biomarkers_md, collapse = ", "), ")")
     }  
 
-    exposures_md <- modeldefinition$exposureTypes
-    exposures_obs <- modeldefinition$abkineticsModel$model %>% map(~.x$exposureType) %>% unlist %>% unique
+    exposures_md <- exposureTypes
+    exposures_obs <- abkineticsModel$model %>% map(~.x$exposureType) %>% unlist %>% unique
 
     # CHECK EXPSURETYPES ARE WELL DEFINED
     if (!is.null(data_known)) {
         exposure_type_names <- data_known$exposure_type %>% unique
         for(e in exposure_type_names) {
             if(!e %in% exposures_md) {
                 stop("Exposure type, ", e, ", in known exposure data.frame column 'exposure_type' (", paste(exposure_type_names, collapse = ", "),
-                    "is not defined in `modeldefinition$exposureTypes` (", paste(exposures_md, collapse = ", "), ")")
+                    "is not defined in `exposureTypes` (", paste(exposures_md, collapse = ", "), ")")
             }
         }
     }
     if(!identical(exposures_md, exposures_obs) ) {
-        stop("Exposure types in abkineticsModel (", paste(exposures_obs, collapse = ", "), ") do not match exposure types in `modeldefinition$exposureTypes` (", paste(exposures_md, collapse = ", "), ")")
+        stop("Exposure types in abkineticsModel (", paste(exposures_obs, collapse = ", "), ") do not match exposure types in `exposureTypes` (", paste(exposures_md, collapse = ", "), ")")
     }
-    if(is.null(modeldefinition$exposureFitted)) {
-        stop("`modeldefinition$exposureFitted` is NULL, please define a biomarker to fit.")
+    if(is.null(exposureFitted)) {
+        stop("`exposureFitted` is NULL, please define a biomarker to fit.")
     }
 
-    exposure_fitted <- modeldefinition$exposureFitted
+    exposure_fitted <- exposureFitted
     if(!exposure_fitted %in% exposures_md) {
-        stop("The fitted exposure type, ", exposure_fitted, ", is not defined in, `modeldefinition$exposureTypes`: ", paste(exposures_md, collapse = ", "))
+        stop("The fitted exposure type, ", exposure_fitted, ", is not defined in, `exposureTypes`: ", paste(exposures_md, collapse = ", "))
     }
 
-    names_obs <- modeldefinition$observationalModel$model %>% map(~.x$name) %>% unlist
-    names_abkin <- modeldefinition$abkineticsModel$model %>% map(~.x$name) %>% unlist %>% unique
+    names_obs <- observationalModel$model %>% map(~.x$name) %>% unlist
+    names_abkin <-abkineticsModel$model %>% map(~.x$name) %>% unlist %>% unique
     # Read out into console
     cat("There are ", length(biomarkers_md), " measured biomarkers: ", paste(biomarkers_md, collapse = ", "), "\n")
     cat("There are ", length(exposures_md), " exposure types in the study period: ", paste(exposures_md, collapse = ", "), "\n")
-    cat("The fitted exposure type is ", modeldefinition$exposureFitted, "\n")
+    cat("The fitted exposure type is ", exposureFitted, "\n")
 }
 
 
 
-check_priors <- function(modeldefinition) {
+check_priors <- function(observationalModel, abkineticsModel) {
     priors <- bind_rows(
-        modeldefinition$observationalModel$prior,
-        modeldefinition$abkineticsModel$prior,
-        modeldefinition$copModel$prior
+        observationalModel$prior,
+        abkineticsModel$prior,
     )
     if(any(duplicated(priors$par_name))) {
         stop("Priors: ", paste0(priors$par_name[duplicated(priors$par_name)], collapse = ", "), " are duplicated, please assign original names to each prior")
@@ -139,7 +138,7 @@ check_priors <- function(modeldefinition) {
     }
 
     cat("PRIOR DISTRIBUTIONS", "\n")
-    cat("Prior parameters of observationalModel are: ", paste(modeldefinition$observationalModel$prior$par_name, collapse = ", "), "\n")
-    cat("Prior parameters of abkineticsModel are: ", paste(modeldefinition$abkineticsModel$prior$par_name, collapse = ", "), "\n")
+    cat("Prior parameters of observationalModel are: ", paste(observationalModel$prior$par_name, collapse = ", "), "\n")
+    cat("Prior parameters of abkineticsModel are: ", paste(abkineticsModel$prior$par_name, collapse = ", "), "\n")
 
 }

R/postprocess.R

@@ -53,7 +53,6 @@ postprocess_fit <- function(model_fit) {
 
     n_post <- post$mcmc[[1]] %>% nrow
     post_exp_combine <- post$jump 
-    post_inf_combine <- post$inf 
 
     post_infexp <- 1:n_chains %>% map_df(
         function(i) {

R/utils_data.R

@@ -5,8 +5,8 @@
 #' \describe{
 #'   \item{id}{Unique identifier for each biomarker.}
 #'   \item{time}{Description of the model using the \code{makeModel} and \code{addObservationalModel} functions.}
-#'   \item{titre}{Description of the prior distribution using the \code{add_par_df} function.}
-#'   \item{biomarker}{Description of the prior distribution using the \code{add_par_df} function.}
+#'   \item{titre}{Description of the prior distribution using the \code{addPrior} function.}
+#'   \item{biomarker}{Description of the prior distribution using the \code{addPrior} function.}
 #' }
 #' 
 #' @details Add more data about this model here.
@@ -85,12 +85,12 @@ data_known_exposures_ex1 <- data.frame(
 #' \describe{
 #'   \item{names}{Unique identifier for each exposure type.}
 #'   \item{model}{Description of the model using the \code{makeModel} and \code{addAbkineticsModel} functions.}
-#'   \item{prior}{Description of the prior distribution using the \code{add_par_df} function.}
+#'   \item{prior}{Description of the prior distribution using the \code{addPrior} function.}
 #' }
 #' 
 #' @details Add more data about this model here.
 #'
-#' @seealso \code{\link{makeModel}}, \code{\link{addAbkineticsModel}}, \code{\link{add_par_df}},  for related functions.
+#' @seealso \code{\link{makeModel}}, \code{\link{addAbkineticsModel}},  for related functions.
 #' 
 #' @examples
 #' # Example usage. This describes the antibody kinetics for a SARS-CoV-2 delta wave using IgG data. First define the antibody kinetics function:
@@ -115,9 +115,9 @@ data_known_exposures_ex1 <- data.frame(
 #'            addAbkineticsModel("IgG", "delta", TRUE, c("a_d", "b_d", "c_d"), infSerumKinetics)
 #'        ),
 #'    prior = dplyr::bind_rows(
-#'        add_par_df("a_d", -2, 2, "norm",  0, 1), # ab kinetics
-#'        add_par_df("b_d", 0, 1, "norm",  0.3, 0.05), # ab kinetics
-#'        add_par_df("c_d", 0, 4, "unif", 0,  4) # ab kinetics
+#'        addPrior("a_d", -2, 2, "norm",  0, 1), # ab kinetics
+#'        addPrior("b_d", 0, 1, "norm",  0.3, 0.05), # ab kinetics
+#'        addPrior("c_d", 0, 4, "unif", 0,  4) # ab kinetics
 #'    )
 #' )
 #' 
@@ -134,12 +134,12 @@ NULL
 #' \describe{
 #'   \item{names}{Unique identifier for each biomarker.}
 #'   \item{model}{Description of the model using the \code{makeModel} and \code{addObservationalModel} functions.}
-#'   \item{prior}{Description of the prior distribution using the \code{add_par_df} function.}
+#'   \item{prior}{Description of the prior distribution using the \code{addPrior} function.}
 #' }
 #' 
 #' @details Add more data about this model here.
 #'
-#' @seealso \code{\link{makeModel}}, \code{\link{addObservationalModel}}, \code{\link{add_par_df}},  for related functions.
+#' @seealso \code{\link{makeModel}}, \code{\link{addObservationalModel}}, \code{\link{addPrior}},  for related functions.
 #' 
 #' @examples
 #' # Example usage. This describes the observation model for a SARS-CoV-2 delta wave using IgG data. First define the log likelihood function, which is cauchy, with a LOD at a titre value of log10(40):
@@ -157,7 +157,7 @@ NULL
 #' observationalModel <- list(
 #'     names = c("IgG"),
 #'     model = makeModel(addObservationalModel("IgG", c("sigma"), obsFunction)),
-#'     prior = add_par_df("sigma", 0.0001, 4, "unif", 0.0001, 4) # observational model,
+#'     prior = addPrior("sigma", 0.0001, 4, "unif", 0.0001, 4) # observational model,
 #' )
 #' 
 #'