WIP S4 conversion

SciomeLLC · Jan 17, 2025 · 7ba4f4e · 7ba4f4e
1 parent 50f567d
commit 7ba4f4e
Show file tree

Hide file tree

Showing 11 changed files with 2,351 additions and 92 deletions.
diff --git a/R/Build_Priors4.R b/R/Build_Priors4.R
@@ -0,0 +1,195 @@
+setGeneric("parse_prior", function(prior) standardGeneric("parse_prior"))
+setGeneric("create_continuous_prior",
+  function(prior_list, model, distribution, deg = 2) {
+    standardGeneric("create_continuous_prior")
+  }
+)
+setGeneric("create_dichotomous_prior", function(prior, model) {
+  standardGeneric("create_dichotomous_prior")
+})
+
+setMethod("parse_prior", "list", function(prior) {
+  # Build the distribution from the string
+  rV <- list()
+  rV$prior <- prior$prior
+
+  temp_a <- regexpr("[[][a-zA-Z]+-*[a-zA-Z]+[]]", prior$model)
+  start <- temp_a[1] + 1
+  end <- start + attr(temp_a, "match.length") - 3
+  if (temp_a == -1) {
+    stop("Could not find a distribution for analysis.")
+  }
+  rV$distribution <- substr(prior$model, start, end)
+  rV$model <- prior$mean
+  BMD_Bayes_continuous_model(
+    prior = rV$prior,
+    model = rV$model,
+    distribution = rV$distribution,
+    parameters = character(),
+    mean = rV$model
+  )
+})
+
+setMethod("create_continuous_prior",
+  signature(
+    prior_list = "ANY",
+    model = "character",
+    distribution = "character"
+  ),
+  function(prior_list, model, distribution, deg = 2) {
+    if (!("BMDmodelprior" %in% class(prior_list))) {
+      stop("Prior is not of a 'BMDmodelprior' class. 
+           Define the prior with `create_prior_list()` or similar.")
+    }
+
+    p <- NA
+    if (grepl("aerts", model)) {
+      if (model %in% c("exp-aerts", "invexp-aerts", "hill-aerts", "lognormal-aerts")) {
+        p <- .check_4param(prior_list, distribution)
+      } else if (model %in% c("logistic-aerts", "probit-aerts")) {
+        p <- .check_4param_sigmoidal(prior_list, distribution)
+      } else if (model %in% c("gamma-aerts", "invgamma-aerts", "lomax-aerts",
+                              "invlomax-aerts", "logskew-aerts", "invlogskew-aerts")) {
+        p <- .check_5param(prior_list, distribution)
+      }
+      if (model %in% c("gamma-aerts", "invgamma-aerts")) {
+        p <- .check_5param_gamma(prior_list, distribution)
+      }
+      p$mean <- model
+    } else {
+      # handle the non-aerts
+      if (model == "LMS") {
+        p <- .check_4param(prior_list, distribution)
+        p$mean <- model
+      } else if (model == "gamma-efsa") {
+        p <- .check_4param_gamma(prior_list, distribution)
+        p$mean <- model
+      } else if (model == "hill") {
+        p <- .check_hill(prior_list, distribution)
+      } else if (model == "FUNL") {
+        p <- .check_FUNL(prior_list, distribution)
+      } else if (model == "exp-5") {
+        p <- .check_exp5(prior_list, distribution)
+      } else if (model == "exp-3") {
+        p <- .check_exp3(prior_list, distribution)
+      } else if (model == "polynomial") {
+        p <- .check_polynomial(prior_list, distribution)
+      } else if (model == "power") {
+        p <- .check_power(prior_list, distribution)
+      }
+    }
+
+    ret_obj <- BMD_Bayes_continuous_model(
+      prior = p$prior,
+      distribution = distribution,
+      model = p$model,
+      parameters = p$parameters,
+      mean = p$mean,
+      degree = ifelse(!is.null(p$degree), p$degree, deg)
+    )
+    return(ret_obj)
+  }
+)
+
+ma_build_priors <- function(Y = NA, model_list = NA, distribution_list = NA, ma_type = "EFSA") {
+  if (is.na(model_list[1])) {
+    if (ma_type == "EFSA") {
+      model_list <- c(
+        rep("exp-aerts", 2), rep("invexp-aerts", 2), rep("hill-aerts", 2),
+        rep("lognormal-aerts", 2), rep("gamma-efsa", 2), rep("LMS", 2),
+        rep("probit-aerts", 2), rep("logistic-aerts", 2)
+      )
+      distribution_list <- rep(c("normal", "lognormal"), 8)
+    }
+    else if (ma_type == "original") {
+      model_list <- c(
+        rep("hill", 2), rep("exp-3", 3),
+        rep("exp-5", 3), rep("power", 2)
+      )
+      distribution_list <- c(
+        "normal", "normal-ncv",
+        rep(c("normal", "normal-ncv", "lognormal"), 2),
+        "normal", "normal-ncv"
+      )
+    }
+    else if (ma_type == "all") {
+      model_list <- rep(.continuous_models, each = length(.continuous_distributions))
+      distribution_list <- rep(.continuous_distributions, length(.continuous_models))
+      badd_inds <- c(3, 12, 13:18)
+      model_list <- model_list[-badd_inds]
+      distribution_list <- distribution_list[-badd_inds]
+    }
+  } else {
+    # check that the models are valid
+    if (!all(model_list %in% .continuous_models[-(5:6)])) {
+      stop("Please specify only the following model types: \n
+            \"hill\",\"exp-3\",\"exp-5\",\"power\", \"exp-aerts\", \"invexp-aerts\", \"gamma-aerts\", \"invgamma-aerts\", \"hill-aerts\",
+            \"lomax-aerts\", \"invlomax-aerts\", \"lognormal-aerts\", \"logskew-aerts\", \"invlogskew-aerts\", \"logistic-aerts\", \"probit-aerts\", \"LMS\",
+            \"gamma-efsa\"")
+    }
+    # if no distributions, default to normal-ncv
+    if (is.na(distribution_list[1])) {
+      temp <- unique(model_list)
+      if (length(temp) != length(model_list)) {
+        warning("Removing duplicate models. Please specify distribution_list to avoid this behavior.")
+        model_list <- temp
+      }
+      distribution_list <- rep("normal-ncv", length(temp))
+    } else {
+      # check size compatibility
+      if (length(model_list) != length(distribution_list)) {
+        stop("Please specify a distribution for each model.")
+      }
+      # check distribution validity
+      if (!all(distribution_list %in% .continuous_distributions)) {
+        stop("Please specify only the following distribution types: \"normal\", \"normal-ncv\", \"lognormal\"")
+      }
+    }
+  }
+
+  if (!is.na(Y[1])) {
+    is_neg <- .check_negative_response(Y)
+    if (is_neg) {
+      tmp_idx <- which(distribution_list == "lognormal")
+      model_list <- model_list[-tmp_idx]
+      distribution_list <- distribution_list[-tmp_idx]
+      if (identical(tmp_idx, integer(0))) {
+        warning("Negative response values found. All lognormal models were removed.")
+      }
+    }
+  }
+
+  prior_list <- list()
+  for (ii in seq_along(model_list)) {
+    PR <- .bayesian_prior_continuous_default(model_list[ii], distribution_list[ii], 2)
+
+    if (!is.na(Y[1])) {
+      if (distribution_list[ii] == "lognormal") {
+        if (ncol(Y) > 1) {
+          PR$priors[nrow(PR$priors), 2] <- log(mean(Y[, 3]) ^ 2)
+        } else {
+          PR$priors[nrow(PR$priors), 2] <- log(var(log(Y)))
+        }
+      } else {
+        if (ncol(Y) > 1) {
+          if (distribution_list[ii] == "normal") {
+            PR$priors[nrow(PR$priors), 2] <- log(mean(Y[, 3]) ^ 2)
+          } else {
+            PR$priors[nrow(PR$priors), 2] <- log(abs(mean(Y[1,])) / mean(Y[, 3]) ^ 2)
+          }
+        } else {
+          if (distribution_list[ii] == "normal") {
+            PR$priors[nrow(PR$priors), 2] <- log(var(Y))
+          } else {
+            PR$priors[nrow(PR$priors), 2] <- log(abs(mean(Y)) / var(Y))
+          }
+        }
+      }
+    }
+
+    t_prior_result <- create_continuous_prior(PR, model_list[ii], distribution_list[ii], 2)
+    prior_list[[ii]] <- t_prior_result
+  }
+
+  return(prior_list)
+}
diff --git a/R/NTP.R b/R/NTP.R
@@ -136,6 +136,63 @@ ntp_polyk <- function(dose, tumor, daysOnStudy) {
   cat(sprintf("Poly-3   P-value = %1.4f\n", result[2]))
   cat(sprintf("Poly-6   P-value = %1.4f\n", result[3]))
 }
+
+setClass(
+  "NtpPolyk",
+  slots = c(
+    dose         = "numeric",
+    tumor        = "numeric",
+    daysOnStudy  = "numeric",
+    result       = "matrix"
+  )
+)
+
+ntp_polyk <- function(dose, tumor, daysOnStudy) {
+  # Perform your input checks:
+  if (sum(tumor > 1) > 0) {
+    stop("Tumors need to be 0 or 1.")
+  }
+  if (sum(tumor < 0) > 0) {
+    stop("Tumors need to be 0 or 1.")
+  }
+  if (sum(daysOnStudy < 0) > 0) {
+    stop("Cannot have negative days on study.")
+  }
+  if ((length(dose) != length(tumor)) ||
+      (length(tumor) != length(daysOnStudy))) {
+    stop("All arrays are not of the same length.")
+  }
+  if ((sum(is.na(dose)) + sum(is.na(tumor)) + sum(is.na(daysOnStudy))) > 0) {
+    stop("There is an NA in the data.")
+  }
+
+  # Call your C++ function that does the Poly-K calculation:
+  result <- .polykCPP(dose, tumor, daysOnStudy)
+
+  # Convert to a matrix, set row names (as you did before):
+  result <- as.matrix(result)
+  row.names(result) <- c("Poly 1.5", "Poly-3", "Poly-6")
+
+  # Return a new S4 object
+  new("NtpPolyk",
+      result      = result
+      # If you want to store inputs as well:
+      # dose        = dose,
+      # tumor       = tumor,
+      # daysOnStudy = daysOnStudy
+  )
+}
+
+setMethod(
+  f = "show",
+  signature = "NtpPolyk",
+  definition = function(object) {
+    cat("The results of the Poly-K test for trend.\n")
+    cat(sprintf("Poly-1.5 P-value = %1.4f\n", object@result[1]))
+    cat(sprintf("Poly-3   P-value = %1.4f\n", object@result[2]))
+    cat(sprintf("Poly-6   P-value = %1.4f\n", object@result[3]))
+  }
+)
 ## -----------------------------------------------------------
 ## JONCKHEERE'S TEST
 ## ----------------Changelog----------------------------------

diff --git a/R/RcppExports.bak b/R/RcppExports.bak
@@ -0,0 +1,47 @@
+# Generated by using Rcpp::compileAttributes() -> do not edit by hand
+# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
+
+.owenst_fn <- function(x, fx) {
+    .Call(`_ToxicR_owenst_fn`, x, fx)
+}
+
+.run_single_dichotomous <- function(model, data, pr, options1, options2) {
+    .Call(`_ToxicR_run_single_dichotomous`, model, data, pr, options1, options2)
+}
+
+.run_continuous_single <- function(model, Y, X, prior, options, dist_type) {
+    .Call(`_ToxicR_run_continuous_single`, model, Y, X, prior, options, dist_type)
+}
+
+.run_continuous_ma_laplace <- function(model_priors, model_type, dist_type, Y, X, options) {
+    .Call(`_ToxicR_run_continuous_ma_laplace`, model_priors, model_type, dist_type, Y, X, options)
+}
+
+.run_continuous_ma_mcmc <- function(model_priors, model_type, dist_type, Y, X, options) {
+    .Call(`_ToxicR_run_continuous_ma_mcmc`, model_priors, model_type, dist_type, Y, X, options)
+}
+
+.run_ma_dichotomous <- function(data, priors, models, model_p, is_MCMC, options1, options2) {
+    .Call(`_ToxicR_run_ma_dichotomous`, data, priors, models, model_p, is_MCMC, options1, options2)
+}
+
+.run_dichotomous_single_mcmc <- function(model, Y, D, pr, options) {
+    .Call(`_ToxicR_run_dichotomous_single_mcmc`, model, Y, D, pr, options)
+}
+
+.run_continuous_single_mcmc <- function(model, Y, D, priors, options, is_logNormal, suff_stat) {
+    .Call(`_ToxicR_run_continuous_single_mcmc`, model, Y, D, priors, options, is_logNormal, suff_stat)
+}
+
+.polykCPP <- function(dose, tumor, daysOnStudy) {
+    .Call(`_ToxicR_polyk`, dose, tumor, daysOnStudy)
+}
+
+.setseedGSL <- function(s) {
+    invisible(.Call(`_ToxicR_setseedGSL`, s))
+}
+
+.set_threads <- function(num_threads) {
+    invisible(.Call(`_ToxicR_set_threads`, num_threads))
+}
+