updated vignette to include data on analysis with replicates

DESCRIPTION

@@ -33,7 +33,8 @@ Suggests:
     ggplot2,
     ggforce,
     gridExtra,
-    ggrepel
+    ggrepel,
+    patchwork
 LinkingTo: 
     BH (>= 1.66.0),
     Rcpp (>= 0.12.0),

NEWS

@@ -1,5 +1,6 @@
 Changes in version 1.17
-+ model explanded for designs with biological replicates 
++ model expanded for designs with biological replicates 
++ model expanded for paired sample designs (also with replicates)
 
 Changes in version 1.15
 + model for gene usage (GU) analysis

R/dgu.R

@@ -1,20 +1,22 @@
 
 DGU <- function(ud,
-               mcmc_warmup = 500,
-               mcmc_steps = 1500,
-               mcmc_chains = 4,
-               mcmc_cores = 1,
-               hdi_lvl = 0.95,
-               adapt_delta = 0.95,
-               max_treedepth = 12) {
+                mcmc_warmup = 500,
+                mcmc_steps = 1500,
+                mcmc_chains = 4,
+                mcmc_cores = 1,
+                hdi_lvl = 0.95,
+                adapt_delta = 0.95,
+                max_treedepth = 12,
+                paired = FALSE) {
 
   # check inputs
   check_dgu_input(ud = ud,
-                  mcmc_chains = base::as.integer(x = mcmc_chains),
-                  mcmc_cores = base::as.integer(x = mcmc_cores),
-                  mcmc_steps = base::as.integer(x = mcmc_steps),
-                  mcmc_warmup = base::as.integer(x = mcmc_warmup),
-                  hdi_lvl = hdi_lvl)
+                  mcmc_chains = as.integer(x = mcmc_chains),
+                  mcmc_cores = as.integer(x = mcmc_cores),
+                  mcmc_steps = as.integer(x = mcmc_steps),
+                  mcmc_warmup = as.integer(x = mcmc_warmup),
+                  hdi_lvl = hdi_lvl,
+                  paired = paired)
 
   udr <- ud
   ud <- get_usage(u = udr)
@@ -25,19 +27,20 @@ DGU <- function(ud,
 
   # get model
   m <- get_model(has_conditions = ud$has_conditions, 
-                 has_replicates = ud$has_replicates)
+                 has_replicates = ud$has_replicates,
+                 paired = paired)
 
   # fit model
-  glm <- rstan::sampling(object = m$model,
-                         data = ud,
-                         chains = mcmc_chains,
-                         cores = mcmc_cores,
-                         iter = mcmc_steps,
-                         warmup = mcmc_warmup,
-                         algorithm = "NUTS",
-                         control = control_list,
-                         pars = m$pars,
-                         refresh = 50)
+  glm <- sampling(object = m$model,
+                  data = ud,
+                  chains = mcmc_chains,
+                  cores = mcmc_cores,
+                  iter = mcmc_steps,
+                  warmup = mcmc_warmup,
+                  algorithm = "NUTS",
+                  control = control_list,
+                  pars = m$pars,
+                  refresh = 50)
 
   message("Computing summaries ... \n")
   gu <- get_condition_prop(glm = glm, hdi_lvl = hdi_lvl, ud = ud, 
@@ -47,7 +50,7 @@ DGU <- function(ud,
   dgu_prob <- get_dgu_prob(glm = glm, hdi_lvl = hdi_lvl, ud = ud, 
                            model_type = m$model_type)
   theta <- get_sample_prop_gu(glm = glm, hdi_lvl = hdi_lvl, ud = ud)
-
+  
 
   # ppc
   message("Computing posterior predictions ... \n")
@@ -56,11 +59,11 @@ DGU <- function(ud,
     ppc_condition = get_ppc_condition(glm = glm, ud = ud, hdi_lvl = hdi_lvl))
 
   # result pack
-  return (list(dgu = dgu,
-               dgu_prob = dgu_prob,
-               gu = gu,
-               theta = theta,
-               ppc = ppc,
-               ud = ud,
-               fit = glm))
+  return(list(dgu = dgu,
+              dgu_prob = dgu_prob,
+              gu = gu,
+              theta = theta,
+              ppc = ppc,
+              ud = ud,
+              fit = glm))
 }

R/loo.R

@@ -8,15 +8,17 @@ LOO <- function(ud,
                 mcmc_cores = 1,
                 hdi_lvl = 0.95,
                 adapt_delta = 0.95,
-                max_treedepth = 12) {
+                max_treedepth = 12,
+                paired = FALSE) {
 
   # check inputs
   check_dgu_input(ud = ud,
                   mcmc_chains = as.integer(x = mcmc_chains),
                   mcmc_cores = as.integer(x = mcmc_cores),
                   mcmc_steps = as.integer(x = mcmc_steps),
                   mcmc_warmup = as.integer(x = mcmc_warmup),
-                  hdi_lvl = hdi_lvl)
+                  hdi_lvl = hdi_lvl,
+                  paired = paired)
 
   # process data
   udp <- get_usage(u = ud)
@@ -57,7 +59,8 @@ LOO <- function(ud,
                mcmc_cores = mcmc_cores,
                hdi_lvl = hdi_lvl,
                adapt_delta = adapt_delta,
-               max_treedepth = max_treedepth)
+               max_treedepth = max_treedepth,
+               paired = paired)
 
     if(is.data.frame(out$gu)==TRUE) {
       out$gu$loo_id <- rs[r]

R/utils_input.R

@@ -8,14 +8,16 @@ check_dgu_input <- function(ud,
                             mcmc_cores,
                             mcmc_steps,
                             mcmc_warmup,
-                            hdi_lvl) {
+                            hdi_lvl,
+                            paired) {
 
   if(missing(ud) || is.null(ud) ||
      missing(mcmc_chains) || is.null(mcmc_chains) ||
      missing(mcmc_steps) || is.null(mcmc_steps) ||
      missing(mcmc_warmup) || is.null(mcmc_warmup) ||
      missing(mcmc_cores) || is.null(mcmc_cores) ||
-     missing(hdi_lvl) || is.null(hdi_lvl)) {
+     missing(hdi_lvl) || is.null(hdi_lvl) ||
+     missing(paired) || is.null(paired)) {
     stop("arguments must be specified")
   }
 
@@ -25,6 +27,7 @@ check_dgu_input <- function(ud,
   check_mcmc_chains(mcmc_chains = mcmc_chains)
   check_mcmc_cores(mcmc_cores = mcmc_cores)
   check_hdi(hdi_lvl = hdi_lvl)
+  check_paired(paired = paired)
 }
 
 
@@ -168,3 +171,17 @@ check_ud <- function(ud) {
     }
   }
 }
