Merge pull request #21 from calico/issue-19

Issue 19 - add `n_iterations` as a parameter in `estimate_timecourse_params_tf()`

.github/workflows/R-CMD-check.yaml

@@ -27,9 +27,9 @@ jobs:
       GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }}
 
     steps:
-      - uses: actions/checkout@v2
+      - uses: actions/checkout@v4
 
-      - uses: r-lib/actions/setup-r@v1
+      - uses: r-lib/actions/setup-r@v2
         with:
           r-version: ${{ matrix.config.r }}
 

R/data.R

@@ -0,0 +1,11 @@
+#' Example Timecourse
+#'
+#' An example timecourse of fluorescent reporter dynamics in yeast.
+#'
+#' @format A tibble with 24 rows and 3 columns:
+#' \describe{
+#'   \item{tc_id}{A number representing a distinct timecourse}
+#'   \item{time}{Time (in divisions)}
+#'   \item{abundance}{Relative expression of the reporter}
+#'   }
+"example_timecourse"

R/impulse_fitting.R

@@ -28,6 +28,7 @@
 #'   }
 #' @param fit_intercept If TRUE, the intercept will be fit, if FALSE, the intercept will be constrainted to zero
 #' @param learning_rate learning rate for the Adams optimizer
+#' @param n_iterations the number of iterations to run the optimizer
 #' @param verbose if TRUE then print additional information
 #'
 #' @return a timecourse list:
@@ -77,6 +78,7 @@ estimate_timecourse_params_tf <- function(
     "time_scale" = 15),
   fit_intercept = FALSE,
   learning_rate = 0.1,
+  n_iterations = 100,
   verbose = FALSE
   ) {
 
@@ -105,6 +107,7 @@ estimate_timecourse_params_tf <- function(
   checkmate::assertLogical(use_prior, len = 1)
   checkmate::assertLogical(fit_intercept, len = 1)
   checkmate::assertNumber(learning_rate, lower = 0)
+  checkmate::assertNumber(n_iterations, lower = 1)
   checkmate::assertLogical(verbose, len = 1)
 
   # test parameters supplied for parameter initialization / priors
@@ -153,6 +156,7 @@ estimate_timecourse_params_tf <- function(
         prior_pars = prior_pars,
         fit_intercept = fit_intercept,
         learning_rate = learning_rate,
+        n_iterations = n_iterations,
         verbose = verbose
       )
 
@@ -197,6 +201,7 @@ estimate_one_timecourse_params_tf <- function (
   prior_pars = NULL,
   fit_intercept = FALSE,
   learning_rate = 0.1,
+  n_iterations = 100,
   verbose = FALSE
   ) {
 
@@ -207,6 +212,7 @@ estimate_one_timecourse_params_tf <- function (
   checkmate::assertNumber(learning_rate, lower = 0)
   checkmate::assertLogical(verbose, len = 1)
   checkmate::assertLogical(fit_intercept, len = 1)
+  checkmate::assertNumber(n_iterations, lower = 1)
 
   tfp <- reticulate::import("tensorflow_probability")
 
@@ -368,6 +374,12 @@ estimate_one_timecourse_params_tf <- function (
       )
     } else {
       # constant noise
+
+      noise_value <- ifelse(use_prior, 0.2, 1)
+      if (verbose) {
+        cli::cli_alert("No noise variable is present so setting noise to {.field {noise_value}}")
+      }
+
       noise_entry <- matrix(
         ifelse(use_prior, 0.2, 1),
         nrow = nrow(one_timecourse),
@@ -407,7 +419,7 @@ estimate_one_timecourse_params_tf <- function (
   continue <- TRUE
   while (continue) {
     # train
-    for (i in 1:100) {
+    for (i in 1:n_iterations) {
       sess$run(
         train,
         feed_dict = timecourse_dict

tests/testthat/test-likelihood.R → tests/testthat/test-impulse_fitting.R

@@ -59,7 +59,7 @@ test_that("prior-free minimizing weighted-MSE", {
     tolerance = 0.0001
   )
 
-  })
+})
 
 test_that("Bayesian impulse minimizes MAP estimate", {
 
@@ -107,9 +107,55 @@ test_that("Bayesian impulse minimizes MAP estimate", {
     object = abund_vs_fit$eval_logLik,
     expected = as.numeric(abund_vs_fit$logLik),
     tolerance = 0.001
-    )
+  )
+})
+
+test_that("Sigmoid Likelihood is correct", {
+
+  fitted_kinetics <- impulse::estimate_timecourse_params_tf(
+    measurements = example_timecourse,
+    prior_pars = c(v_sd = 1.2, rate_shape = 2, rate_scale = 0.25, time_shape = 2, time_scale = 15),
+    model = "sigmoid",
+    use_prior = TRUE,
+    n_initializations = 300L,
+    learning_rate = 0.2,
+    n_iterations = 1000,
+    fit_intercept = TRUE
+  )
+
+  # select a consensus parameter set
+  consensus_kinetics <- impulse::reduce_best_timecourse_params(fitted_kinetics)
+
+  # calculate logLik
 
-  ### TO DO - compare MAP parameters to prior estimates
+  fitted_timecourses <- fit_timecourse(
+    consensus_kinetics$parameters %>%
+      dplyr::select(variable, value) %>%
+      tidyr::spread(variable, value),
+    timepts = unique(example_timecourse$time),
+    model = "sigmoid"
+  )
+
+  combined_data <- example_timecourse %>%
+    dplyr::left_join(fitted_timecourses, by = "time")
+
+  #ggplot(combined_data, aes(x = time, y = abundance)) +
+  #  geom_point() +
+  #  geom_line(aes(y = fit), color = "red") +
+  #  ggtitle("Example timecourse with sigmoid fit") +
+  #  theme_minimal()
+
+  timecourse_logLik <- combined_data %>%
+    dplyr::mutate(
+      normal_logLik = dnorm(abundance, fit, sd = 0.2, log = TRUE)
+    ) %>%
+    dplyr::summarize(logLik = sum(normal_logLik))
+
+  testthat::expect_equal(
+    object = timecourse_logLik$logLik,
+    expected = consensus_kinetics$loss$logLik[1],
+    tolerance = 0.01
+    )
 
 })