working on small reviewer changes

ropensci · Dec 31, 2024 · 3db3f78 · 3db3f78
1 parent 8e95efc
commit 3db3f78
Show file tree

Hide file tree

Showing 18 changed files with 543 additions and 528 deletions.
diff --git a/.Rhistory b/.Rhistory
diff --git a/R/detectionPlot.R b/R/detectionPlot.R
@@ -124,8 +124,7 @@ detectionPlot <- function(modelfit, mu.min, mu.max, cov.val = NULL,
   out <- detectionCalculate(modelfit, mu, cov.val, probability, qPCR.N)
 
   # convert to long df
-  '%>%' <- magrittr::`%>%`
-  out_long <- as.data.frame(out) %>%
+  out_long <- as.data.frame(out) |>
     tidyr::pivot_longer(cols =! mu, names_to = 'survey_type')
 
   # get full names

diff --git a/R/globals.R b/R/globals.R
@@ -1 +1 @@
-utils::globalVariables(c("L","median","survey_type","inits","value"))
+utils::globalVariables(c("L_ind","median","survey_type","inits","value"))
diff --git a/R/jointModel.R b/R/jointModel.R
@@ -15,7 +15,7 @@
 #' \href{https://ednajoint.netlify.app/}{Package
 #' Vignette}.
 #'
-#' @srrstats {G1.0,G1.1} This software makes available a novel algorithm/model
+#' @srrstats {G1.0,G1.1} This software makes available a algorithm/model
 #'   that was previously published in the literature. The literature reference
 #'   for the joint model is provided here.
 #'
@@ -78,9 +78,10 @@
 #' @srrstats {BS1.3} Description of parameters used in the computational process
 #'   begins here
 #' @param n.chain Number of MCMC chains. Default value is 4.
-#' @param n.iter.burn Number of warm-up MCMC iterations. Default value is 500.
-#' @param n.iter.sample Number of sampling MCMC iterations. Default value is
-#'   2500.
+#' @param n.warmup A positive integer specifying the number of warm-up MCMC
+#'   iterations. Default value is 500.
+#' @param n.iter A positive integer specifying the number of iterations for each
+#'   chain (including warmup). Default value is 3000.
 #' @param thin A positive integer specifying the period for saving samples.
 #'   Default value is 1.
 #' @param adapt_delta Numeric value between 0 and 1 indicating target average
@@ -171,17 +172,17 @@
 #' fit.q <- jointModel(data = greencrabData, cov = NULL, family = "negbin",
 #'                     p10priors = c(1,20), q = TRUE, phipriors = c(0.25,0.25),
 #'                     multicore = FALSE, initial_values = NULL,
-#'                     n.chain = 4, n.iter.burn = 500,
-#'                     n.iter.sample = 2500, thin = 1, adapt_delta = 0.9,
+#'                     n.chain = 4, n.warmup = 500,
+#'                     n.iter = 3000, thin = 1, adapt_delta = 0.9,
 #'                     verbose = TRUE, seed = 123)
 #' }
 #'
 
 jointModel <- function(data, cov = NULL, family = 'poisson',
                        p10priors = c(1,20),
                        q = FALSE, phipriors = NULL, multicore = FALSE,
-                       initial_values = NULL, n.chain = 4, n.iter.burn = 500,
-                       n.iter.sample = 2500, thin = 1, adapt_delta = 0.9,
+                       initial_values = NULL, n.chain = 4, n.warmup = 500,
+                       n.iter = 3000, thin = 1, adapt_delta = 0.9,
                        verbose = TRUE, seed = NULL) {
 
 
@@ -218,8 +219,8 @@ jointModel <- function(data, cov = NULL, family = 'poisson',
   }
 
   # all models
-  all_checks(data,cov,family,p10priors,phipriors,n.chain,n.iter.burn,
-             n.iter.sample,thin,adapt_delta,seed)
+  all_checks(data,cov,family,p10priors,phipriors,n.chain,n.warmup,
+             n.iter,thin,adapt_delta,seed)
 
   # initial value checks
   if(all(!is.null(initial_values))){
@@ -233,20 +234,18 @@ jointModel <- function(data, cov = NULL, family = 'poisson',
   ###
   #convert data to long format
 
-  '%>%' <- magrittr::`%>%`
-
   # convert qPCR data to long format
   #' @srrstats {G2.7} Use as.data.frame() to allow input list of any tabular
   #'   form (i.e., matrix, etc.)
-  qPCR_all <- as.data.frame(data$qPCR.N) %>%
-    dplyr::mutate(L_ind = 1:dim(data$qPCR.N)[1]) %>%
-    tidyr::pivot_longer(cols =! L_ind, values_to = 'n_N') %>%
+  qPCR_all <- as.data.frame(data$qPCR.N) |>
+    dplyr::mutate(L_ind = 1:dim(data$qPCR.N)[1]) |>
+    tidyr::pivot_longer(cols =! L_ind, values_to = 'n_N') |>
     #' @srrstats {G2.15} Software does not assume non-missingness and actually
     #'   expects it if the number of observations across sites is unequal
     tidyr::drop_na()
-  qPCR.K_df <- as.data.frame(data$qPCR.K) %>%
-    dplyr::mutate(L_ind = 1:dim(data$qPCR.K)[1]) %>%
-    tidyr::pivot_longer(cols =! L_ind, values_to = 'n_K') %>%
+  qPCR.K_df <- as.data.frame(data$qPCR.K) |>
+    dplyr::mutate(L_ind = 1:dim(data$qPCR.K)[1]) |>
+    tidyr::pivot_longer(cols =! L_ind, values_to = 'n_K') |>
     #' @srrstats {G2.15} Software does not assume non-missingness and actually
     #'   expects it if the number of observations across sites is unequal
     tidyr::drop_na()
@@ -255,9 +254,9 @@ jointModel <- function(data, cov = NULL, family = 'poisson',
   # convert count data to long format
   #' @srrstats {G2.7} Use as.data.frame() to allow input list of any tabular
   #'   form (i.e., matrix, etc.)
-  count_all <- as.data.frame(data$count) %>%
-    dplyr::mutate(L_ind = 1:dim(data$count)[1]) %>%
-    tidyr::pivot_longer(cols =! L_ind, values_to = 'count') %>%
+  count_all <- as.data.frame(data$count) |>
+    dplyr::mutate(L_ind = 1:dim(data$count)[1]) |>
+    tidyr::pivot_longer(cols =! L_ind, values_to = 'count') |>
     #' @srrstats {G2.15} Software does not assume non-missingness and actually
     #' expects it if the number of observations across sites is unequal
     tidyr::drop_na()
@@ -301,9 +300,9 @@ jointModel <- function(data, cov = NULL, family = 'poisson',
     q_ref <- 1
     #' @srrstats {G2.7} Use as.data.frame() to allow input list of any tabular
     #'   form (i.e., matrix, etc.)
-    count.type_df <- as.data.frame(data$count.type) %>%
-      dplyr::mutate(L_ind = 1:dim(data$count.type)[1]) %>%
-      tidyr::pivot_longer(cols =! L_ind, values_to = 'count.type') %>%
+    count.type_df <- as.data.frame(data$count.type) |>
+      dplyr::mutate(L_ind = 1:dim(data$count.type)[1]) |>
+      tidyr::pivot_longer(cols =! L_ind, values_to = 'count.type') |>
       tidyr::drop_na()
     count_all$count.type <- count.type_df$count.type
 
@@ -361,6 +360,7 @@ jointModel <- function(data, cov = NULL, family = 'poisson',
     nsitecov = length(cov)+1,
     mat_site = as.matrix(site_mat),
     p10priors = c(mu_ln,sigma_ln),
+    alphapriors = c(0,10),
     control = list(adapt_delta = adapt_delta)
   )
 
@@ -393,6 +393,12 @@ jointModel <- function(data, cov = NULL, family = 'poisson',
       phipriors = c(1,1),
       negbin = 0
     )
+  } else if(get_family_index(family)==3){
+    model_data <- rlist::list.append(
+      model_data,
+      bgammapriors = c(0.25,0.25),
+      agammapriors = c(0.01,0.01)
+    )
   }
 
   # set up core number
@@ -428,10 +434,8 @@ jointModel <- function(data, cov = NULL, family = 'poisson',
     #'   integers for sampling arguments
     chains = as.integer(n.chain),
     thin = as.integer(thin),
-    warmup = as.integer(n.iter.burn),
-    iter = (
-      as.integer(n.iter.burn) + as.integer(n.iter.sample)
-    ),
+    warmup = as.integer(n.warmup),
+    iter = as.integer(n.iter),
     init = inits,
     refresh = ifelse(verbose == TRUE,500,0)
   )

diff --git a/R/traditionalModel.R b/R/traditionalModel.R
@@ -54,9 +54,10 @@
 #' @srrstats {BS1.3} Description of parameters used in the computational
 #'   process begins here
 #' @param n.chain Number of MCMC chains. Default value is 4.
-#' @param n.iter.burn Number of warm-up MCMC iterations. Default value is 500.
-#' @param n.iter.sample Number of sampling MCMC iterations. Default value is
-#'   2500.
+#' @param n.warmup A positive integer specifying the number of warm-up MCMC
+#'   iterations. Default value is 500.
+#' @param n.iter A positive integer specifying the number of iterations for each
+#'   chain (including warmup). Default value is 3000.
 #' @param thin A positive integer specifying the period for saving samples.
 #'   Default value is 1.
 #' @param adapt_delta Numeric value between 0 and 1 indicating target average
@@ -120,16 +121,16 @@
 #' fit.q <- traditionalModel(data = greencrabData, family = "negbin", q = TRUE,
 #'                           phipriors = c(0.25,0.25), multicore = FALSE,
 #'                           initial_values = NULL, n.chain = 4,
-#'                           n.iter.burn = 500, n.iter.sample = 2500, thin = 1,
+#'                           n.warmup = 500, n.iter = 3000, thin = 1,
 #'                           adapt_delta = 0.9, verbose = TRUE, seed = 123)
 #' }
 #'
 
 traditionalModel <- function(data, family = 'poisson',
                              q = FALSE, phipriors = NULL,
                              multicore = FALSE, initial_values = NULL,
-                             n.chain = 4, n.iter.burn = 500,
-                             n.iter.sample = 2500, thin = 1,
+                             n.chain = 4, n.warmup = 500,
+                             n.iter = 3000, thin = 1,
                              adapt_delta = 0.9, verbose = TRUE, seed = NULL) {
 
 
@@ -150,7 +151,7 @@ traditionalModel <- function(data, family = 'poisson',
   #' @srrstats {G2.1} Types of inputs are checked/asserted using this helper
   #'   function
   traditionalModel_input_checks(data, family, q, phipriors, n.chain,
-                                n.iter.burn, n.iter.sample,
+                                n.warmup, n.iter,
                                 thin, adapt_delta, seed)
 
   # initial value checks
@@ -163,15 +164,12 @@ traditionalModel <- function(data, family = 'poisson',
   }
 
   ###
-  #convert data to long format
-  '%>%' <- magrittr::`%>%`
-
   #convert count data to long format
   #' @srrstats {G2.7} Use as.data.frame() to allow input list of any tabular
   #'   form (i.e., matrix, etc.)
-  count_all <- as.data.frame(data$count) %>%
-    dplyr::mutate(L_ind = 1:dim(data$count)[1]) %>%
-    tidyr::pivot_longer(cols=!L_ind,values_to = 'count') %>%
+  count_all <- as.data.frame(data$count) |>
+    dplyr::mutate(L_ind = 1:dim(data$count)[1]) |>
+    tidyr::pivot_longer(cols=!L_ind,values_to = 'count') |>
     #' @srrstats {G2.15} Software does not assume non-missingness and actually
     #'   expects it if the number of observations across sites is unequal
     tidyr::drop_na()
@@ -181,9 +179,9 @@ traditionalModel <- function(data, family = 'poisson',
     q_ref <- 1
     #' @srrstats {G2.7} Use as.data.frame() to allow input list of any tabular
     #'   form (i.e., matrix, etc.)
-    count.type_df <- as.data.frame(data$count.type) %>%
-      dplyr::mutate(L_ind = 1:dim(data$count.type)[1]) %>%
-      tidyr::pivot_longer(cols=!L_ind,values_to = 'count.type') %>%
+    count.type_df <- as.data.frame(data$count.type) |>
+      dplyr::mutate(L_ind = 1:dim(data$count.type)[1]) |>
+      tidyr::pivot_longer(cols=!L_ind,values_to = 'count.type') |>
       #' @srrstats {G2.15} Software does not assume non-missingness and
       #'   actually expects it if the number of observations across sites is
       #'   unequal
@@ -244,6 +242,12 @@ traditionalModel <- function(data, family = 'poisson',
       phipriors = c(1,1),
       negbin = 0
     )
+  } else if(family == 'gamma'){
+    model_data <- rlist::list.append(
+      model_data,
+      betapriors = c(0.25,0.25),
+      alphapriors = c(0.01,0.01)
+    )
   }
 
   # get stan model
@@ -273,10 +277,8 @@ traditionalModel <- function(data, family = 'poisson',
     #'   integers for sampling arguments
     chains = as.integer(n.chain),
     thin = as.integer(thin),
-    warmup = as.integer(n.iter.burn),
-    iter = (
-      as.integer(n.iter.burn) + as.integer(n.iter.sample)
-    ),
+    warmup = as.integer(n.warmup),
+    iter = as.integer(n.iter),
     init = inits,
     refresh = ifelse(verbose == TRUE,500,0)
   )

diff --git a/R/utils-joint.R b/R/utils-joint.R
@@ -380,7 +380,7 @@ no_catchability_checks <- function(data,cov){
 #' @noRd
 # input checks for all variations
 all_checks <- function(data, cov, family, p10priors, phipriors, n.chain,
-                       n.iter.burn, n.iter.sample, thin, adapt_delta, seed){
+                       n.warmup, n.iter, thin, adapt_delta, seed){
 
 
   ## make sure count, qPCR.N, and qPCR.K are not zero-length
@@ -584,15 +584,15 @@ all_checks <- function(data, cov, family, p10priors, phipriors, n.chain,
     stop(errMsg)
   }
 
-  ## check length and range of n.iter.sample
-  if(any(length(as.integer(n.iter.sample)) > 1 | n.iter.sample < 1)){
-    errMsg <- "n.iter.sample should be an integer > 0 and of length 1."
+  ## check length and range of n.iter
+  if(any(length(as.integer(n.iter)) > 1 | n.iter < 1)){
+    errMsg <- "n.iter should be an integer > 0 and of length 1."
     stop(errMsg)
   }
 
-  ## check length and range of n.iter.burn
-  if(any(length(as.integer(n.iter.burn)) > 1 | n.iter.burn < 1)){
-    errMsg <- "n.iter.burn should be an integer > 0 and of length 1."
+  ## check length and range of n.warmup
+  if(any(length(as.integer(n.warmup)) > 1 | n.warmup < 1)){
+    errMsg <- "n.warmup should be an integer > 0 and of length 1."
     stop(errMsg)
   }
 

diff --git a/R/utils-traditional.R b/R/utils-traditional.R
@@ -73,7 +73,7 @@ init_trad <- function(n.chain, count_all, initial_values){
 #'   messages
 #' @noRd
 traditionalModel_input_checks <- function(data, family, q, phipriors, n.chain,
-                                          n.iter.burn, n.iter.sample,
+                                          n.warmup, n.iter,
                                           thin, adapt_delta, seed){
 
   ## make sure all data tags are valid -- if q == TRUE
@@ -233,15 +233,15 @@ traditionalModel_input_checks <- function(data, family, q, phipriors, n.chain,
     stop(errMsg)
   }
 
-  ## check length and range of n.iter.sample
-  if(any(length(as.integer(n.iter.sample)) > 1 | n.iter.sample < 1)){
-    errMsg <- "n.iter.sample should be an integer > 0 and of length 1."
+  ## check length and range of n.iter
+  if(any(length(as.integer(n.iter)) > 1 | n.iter < 1)){
+    errMsg <- "n.iter should be an integer > 0 and of length 1."
     stop(errMsg)
   }
 
-  ## check length and range of n.iter.burn
-  if(any(length(as.integer(n.iter.burn)) > 1 | n.iter.burn < 1)){
-    errMsg <- "n.iter.burn should be an integer > 0 and of length 1."
+  ## check length and range of n.warmup
+  if(any(length(as.integer(n.warmup)) > 1 | n.warmup < 1)){
+    errMsg <- "n.warmup should be an integer > 0 and of length 1."
     stop(errMsg)
   }
 

diff --git a/eDNAjoint.Rproj b/eDNAjoint.Rproj
@@ -1,4 +1,5 @@
 Version: 1.0
+ProjectId: 74641fc2-7fcd-46e7-bf18-ec0a939e6539
 
 RestoreWorkspace: No
 SaveWorkspace: No

diff --git a/inst/stan/joint_continuous.stan b/inst/stan/joint_continuous.stan
@@ -22,6 +22,9 @@ data{/////////////////////////////////////////////////////////////////////
     array[S_dna] int<lower=1> N_dna;   // number of qPCR replicates per site of unpaired samples
     array[S_dna] int<lower=0> K_dna; // number of qPCR detections among these replicates of unpaired samples
     array[2] real p10priors; // priors for normal distrib on p10
+    array[2] real alphapriors; // priors for site covariate shrinkage params
+    array[2] real bgammapriors; // priors for beta param in expected catch rate
+    array[2] real agammapriors; // priors for alpha param in expected catch rate
     int<lower=0> nparams;  // number of gear types
     array[n_C] int<lower=1> mat;  // vector of gear type integers
     int<lower=0> nsitecov;  // number of site-level covariates
@@ -81,9 +84,9 @@ model{/////////////////////////////////////////////////////////////////////
 
   //priors
   log_p10 ~ normal(p10priors[1], p10priors[2]); // p10 prior
-  alpha ~ normal(0,10); // sitecov shrinkage priors
-  beta_gamma ~ gamma(0.25,0.25);
-  alpha_gamma ~ gamma(0.01,0.01);
+  alpha ~ normal(alphapriors[1], alphapriors[2]); // sitecov shrinkage priors
+  beta_gamma ~ gamma(bgammapriors[1], bgammapriors[2]);
+  alpha_gamma ~ gamma(agammapriors[1], agammapriors[2]);
 
 }
 

diff --git a/inst/stan/joint_count.stan b/inst/stan/joint_count.stan
@@ -21,6 +21,7 @@ data{/////////////////////////////////////////////////////////////////////
     array[S_dna] int<lower=1> N_dna;   // number of qPCR replicates per site of unpaired samples
     array[S_dna] int<lower=0> K_dna; // number of qPCR detections among these replicates of unpaired samples
     array[2] real p10priors; // priors for normal distrib on p10
+    array[2] real alphapriors; // priors for site covariate shrinkage params
     int<lower=0> nparams;  // number of gear types
     array[n_C] int<lower=1> mat;  // vector of gear type integers
     int<lower=0> nsitecov;  // number of site-level covariates
@@ -82,7 +83,7 @@ model{/////////////////////////////////////////////////////////////////////
 
   //priors
   log_p10 ~ normal(p10priors[1], p10priors[2]); // p10 prior
-  alpha ~ normal(0,10); // sitecov shrinkage priors
+  alpha ~ normal(alphapriors[1], alphapriors[2]); // sitecov shrinkage priors
   if(negbin == 1)
     phi ~ gamma(phipriors[1], phipriors[2]); // phi prior
 

diff --git a/inst/stan/traditional_continuous.stan b/inst/stan/traditional_continuous.stan
@@ -11,6 +11,9 @@ data{/////////////////////////////////////////////////////////////////////
     int<lower=0> nparams;  // number of gear types
     array[n_C] int<lower=1> mat;  // vector of gear type integers
     int<lower=0,upper=1> ctch; // binary indicator of presence of catchability coefficient
+    array[2] real betapriors; // priors for beta param in expected catch rate
+    array[2] real alphapriors; // priors for alpha param in expected catch rate
+
 
 }
 
@@ -43,8 +46,8 @@ model{/////////////////////////////////////////////////////////////////////
       E_trans[j] ~ gamma(lambda[j], beta[R_ind[j]]);  // Eq. 1.1
     }
 
-    beta ~ gamma(0.25,0.25);
-    alpha ~ gamma(0.01,0.01);
+    alpha ~ gamma(alphapriors[1], alphapriors[2]);
+    beta ~ gamma(betapriors[1], betapriors[2]);
 
 }
Original file line number	Diff line number	Diff line change
		@@ -1 +1 @@
		utils::globalVariables(c("L","median","survey_type","inits","value"))
		utils::globalVariables(c("L_ind","median","survey_type","inits","value"))