Release 3.7.0: seedable MCMC runs

.Rbuildignore

@@ -24,3 +24,5 @@ scratch\.R
 ^debug-data$
 ^cpp$
 ^CONTEXT\.md$
+^cran-comments\.md$
+^CRAN-SUBMISSION$

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: moire
 Title: Multiplicity of Infection and Allele Frequency Recovery from Noisy Polyallelic Genetics Data
-Version: 3.6.0
+Version: 3.7.0
 Authors@R: c(
     person(given = "Maxwell",
            family = "Murphy",
@@ -12,8 +12,11 @@ Authors@R: c(
            role = c("aut", "ths"),
            email = "bryan.greenhouse@ucsf.edu",
            comment = c(ORCID = "0000-0003-0287-9111")))
-Description: A Markov Chain Monte Carlo (MCMC) based approach to Bayesian estimation of individual level multiplicity of infection, within host relatedness, and 
-    population allele frequencies from polyallelic genetic data.
+Description: A Markov Chain Monte Carlo (MCMC) based approach to Bayesian
+    estimation of individual level multiplicity of infection, within host
+    relatedness, and population allele frequencies from polyallelic genetic
+    data. Implements the model described in Murphy and Greenhouse (2024)
+    <doi:10.1093/bioinformatics/btae619>.
 License: GPL (>= 3)
 Encoding: UTF-8
 LazyData: true
@@ -31,6 +34,7 @@ Imports:
     dplyr,
     tidyr,
     stats,
+    parallel,
     purrr,
     rlang,
     ggplot2
@@ -41,7 +45,6 @@ RoxygenNote: 7.3.3
 Suggests: 
     knitr,
     rmarkdown,
-    markdown,
     forcats,
     testthat (>= 3.0.0),
     parallelly

NEWS.md

@@ -1,4 +1,4 @@
-# moire (development version)
+# moire 3.7.0
 
 ## New Features
 
@@ -8,6 +8,20 @@
   `$seed` (printed if the run errors), with offsets as `$chain_seeds`. Replay
   needs the same sampler settings; `num_cores` does not affect the draws.
 
+## Other Changes
+
+* Every exported function now has runnable examples and documents its return
+  value, in preparation for CRAN submission.
+* Corrected the swapped `epsilon_pos` / `epsilon_neg` descriptions in
+  `simulate_observed_allele()`.
+* `parallel` is now declared in `Imports`; `markdown` was dropped from
+  `Suggests` as nothing used it.
+* Parallel tempering temperature adaptation now interpolates temperature as a
+  function of cumulative communication barrier with a monotone cubic Hermite
+  scheme (`src/monotone_interpolator.h`), replacing the bundled third-party
+  `spline.h` and its root solving. Adapted ladders differ slightly from
+  earlier versions; if no barrier has been observed yet the update is skipped.
+
 # moire 3.6.0
 
 ## New Features

R/data.R

@@ -2,10 +2,18 @@
 #'
 #' A simulated dataset created using `simulate_data()`
 #'
+#' @format A list as returned by [simulate_data()], with 100 samples at 30
+#'  loci (15 with 5 alleles, 15 with 10) and a mean COI of 5: observed `data`,
+#'  `sample_ids`, `loci`, `is_missing`, and the simulated truth
+#'  (`allele_freqs`, `sample_cois`, `sample_relatedness`, `true_genotypes`).
 "simulated_data"
 
 #' MCMC results from using the packaged simulated data and calling `run_mcmc()`
 #'
+#' @format A list as returned by [run_mcmc()] on [simulated_data]: a single
+#'  chain with 1000 burnin and 1000 sampling iterations, using 80 parallel
+#'  tempering replicas with adaptive temperatures. The `mcmc_demo` vignette
+#'  shows the call.
 "mcmc_results"
 
 #' Genetic and epidemiological data from Namibia
@@ -22,7 +30,7 @@
 #'   \item{locus}{Genetic locus}
 #'   \item{allele}{Allele observed}
 #' }
-#' @source \url{https://doi.org/10.7554/eLife.43510.018}
+#' @source \doi{10.7554/eLife.43510.018}
 "namibia_data"
 
 #' Allele frequencies for different regions

R/mcmc.R

@@ -79,6 +79,41 @@
 #' and whose `$diagnostics` field has Failure-locus counts and classification
 #' (Consistent failure cause vs Hard starting set). Inspect with
 #' `tryCatch(run_mcmc(...), error = function(e) e$diagnostics)`.
+#'
+#' @return List with elements:
+#' \describe{
+#'   \item{chains}{List with one entry per chain. Each holds the sampled
+#'     `coi`, `allele_freqs`, `eps_pos`, `eps_neg`, and `relatedness`, along
+#'     with log-likelihood traces, acceptance rates, and parallel tempering
+#'     diagnostics. Pass the whole result to the `summarize_*()` functions.}
+#'   \item{args}{The arguments the sampler was called with, including `data`.}
+#'   \item{seed}{Integer seed used for the run.}
+#'   \item{chain_seeds}{Integer seed used for each chain.}
+#'   \item{runtime}{Wall-clock time of the run as a `difftime`.}
+#' }
+#'
+#' @examples
+#' sim <- simulate_data(
+#'   mean_coi = 2,
+#'   num_samples = 10,
+#'   epsilon_pos = 0.01,
+#'   epsilon_neg = 0.1,
+#'   locus_freq_alphas = list(rep(1, 4), rep(1, 4), rep(1, 4))
+#' )
+#'
+#' # A short run for illustration. Real analyses should use far more burnin
+#' # and samples (see the defaults) and parallel tempering (pt_chains) to
+#' # improve mixing.
+#' res <- run_mcmc(
+#'   sim,
+#'   is_missing = sim$is_missing,
+#'   burnin = 100,
+#'   samples_per_chain = 100,
+#'   verbose = FALSE,
+#'   seed = 1
+#' )
+#' res$seed
+#' summarize_coi(res)
 run_mcmc <-
   function(data,
            is_missing = FALSE,

R/simulator.R

@@ -7,6 +7,12 @@
 #'
 #' @param n total number of draws
 #' @param alpha vector controlling the concentration of simplex
+#'
+#' @return Numeric matrix with `n` rows and `length(alpha)` columns; each row
+#'  is a draw from the simplex.
+#'
+#' @examples
+#' rdirichlet(3, alpha = c(1, 1, 1))
 rdirichlet <- function(n, alpha) {
   len_alpha <- length(alpha)
   d <- matrix(rgamma(len_alpha * n, alpha), ncol = len_alpha, byrow = TRUE)
@@ -24,6 +30,12 @@ rdirichlet <- function(n, alpha) {
 #'
 #' @param alpha vector parameter controlling the Dirichlet distribution
 #' @param num_loci total number of loci to draw
+#'
+#' @return Numeric matrix with `length(alpha)` rows and `num_loci` columns;
+#'  each column is an allele frequency vector for one locus.
+#'
+#' @examples
+#' simulate_allele_frequencies(alpha = c(1, 1, 1, 1), num_loci = 3)
 simulate_allele_frequencies <- function(alpha, num_loci) {
   dists <- rdirichlet(num_loci, alpha)
   sapply(seq_len(num_loci), function(x) {
@@ -41,6 +53,12 @@ simulate_allele_frequencies <- function(alpha, num_loci) {
 #'
 #' @param num_samples the total number of biological samples to simulate
 #' @param mean_coi mean multiplicity of infection
+#'
+#' @return Integer vector of length `num_samples` with a COI of at least 1
+#'  for each sample.
+#'
+#' @examples
+#' simulate_sample_coi(num_samples = 10, mean_coi = 2)
 simulate_sample_coi <- function(num_samples, mean_coi) {
   qpois(runif(num_samples, dpois(0, mean_coi), 1), mean_coi)
 }
@@ -59,7 +77,16 @@ simulate_sample_coi <- function(num_samples, mean_coi) {
 #'   multinomial distribution
 #' @param internal_relatedness numeric 0-1 indicating the probability for a
 #'   strain's allele to come from an existing lineage within host
-#' @export
+#'
+#' @return List with one element per sample, each a 1-row integer matrix
+#'  counting the number of distinct strains carrying each allele.
+#'
+#' @examples
+#' simulate_sample_genotype(
+#'   sample_cois = c(1, 2, 3),
+#'   locus_allele_dist = c(0.5, 0.3, 0.2),
+#'   internal_relatedness = c(0, 0, 0.5)
+#' )
 simulate_sample_genotype <- function(sample_cois, locus_allele_dist, internal_relatedness) {
   purrr::map2(sample_cois, internal_relatedness, function(coi, r) {
     related_strains <- rbinom(1, coi - 1, r)
@@ -79,9 +106,18 @@ simulate_sample_genotype <- function(sample_cois, locus_allele_dist, internal_re
 #'
 #' @param alleles A numeric vector representing the number of strains
 #'  contributing each allele
-#' @param epsilon_pos expected number of false negatives
-#' @param epsilon_neg expected number of false positives
+#' @param epsilon_pos expected number of false positives
+#' @param epsilon_neg expected number of false negatives
 #' @param missingness probability that the data is missing
+#'
+#' @return Binary numeric vector the same length as `alleles` indicating
+#'  which alleles were observed. All zeros if the observation is missing.
+#'
+#' @examples
+#' simulate_observed_allele(
+#'   alleles = c(2, 0, 1),
+#'   epsilon_pos = 0.01, epsilon_neg = 0.1, missingness = 0
+#' )
 simulate_observed_allele <- function(alleles, epsilon_pos, epsilon_neg, missingness) {
   # scale eps to the number of alleles so that given a fixed COI, there is a fixed
   # number of expected false positives or negatives across loci of varying
@@ -118,6 +154,16 @@ simulate_observed_allele <- function(alleles, epsilon_pos, epsilon_neg, missingn
 #' @param epsilon_pos expected number of false positives
 #' @param epsilon_neg expected number of false negatives
 #' @param missingness probability of data being missing
+#'
+#' @return List the same length as `true_genotypes`, each element a binary
+#'  vector of observed alleles as returned by [simulate_observed_allele()].
+#'
+#' @examples
+#' true_genotypes <- list(c(1, 0, 1), c(0, 2, 0))
+#' simulate_observed_genotype(
+#'   true_genotypes,
+#'   epsilon_pos = 0.01, epsilon_neg = 0.1, missingness = 0
+#' )
 simulate_observed_genotype <- function(true_genotypes,
                                        epsilon_pos,
                                        epsilon_neg,
@@ -149,8 +195,21 @@ simulate_observed_genotype <- function(true_genotypes,
 #' @param internal_relatedness List of internal relatedness values to be used
 #'  instead of simulating
 #' @param missingness probability of data being missing
-#' @return Simulated data that is structured to go into the MCMC sampler
+#' @return List structured for [run_mcmc()], with elements `data`,
+#'  `sample_ids`, `loci`, and `is_missing`, along with the simulated truth
+#'  (`allele_freqs`, `sample_cois`, `sample_relatedness`, `true_genotypes`)
+#'  and the `input` arguments.
 #'
+#' @examples
+#' sim <- simulate_data(
+#'   mean_coi = 2,
+#'   num_samples = 10,
+#'   epsilon_pos = 0.01,
+#'   epsilon_neg = 0.1,
+#'   locus_freq_alphas = list(rep(1, 4), rep(1, 4), rep(1, 4))
+#' )
+#' sim$sample_cois
+#' str(sim$data, max.level = 1)
 simulate_data <- function(mean_coi = NULL,
                           num_samples,
                           epsilon_pos,