new error model (#2)

* update documentation

* check if initialized to NaN

* fix output bug

* fix initialization bug

* change false positive model

* track variances

* fix adaptive sampling bug

* make eps scale with number of alleles

* implement constrained sampling and remove exact sampling

* remove complexity limit

* update simulator to take allele_freq specification

* fix bug in simulator

* parameterize epsilon as draw from beta

* missed semicolon

* fix sampling bounds

* update simulator and documentation

Author: m-murphy

.gitignore

@@ -27,3 +27,7 @@ src/Makevars
 *.tmp
 
 compile_flags.txt
+
+manuscript/
+
+*.rds

NAMESPACE

@@ -13,9 +13,12 @@ export(simulate_allele_frequencies)
 export(simulate_data)
 export(simulate_observed_genotype)
 export(simulate_sample_coi)
+export(simulate_sample_genotype)
 export(summarize_allele_freq_fn)
 export(summarize_allele_freqs)
 export(summarize_coi)
+export(summarize_epsilon_neg)
+export(summarize_epsilon_pos)
 export(summarize_he)
 importFrom(Rcpp,sourceCpp)
 importFrom(magrittr,"%>%")

R/mcmc.R

@@ -17,23 +17,23 @@
 #'  discard as burnin
 #' @param samples Positive Integer. Number of samples to take
 #'  after burnin
-#' @param complexity_limit Limit on the number of alleles possible
-#'  before augmenting with a latent genetic state representation.
 #' @param verbose Logical indicating if progress is printed
-#' @param eps_pos_0 0-1 Numeric. Initial eps_pos value
+#' @param allele_freq_threshold 0-1 Numeric. Lowest allowed value for an
+#'  allele frequency.
+#' @param eps_pos_0 Numeric. Initial eps_pos value
 #' @param eps_pos_var Numeric. Variance used in sampling eps_pos
-#' @param eps_pos_alpha Positive Numeric. Alpha parameter in
-#'  Beta distribution for eps_pos prior
-#' @param eps_pos_beta Positive Numeric. Beta parameter in
-#'  Beta distribution for eps_pos prior
-#' @param eps_neg_0 0-1 Numeric. Initial eps_neg value
+#' @param eps_pos_shape Positive Numeric. Shape parameter in
+#'  Gamma distribution for eps_pos prior
+#' @param eps_pos_scale Positive Numeric. Scale parameter in
+#'  Gamma distribution for eps_pos prior
+#' @param eps_neg_0 Numeric. Initial eps_neg value
 #' @param eps_neg_var Numeric. Variance used in sampling eps_neg
-#' @param eps_neg_alpha Positive Numeric. Alpha parameter in
-#'  Beta distribution for eps_neg prior
-#' @param eps_neg_beta Positive Numeric. Beta parameter in
-#'  Beta distribution for eps_neg prior
-#' @param max_eps_pos 0-1 Numeric. Maximum allowed value for eps_pos
-#' @param max_eps_neg 0-1 Numeric. Maximum allowed value for eps_neg
+#' @param eps_neg_shape Positive Numeric. Shape parameter in
+#'  Gamma distribution for eps_neg prior
+#' @param eps_neg_scale Positive Numeric. Scale parameter in
+#'  Gamma distribution for eps_neg prior
+#' @param max_eps_pos Numeric. Maximum allowed value for eps_pos
+#' @param max_eps_neg Numeric. Maximum allowed value for eps_neg
 #' @param max_coi Positive Numeric. Maximum allowed complexity of infection
 #' @param allele_freq_vars Positive Numeric. Variance used in sampling allele
 #'  frequencies
@@ -47,21 +47,21 @@ run_mcmc <-
            thin = 1,
            burnin = 1e4,
            samples = 1e4,
-           complexity_limit = 5,
            verbose = TRUE,
-           eps_pos_0 = .01,
-           eps_pos_var = .001,
-           eps_pos_alpha = 1,
-           eps_pos_beta = 99,
-           eps_neg_0 = .05,
-           eps_neg_var = .005,
-           eps_neg_alpha = 5,
-           eps_neg_beta = 95,
-           max_eps_pos = .5,
-           max_eps_neg = .5,
+           allele_freq_threshold = 1e-5,
+           eps_pos_0 = .005,
+           eps_pos_var = 1,
+           eps_pos_alpha = .5,
+           eps_pos_beta = 99.5,
+           eps_neg_0 = .005,
+           eps_neg_var = 1,
+           eps_neg_alpha = .5,
+           eps_neg_beta = 99.5,
+           max_eps_pos = 2,
+           max_eps_neg = 2,
            max_coi = 20,
-           allele_freq_vars = .1,
-           adapt_allele_freq_vars = FALSE) {
+           allele_freq_vars = 1,
+           adapt_allele_freq_vars = TRUE) {
     args <- as.list(environment())
 
     ## if is_missing == FALSE, then generate a default FALSE matrix

R/simulator.R

@@ -25,7 +25,7 @@ rdirichlet <- function(n, alpha) {
 #' @param num_loci total number of loci to draw
 simulate_allele_frequencies <- function(alpha, num_loci) {
   dists <- rdirichlet(num_loci, alpha)
-  lapply(seq_len(num_loci), function(x) {
+  sapply(seq_len(num_loci), function(x) {
     dists[x, ]
   })
 }
@@ -44,15 +44,32 @@ simulate_sample_coi <- function(num_samples, mean_coi) {
 
 #' Simulate sample genotype
 #' @details Simulates sampling the genetics at a single locus given an allele
-#'  frequency distribution and a vector of sample COIs
+#'   frequency distribution and a vector of sample COIs
 #'
 #' @param sample_cois Numeric vector indicating the multiplicity of infection
-#'  for each biological sample
+#'   for each biological sample
 #' @param locus_allele_dist Allele frequencies -- simplex parameter of a
-#'  multinomial distribution
-simulate_sample_genotype <- function(sample_cois, locus_allele_dist) {
+#'   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
+simulate_sample_genotype <- function(sample_cois, locus_allele_dist, internal_relatedness) {
   lapply(sample_cois, function(coi) {
-    rmultinom(1, coi, locus_allele_dist)
+    genotypes <- matrix(nrow = coi, ncol = length(locus_allele_dist))
+    for (i in 1:coi) {
+      if (i == 1) {
+        genotypes[i,] = rmultinom(1, 1, locus_allele_dist)
+      } else {
+        sample_internal = as.logical(rbinom(1, 1, internal_relatedness))
+        if (sample_internal) {
+          genotypes[i,] <- genotypes[sample(1:(i - 1), 1),]
+        } else {
+          genotypes[i,] <- rmultinom(1, 1, locus_allele_dist)
+        }
+      }
+    }
+    g <- colSums(genotypes)
+    g
   })
 }
 
@@ -64,16 +81,27 @@ simulate_sample_genotype <- function(sample_cois, locus_allele_dist) {
 #'
 #' @param alleles A numeric vector representing the number of strains
 #'  contributing each allele
-#' @param epsilon_pos false positive rate
-#' @param epsilon_neg false negative rate
+#' @param epsilon_pos expected number of false negatives
+#' @param epsilon_neg expected number of false positives
 simulate_observed_allele <- function(alleles, epsilon_pos, epsilon_neg) {
-  sapply(alleles, function(allele) {
+  positive_indices <- which(as.logical(alleles)) # True Positives
+  negative_indices <- which(!as.logical(alleles)) # True Negatives
+
+  # 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
+  # diversity.
+  eps_pos_prob = epsilon_pos / length(alleles)
+  eps_neg_prob = epsilon_neg / length(alleles)
+
+  alleles <- sapply(alleles, function(allele) {
     if (allele > 0) {
-      rbinom(1, 1, prob = 1 - epsilon_neg)
+      rbinom(1, 1, prob = 1 - eps_neg_prob)
     } else {
-      rbinom(1, 1, epsilon_pos)
+      rbinom(1, 1, prob = eps_pos_prob)
     }
   })
+
+  return(alleles)
 }
 
 #' Simulate observed genotypes
@@ -85,8 +113,8 @@ simulate_observed_allele <- function(alleles, epsilon_pos, epsilon_neg) {
 #'
 #' @param true_genotypes a list of numeric vectors that are input
 #'  to sim_observed_allele
-#' @param epsilon_pos false positive rate
-#' @param epsilon_neg false negative rate
+#' @param epsilon_pos expected number of false positives
+#' @param epsilon_neg expected number of false negatives
 simulate_observed_genotype <- function(true_genotypes,
                                        epsilon_pos,
                                        epsilon_neg) {
@@ -103,28 +131,31 @@ simulate_observed_genotype <- function(true_genotypes,
 #' @param locus_freq_alphas List of alpha vectors to be used to simulate
 #'  from a Dirichlet distribution to generate allele frequencies.
 #' @param num_samples Total number of biological samples to simulate
-#' @param epsilon_pos False positive rate, between 0 and 1
-#' @param epsilon_neg False negative rate, between 0 and 1
+#' @param epsilon_pos False positive rate, expected number of false positives
+#' @param epsilon_neg False negative rate, expected number of false negatives
+#' @param allele_freqs List of allele frequencies to be used instead of
+#'  simulating allele frequencies
 #' @return Simulated data that is structured to go into the MCMC sampler
 #'
 simulate_data <- function(mean_coi,
-                          locus_freq_alphas,
                           num_samples,
                           epsilon_pos,
-                          epsilon_neg) {
-  allele_freq_dists <- c()
-
-  for (alpha in locus_freq_alphas) {
-    allele_freq_dists <- c(
-      allele_freq_dists,
-      simulate_allele_frequencies(alpha, 1)
-    )
+                          epsilon_neg,
+                          locus_freq_alphas = NULL,
+                          allele_freqs = NULL,
+                          internal_relatedness = 0) {
+  if(is.null(allele_freqs)) {
+    allele_freqs <- list()
+    for (i in 1:length(locus_freq_alphas)) {
+      allele_freqs[[i]] <- simulate_allele_frequencies(locus_freq_alphas[[i]], 1)
+    }
   }
 
+
   sample_cois <- simulate_sample_coi(num_samples, mean_coi)
 
-  true_sample_genotypes <- lapply(allele_freq_dists, function(dist) {
-    simulate_sample_genotype(sample_cois, dist)
+  true_sample_genotypes <- lapply(allele_freqs, function(dist) {
+    simulate_sample_genotype(sample_cois, dist, internal_relatedness)
   })
 
   observed_sample_genotypes <- lapply(
@@ -136,13 +167,14 @@ simulate_data <- function(mean_coi,
   list(
     data = observed_sample_genotypes,
     sample_ids = paste0("S", seq.int(1, num_samples)),
-    loci = paste0("L", seq.int(1, length(locus_freq_alphas))),
-    allele_freqs = allele_freq_dists,
+    loci = paste0("L", seq.int(1, length(allele_freqs))),
+    allele_freqs = allele_freqs,
     sample_cois = sample_cois,
     true_genotypes = true_sample_genotypes,
     input = list(
       mean_coi = mean_coi,
       locus_freq_alphas = locus_freq_alphas,
+      allele_freqs = allele_freqs,
       num_samples = num_samples,
       epsilon_pos = epsilon_pos,
       epsilon_neg = epsilon_neg

R/summary.R

@@ -122,6 +122,76 @@ summarize_coi <- function(mcmc_results, lower_quantile = .025,
   return(coi_data)
 }
 
+#' Summarize epsilon_neg
+#'
+#' @details Summarize epsilon negative results from MCMC. Returns
+#'  a dataframe that contains summaries of the posterior
+#'  distribution of epsilon negative for each biological sample.
+#'
+#' @export
+#'
+#' @param mcmc_results Result of calling run_mcmc()
+#' @param lower_quantile The lower quantile of the posterior
+#'  distribution to return
+#' @param upper_quantile The upper quantile of the posterior
+#'  distribution to return
+summarize_epsilon_neg <- function(mcmc_results, lower_quantile = .025, upper_quantile = .975) {
+  epsilon_neg <- mcmc_results$eps_neg
+  post_eps_neg_lower <- sapply(epsilon_neg, function(x) {
+    quantile(x, lower_quantile)
+  })
+  post_eps_neg_med <- sapply(epsilon_neg, function(x) {
+    quantile(x, .5)
+  })
+  post_eps_neg_upper <- sapply(epsilon_neg, function(x) {
+    quantile(x, upper_quantile)
+  })
+  post_eps_neg_mean <- sapply(epsilon_neg, mean)
+
+  eps_neg_data <- data.frame(
+    sample_id = mcmc_results$args$sample_ids,
+    post_eps_neg_lower, post_eps_neg_med, post_eps_neg_upper, post_eps_neg_mean
+  )
+
+  return(eps_neg_data)
+}
+
+
+#' Summarize epsilon_pos
+#'
+#' @details Summarize epsilon positive results from MCMC. Returns
+#'  a dataframe that contains summaries of the posterior
+#'  distribution of epsilon positive for each biological sample.
+#'
+#' @export
+#'
+#' @param mcmc_results Result of calling run_mcmc()
+#' @param lower_quantile The lower quantile of the posterior
+#'  distribution to return
+#' @param upper_quantile The upper quantile of the posterior
+#'  distribution to return
+summarize_epsilon_pos <- function(mcmc_results, lower_quantile = .025, upper_quantile = .975) {
+  epsilon_pos <- mcmc_results$eps_pos
+  post_eps_pos_lower <- sapply(epsilon_pos, function(x) {
+    quantile(x, lower_quantile)
+  })
+  post_eps_pos_med <- sapply(epsilon_pos, function(x) {
+    quantile(x, .5)
+  })
+  post_eps_pos_upper <- sapply(epsilon_pos, function(x) {
+    quantile(x, upper_quantile)
+  })
+  post_eps_pos_mean <- sapply(epsilon_pos, mean)
+
+  eps_pos_data <- data.frame(
+    sample_id = mcmc_results$args$sample_ids,
+    post_eps_pos_lower, post_eps_pos_med, post_eps_pos_upper, post_eps_pos_mean
+  )
+
+  return(eps_pos_data)
+}
+
+
 #' Summarize Function of Allele Frequencies
 #'
 #' @details General function to summarize the posterior distribution of
@@ -226,7 +296,8 @@ summarize_allele_freqs <- function(mcmc_results,
         post_allele_freqs_lower = post_allele_freqs_lower,
         post_allele_freqs_med = post_allele_freqs_med,
         post_allele_freqs_upper = post_allele_freqs_upper,
-        post_allele_freqs_mean = post_allele_freqs_mean
+        post_allele_freqs_mean = post_allele_freqs_mean,
+        num_alleles = num_alleles
       )
     }
   )

data-raw/mcmc_results.R

@@ -1,33 +1,34 @@
-## ----simulate_data
+## ----load_settings
 set.seed(17325)
 
-mean_moi <- 5
+mean_moi <- 4
 num_biological_samples <- 100
-epsilon_pos <- .01
-epsilon_neg <- .03
+epsilon_pos <- .05
+epsilon_neg <- .05
 
 # Generate the number of alleles at each locus
-allele_counts <- c(rep(5, 15), rep(10, 15), rep(25, 15), rep(50, 15))
+allele_counts <- c(rep(3, 15), rep(5, 15), rep(10, 15))
+
 
 # We'll use flat alpha vectors for our draws from the Dirichlet
 locus_freq_alphas <- lapply(allele_counts, function(allele) rep(1, allele))
 
+## ----simulate_data
 simulated_data <- moire::simulate_data(
-  mean_moi, locus_freq_alphas,
+  mean_moi,
   num_biological_samples,
-  epsilon_pos, epsilon_neg
+  epsilon_pos, epsilon_neg,
+  locus_freq_alphas = locus_freq_alphas
 )
 
 
 ## ----run_mcmc
-burnin <- 1e3
+burnin <- 1e4
 num_samples <- 1e3
 
 mcmc_results <- moire::run_mcmc(
   simulated_data$data, simulated_data$sample_ids, simulated_data$loci,
-  verbose = T, burnin = burnin, samples = num_samples, thin = 1,
-  eps_pos_alpha = 1, eps_pos_beta = 99, complexity_limit = 5,
-  eps_neg_alpha = 10, eps_neg_beta = 990, allele_freq_vars = 1,
+  verbose = T, burnin = burnin, samples = num_samples,
   adapt_allele_freq_vars = TRUE
 )