+
+
+
+# Description:
+# paired input check
+check_paired <- function(paired) {
+  if(length(paired) != 1) {
+    stop('paired must be a logical (TRUE or FALSE)')
+  }
+
+  if(is.logical(paired) == FALSE) {
+    stop('paired must be a logical (TRUE or FALSE)')
+  }
+}

R/utils_usage.R

@@ -93,8 +93,26 @@ get_usage <- function(u) {
               has_conditions = max(condition_ids)>1))
 }
 
-get_model <- function(has_replicates, has_conditions, debug = FALSE) {
+
+
+# Description:
+# checks for structural problems in the data, and mismatch with given inputs
+check_ud_content <- function(ud, paired) {
+  has_conditions <- ud$has_conditions
+  has_replicates <- ud$has_replicates
+  ud <- M$ud$proc_ud
+
+  table(ud$sample_id)
+  diag(table(ud$sample_id, ud$individual_id))
+  identical(diag(table(ud$sample_id, ud$individual_id)))
+}
+
+
+
+
+get_model <- function(has_replicates, has_conditions, paired, debug = FALSE) {
   model_type <- ifelse(test = has_conditions, yes = "DGU", no = "GU")
+
   if(model_type == "GU") {
     if(has_replicates) {
       if(debug) {
@@ -136,6 +154,9 @@ get_model <- function(has_replicates, has_conditions, debug = FALSE) {
                 "Yhat_rep", "Yhat_rep_prop", "Yhat_condition_prop", 
                 "log_lik", "dgu", "dgu_prob", "theta")
       model_name <- "DGU_rep"
+      if(paired) {
+        model_name <- "DGU_rep_paired"
+      }
     } 
     else {
       if(debug) {
@@ -149,6 +170,9 @@ get_model <- function(has_replicates, has_conditions, debug = FALSE) {
                 "Yhat_rep", "Yhat_rep_prop", "Yhat_condition_prop", 
                 "log_lik", "dgu", "dgu_prob", "theta")
       model_name <- "DGU"
+      if(paired) {
+        model_name <- "DGU_paired"
+      }
     }
   }
 
@@ -157,5 +181,6 @@ get_model <- function(has_replicates, has_conditions, debug = FALSE) {
               model_type = model_type,
               pars = pars,
               has_replicates = has_replicates,
-              has_conditions = has_conditions))
+              has_conditions = has_conditions,
+              paired = paired))
 }

man/d_zibb_1.Rd

@@ -8,7 +8,7 @@ A small example dataset that has the following features:
 
   \itemize{
     \item 1 conditions
-    \item 5 replicates (samples) per condition
+    \item 5 samples per condition
     \item 15 Ig genes
   }
 This dataset was simulated from zero-inflated beta-binomial (ZIBB) 

man/d_zibb_2.Rd

@@ -8,7 +8,7 @@ A small example dataset that has the following features:
 
   \itemize{
     \item 2 conditions
-    \item 5 replicates (samples) per condition
+    \item 5 samples per condition
     \item 15 Ig genes
   }
 This dataset was simulated from zero-inflated beta-binomial (ZIBB) 

man/d_zibb_3.Rd

@@ -8,7 +8,7 @@ A small example dataset that has the following features:
 
   \itemize{
     \item 3 conditions
-    \item 5 replicates (samples) per condition
+    \item 5 samples per condition
     \item 15 Ig genes
   }
 This dataset was simulated from zero-inflated beta-binomial (ZIBB)