add error for q < q_obs, add check for this in tests, update package …

…version to 0.0.1
statdivlab · Sep 11, 2019 · bed17dc · bed17dc
1 parent f5824a5
commit bed17dc
Show file tree

Hide file tree

Showing 4 changed files with 26 additions and 14 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: paramedic
 Title: Predicting Absolute and Relative Abundance by Modeling Efficiency to Derive Intervals and Concentrations
-Version: 0.0.0.9000
+Version: 0.0.1
 Authors@R:
     c(person(given = "Brian",
            family = "Williamson",

diff --git a/R/run_paramedic.R b/R/run_paramedic.R
@@ -51,6 +51,8 @@ run_paramedic <- function(W, V,
     } else {
       if (colnames(W)[1] != colnames(V)[1]) stop("W and V must have the same name ")
     }
+    ## error if q < q_obs
+    if (q < q_obs) stop("V must have fewer taxa than W (or the same number of taxa).")
 
     ## ---------------------------
     ## pre-processing and warnings

diff --git a/R/run_paramedic_centered.R b/R/run_paramedic_centered.R
@@ -51,6 +51,8 @@ run_paramedic_centered <- function(W, V,
     } else {
         if (colnames(W)[1] != colnames(V)[1]) stop("W and V must have the same name ")
     }
+    ## error if q < q_obs
+    if (q < q_obs) stop("V must have fewer taxa than W (or the same number of taxa).")
 
     ## ---------------------------
     ## pre-processing and warnings

diff --git a/tests/testthat/test-run_paramedic.R b/tests/testthat/test-run_paramedic.R
@@ -37,29 +37,33 @@ test_that("errors and warnings for run_paramedic work", {
   expect_error(paramedic::run_paramedic(W = example_16S_data[1, 1:10], V = example_qPCR_data,
                                   n_iter = 30, n_burnin = 25, n_chains = 1, stan_seed = 4747,
                                   control = list(max_treedepth = 15)))
+  ## expect error if q < q_obs
+  expect_error(paramedic::run_paramedic(W = example_16S_data[1, 1:3], V = example_qPCR_data,
+                                        n_iter = 30, n_burnin = 25, n_chains = 1, stan_seed = 4747,
+                                        control = list(max_treedepth = 15)))
   ## check to make sure that both W and V have as first column the sample IDs
-  wrong_sample_id_name <- example_16S_data %>% 
-    mutate(SampleID = sample_id) %>% 
+  wrong_sample_id_name <- example_16S_data %>%
+    mutate(SampleID = sample_id) %>%
     select(SampleID, names(example_16S_data)[2:dim(example_16S_data)[2]], -sample_id)
   expect_error(paramedic::run_paramedic(W = wrong_sample_id_name[, 1:10], V = example_qPCR_data,
                                         n_iter = 30, n_burnin = 25, n_chains = 1, stan_seed = 4747,
                                         control = list(max_treedepth = 15)))
   ## check what happens if rows are scrambled
-  scrambled_rows <- example_16S_data %>% 
+  scrambled_rows <- example_16S_data %>%
     arrange(desc(sample_id))
   expect_warning(paramedic::run_paramedic(W = scrambled_rows[, 1:10], V = example_qPCR_data,
                                           n_iter = 30, n_burnin = 25, n_chains = 1, stan_seed = 4747,
                                           control = list(max_treedepth = 15)))
   ## check what happens if columns are scrambled
-  scrambled_cols <- example_16S_data %>% 
+  scrambled_cols <- example_16S_data %>%
     select(sample_id, Lactobacillus.iners, Gardnerella.vaginalis, Lactobacillus.crispatus,
            Lactobacillus.jensenii:Lactobacillus.gasseri)
   expect_warning(paramedic::run_paramedic(W = scrambled_cols, V = example_qPCR_data,
                                           n_iter = 30, n_burnin = 25, n_chains = 1, stan_seed = 4747,
                                           control = list(max_treedepth = 15)))
   ## make sure that code with rows and columns scrambled works
-  scrambled_rows_and_cols <- example_16S_data %>% 
-    arrange(desc(sample_id)) %>% 
+  scrambled_rows_and_cols <- example_16S_data %>%
+    arrange(desc(sample_id)) %>%
     select(sample_id, Lactobacillus.iners, Gardnerella.vaginalis, Lactobacillus.crispatus,
            Lactobacillus.jensenii:Lactobacillus.gasseri)
   expect_warning(mod <- paramedic::run_paramedic(W = scrambled_rows_and_cols, V = example_qPCR_data,
@@ -75,34 +79,38 @@ test_that("errors and warnings for run_paramedic_centered work", {
   expect_error(paramedic::run_paramedic_centered(W = example_16S_data[1, 1:10], V = example_qPCR_data,
                                         n_iter = 30, n_burnin = 25, n_chains = 1, stan_seed = 4747,
                                         control = list(max_treedepth = 15)))
+  ## expect error if q < q_obs
+  expect_error(paramedic::run_paramedic_centered(W = example_16S_data[1, 1:3], V = example_qPCR_data,
+                                        n_iter = 30, n_burnin = 25, n_chains = 1, stan_seed = 4747,
+                                        control = list(max_treedepth = 15)))
   ## check to make sure that both W and V have as first column the sample IDs
-  wrong_sample_id_name <- example_16S_data %>% 
-    mutate(SampleID = sample_id) %>% 
+  wrong_sample_id_name <- example_16S_data %>%
+    mutate(SampleID = sample_id) %>%
     select(SampleID, names(example_16S_data)[2:dim(example_16S_data)[2]], -sample_id)
   expect_error(paramedic::run_paramedic_centered(W = wrong_sample_id_name[, 1:10], V = example_qPCR_data,
                                         n_iter = 30, n_burnin = 25, n_chains = 1, stan_seed = 4747,
                                         control = list(max_treedepth = 15)))
   ## check what happens if rows are scrambled
-  scrambled_rows <- example_16S_data %>% 
+  scrambled_rows <- example_16S_data %>%
     arrange(desc(sample_id))
   expect_warning(paramedic::run_paramedic_centered(W = scrambled_rows[, 1:10], V = example_qPCR_data,
                                           n_iter = 30, n_burnin = 25, n_chains = 1, stan_seed = 4747,
                                           control = list(max_treedepth = 15)))
   ## check what happens if columns are scrambled
-  scrambled_cols <- example_16S_data %>% 
+  scrambled_cols <- example_16S_data %>%
     select(sample_id, Lactobacillus.iners, Gardnerella.vaginalis, Lactobacillus.crispatus,
            Lactobacillus.jensenii:Lactobacillus.gasseri)
   expect_warning(paramedic::run_paramedic_centered(W = scrambled_cols, V = example_qPCR_data,
                                           n_iter = 30, n_burnin = 25, n_chains = 1, stan_seed = 4747,
                                           control = list(max_treedepth = 15)))
   ## make sure that code with rows and columns scrambled works
-  scrambled_rows_and_cols <- example_16S_data %>% 
-    arrange(desc(sample_id)) %>% 
+  scrambled_rows_and_cols <- example_16S_data %>%
+    arrange(desc(sample_id)) %>%
     select(sample_id, Lactobacillus.iners, Gardnerella.vaginalis, Lactobacillus.crispatus,
            Lactobacillus.jensenii:Lactobacillus.gasseri)
   expect_warning(mod <- paramedic::run_paramedic_centered(W = scrambled_rows_and_cols, V = example_qPCR_data,
                                                  n_iter = 30, n_burnin = 25, n_chains = 1, stan_seed = 4747,
                                                  control = list(max_treedepth = 15)))
   mod_summ <- rstan::summary(mod, probs = c(0.025, 0.975))$summary
   expect_equal(mean(mod_summ[grepl("mu", rownames(mod_summ)) & grepl(",1]", rownames(mod_summ)), 1]), mean(example_qPCR_data$Gardnerella.vaginalis), tolerance = 100)
-})
+})