From 3db3f78312041c2e1cac0e1b80529af0e1b24485 Mon Sep 17 00:00:00 2001 From: abigailkeller Date: Tue, 31 Dec 2024 17:00:55 +0000 Subject: [PATCH] working on small reviewer changes --- .Rhistory | 700 +++++++++++------------ R/detectionPlot.R | 3 +- R/globals.R | 2 +- R/jointModel.R | 60 +- R/traditionalModel.R | 42 +- R/utils-joint.R | 14 +- R/utils-traditional.R | 14 +- eDNAjoint.Rproj | 1 + inst/stan/joint_continuous.stan | 9 +- inst/stan/joint_count.stan | 3 +- inst/stan/traditional_continuous.stan | 7 +- tests/testthat/test-detectionCalculate.R | 11 +- tests/testthat/test-detectionPlot.R | 8 +- tests/testthat/test-jointModel.R | 68 +-- tests/testthat/test-jointSelect.R | 15 +- tests/testthat/test-jointSummarize.R | 74 +-- tests/testthat/test-muCritical.R | 20 +- tests/testthat/test-traditionalModel.R | 20 +- 18 files changed, 543 insertions(+), 528 deletions(-) diff --git a/.Rhistory b/.Rhistory index 776a6c6..b5d6471 100644 --- a/.Rhistory +++ b/.Rhistory @@ -1,83 +1,286 @@ -model_data <- moddat[[1]] -inits <- moddat[[2]] -# run model -out <- rstan::sampling( -model, -data = model_data, -cores = 1, -#' @srrstats {G2.4,G2.4a} explicit conversion to -#' integers for sampling arguments -chains = n.chain, -thin = as.integer(thin), -warmup = as.integer(n.iter.burn), -iter = ( -as.integer(n.iter.burn) + as.integer(n.iter.sample) -), -init = inits, -refresh = ifelse(verbose == TRUE,500,0) +#traditional model, catchability coefficient +A <- list() +if(all(!is.null(initial_values))){ +for(i in 1:n.chain){ +A[[i]] <- list( +if('mu' %in% names(initial_values[[i]])){ +mu <- initial_values[[i]]$mu +} else { +mu <- stats::runif(length(unique(count_all$L_ind)), 0.01, 5) +}, +if('q' %in% names(initial_values[[i]])){ +q <- as.data.frame(initial_values[[i]]$q) +} else { +q <- as.data.frame(stats::runif(length(q_names),0.01,1)) +} ) -cov = NULL -family = 'negbin' -p10priors = c(1,20) -q = TRUE -phipriors = NULL -multicore = FALSE -initial_values = inits -n.chain = 1 -n.iter.burn = 100 -n.iter.sample = 75 -thin = 1 -adapt_delta = 0.9 -verbose = TRUE -seed = NULL -# moddat <- get_model_data(data,cov,family,p10priors,q,phipriors, -# multicore,initial_values,n.chain,n.iter.burn, -# n.iter.sample,thin,adapt_delta,verbose,seed) -moddat <- get_model_data_trad(data, family,q,phipriors,multicore, -initial_values,n.chain,n.iter.burn, -n.iter.sample,thin,adapt_delta, -verbose, seed) -inits -initial_values -## 13. -# model includes 'q','phi' (traditional model) -# constants -nsite <- 20 -nobs_count <- 100 -# params -mu <- rlnorm(nsite, meanlog = log(1), sdlog = 1) -phi <- 1.2 -q <- 2 -# traditional type -count_type <- cbind(matrix(1, nrow = nsite, ncol = nobs_count/2), -matrix(2, nrow = nsite, ncol = nobs_count/2)) -# count -count <- matrix(NA, nrow = nsite, ncol = nobs_count) -for(i in 1:nsite){ -for(j in 1:nobs_count){ -if(count_type[i,j]==1){ -count[i,j] <- rnbinom(n = 1, mu = mu[i], size = phi) +names(A[[i]]) <- c('mu','q') +} } else { -count[i,j] <- rnbinom(n = 1, mu = mu[i]*q, size = phi) +for(i in 1:n.chain){ +A[[i]] <- list( +mu <- stats::runif(length(unique(count_all$L_ind)), 0.01, 5), +q <- as.data.frame(stats::runif(length(q_names),0.01,1)) +) +names(A[[i]]) <- c('mu','q') } } +return(A) +} +#' @noRd +init_trad <- function(n.chain, count_all, initial_values){ +#helper function +#traditional model +A <- list() +if(all(!is.null(initial_values))){ +for(i in 1:n.chain){ +A[[i]] <- list( +if('mu' %in% names(initial_values[[i]])){ +mu <- initial_values[[i]]$mu +} else { +mu <- stats::runif(length(unique(count_all$L_ind)), 0.01, 5) } -# collect data -data <- list( -count = count, -count.type = count_type ) -# initial values -inits <- list() -inits[[1]] <- list( -mu = mu, -phi = phi +names(A[[i]]) <- 'mu' +} +} else { +for(i in 1:n.chain){ +A[[i]] <- list( +mu <- stats::runif(length(unique(count_all$L_ind)), 0.01, 5) ) -names(inits[[1]]) <- c('mu','phi') -cov = NULL -family = 'negbin' +names(A[[i]]) <- 'mu' +} +} +return(A) +} +# function for input checks +#' @srrstats {G5.2a} Pre-processing routines to check inputs have unique +#' messages +#' @noRd +traditionalModel_input_checks <- function(data, family, q, phipriors, n.chain, +n.iter.burn, n.iter.sample, +thin, adapt_delta, seed){ +## make sure all data tags are valid -- if q == TRUE +#' @srrstats {G2.13} Pre-processing routines to check for missing data +if (q == TRUE && !all(c('count.type','count') %in% names(data))){ +errMsg <- "Data should include 'count' and 'count.type'." +stop(errMsg) +} +## make sure all data tags are valid -- if q == FALSE +#' @srrstats {G2.13} Pre-processing routines to check for missing data +if (q == FALSE && !all(c('count') %in% names(data))){ +errMsg <- "Data should include 'count'." +stop(errMsg) +} +## make sure count is not zero-length +#' @srrstats {G5.8,G5.8a} Pre-processing routines to check for +#' zero-length data +if (dim(data$count)[1] == 0) { +errMsg <- "count contains zero-length data." +stop(errMsg) +} +## make sure no column is entirely NA in count +#' @srrstats {G5.8,G5.8c} Pre-processing routines to check for column +#' with all NA +if (any(apply(data$count, 2, function(col) all(is.na(col))))) { +errMsg <- "count contains a column with all NA." +stop(errMsg) +} +## make sure dimensions of count and count.type are equal, if count.type is +## present +#' @srrstats {BS2.1,G2.13} Pre-processing routines to ensure all input data +#' is dimensionally commensurate +if (q == TRUE){ +if(dim(data$count)[1] != dim(data$count.type)[1]| +dim(data$count)[2] != dim(data$count.type)[2]) { +errMsg <- "Dimensions of count and count.type do not match." +stop(errMsg) +} +} +## make sure no data are undefined +#' @srrstats {G2.16} Pre-processing routines to check for undefined data +if(any(data$count == Inf,na.rm=TRUE) | any(data$count == -Inf,na.rm=TRUE)){ +errMsg <- "count contains undefined values (i.e., Inf or -Inf)" +stop(errMsg) +} +## make sure all data is numeric -- if q == TRUE +#' @srrstats {BS2.5} Checks of appropriateness of numeric values submitted +#' for distributional parameters (i.e., count data must numeric), +#' implemented prior to analytic routines +#' @srrstats {G5.8,G5.8b} Pre-processing routines to check for data of +#' unsupported type +if (q == TRUE) { +if(is.numeric(data$count) == FALSE | +is.numeric(data$count.type) == FALSE) { +errMsg <- "Data should be numeric." +stop(errMsg) +} +} +## make sure all data is numeric -- if q == FALSE +#' @srrstats {BS2.5} Checks of appropriateness of numeric values submitted +#' for distributional parameters (i.e., count data must positive and +#' numeric), implemented prior to analytic routines +#' @srrstats {G5.8,G5.8b} Pre-processing routines to check for data of +#' unsupported type +if (q == FALSE) { +if(is.numeric(data$count) == FALSE | any(data$count < 0, na.rm=TRUE)) { +errMsg <- "Data should be numeric." +stop(errMsg) +} +} +if(q == TRUE){ +## make sure locations of NAs in count data match locations of NAs in +## count.type data +#' @srrstats {BS2.1,G2.13} Pre-processing routines to ensure all input +#' data is dimensionally commensurate +if(any((which(is.na(data$count.type)) == which(is.na(data$count))) == FALSE) +){ +errMsg <- paste0("Empty data cells (NA) in count data should match ", +"empty data cells (NA) in count.type data.") +stop(errMsg) +} +## make sure count.type is not zero-length +#' @srrstats {G5.8,G5.8a} Pre-processing routines to check for +#' zero-length data +if (dim(data$count.type)[1] == 0) { +errMsg <- "count.type contains zero-length data." +stop(errMsg) +} +## make sure no column is entirely NA in count.type +#' @srrstats {G5.8,G5.8c} Pre-processing routines to check for column +#' with all NA +if (any(apply(data$count.type, 2, function(col) all(is.na(col))))) { +errMsg <- "count.type contains a column with all NA." +stop(errMsg) +} +} +## make sure family is either 'poisson', 'negbin', or 'gamma' +#' @srrstats {G2.3,G2.3a,G2.3b} Permit only expected univariate +#' (case-insensitive) parameter values +if(!c(tolower(family) %in% c('poisson','negbin','gamma'))){ +errMsg <- "Invalid family. Options include 'poisson', 'negbin', or 'gamma'." +stop(errMsg) +} +## the smallest count.type is 1 +if(q == TRUE && min(data$count.type,na.rm=TRUE) != 1){ +errMsg <- paste0("The first gear type should be referenced as 1 in ", +"count.type. Subsequent gear types should be ", +"referenced 2, 3, 4, etc.") +stop(errMsg) +} +## count are integers, if family is poisson or negbin +#' @srrstats {BS2.5} Checks of appropriateness of numeric values submitted +#' for distributional parameters (i.e., count data must be non-negative +#' integers if a poisson or negative binomial distribution is used), +#' implemented prior to analytic routines +if(tolower(family) %in% c('poisson','negbin')){ +if(!all(data$count %% 1 %in% c(0,NA)) | any(data$count < 0,na.rm=TRUE)){ +errMsg <- paste0("All values in count should be non-negative integers. ", +"Use family = 'gamma' if count is continuous.") +stop(errMsg) +} +} +## count.type are integers +#' @srrstats {G5.8,G5.8b} Pre-processing routines to check for data of +#' unsupported type +if(q == TRUE && !all(data$count.type %% 1 %in% c(0,NA))){ +errMsg <- "All values in count.type should be integers." +stop(errMsg) +} +## phipriors is a vector of two numeric values +#' @srrstats {G2.0,BS2.2,BS2.3,BS2.4,BS2.5} Checks of vector length and +#' appropriateness of distributional parameters (i.e., vector of length 2, +#' numeric values > 0), implemented prior to analytic routines +if(family == 'negbin'){ +if(!is.numeric(phipriors) | length(phipriors) != 2 | any(phipriors <= 0)){ +errMsg <- paste0("phipriors should be a vector of two positive ", +"numeric values. ex. c(0.25,0.25)") +stop(errMsg) +} +} +## check length and range of n.chain +if(any(length(as.integer(n.chain)) > 1 | n.chain < 1)){ +errMsg <- "n.chain should be an integer > 0 and of length 1." +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." +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." +stop(errMsg) +} +## check length and range of thin +if(any(length(as.integer(thin)) > 1 | thin < 1)){ +errMsg <- "thin should be an integer > 0 and of length 1." +stop(errMsg) +} +## check length and range of adapt_delta +if(any(length(adapt_delta) > 1 | adapt_delta < 0 | adapt_delta > 1)){ +errMsg <- paste0("adapt_delta should be a numeric value > 0 and < 1 and ", +"of length 1.") +stop(errMsg) +} +## check length of seed +if(!is.null(seed)){ +if(length(as.integer(seed)) > 1){ +errMsg <- "seed should be an integer of length 1." +stop(errMsg) +} +} +} +# checks if initial values are provided +#' @noRd +initial_values_checks_trad <- function(initial_values,data,n.chain){ +## length of initial values is equal to the number of chains +if(length(initial_values) != n.chain){ +errMsg <- paste0("The length of the list of initial values should equal ", +"the number of chains (n.chain, default is 4).") +stop(errMsg) +} +for(i in 1:n.chain){ +## check mu input +if('mu' %in% names(initial_values[[i]])){ +## if mu is numeric +if(any(!is.numeric(initial_values[[i]]$mu)) | +any(initial_values[[i]]$mu < 0)){ +errMsg <- "Initial values for 'mu' should be numeric values > 0." +stop(errMsg) +} +## check mu length +if(length(initial_values[[i]]$mu) != dim(data$count)[1]){ +errMsg <- paste0("The length of initial values for 'mu' should ", +"equal the number of sites.") +stop(errMsg) +} +} +## check q input +if('q' %in% names(initial_values[[i]])){ +## if q is numeric +if(any(!is.numeric(initial_values[[i]]$q)) | +any(initial_values[[i]]$q < 0)){ +errMsg <- "Initial values for 'q' should be numeric." +stop(errMsg) +} +## check q length +if(length(initial_values[[i]]$q) != (length(table(data$count.type))-1)){ +errMsg <- paste0("The length of initial values for 'q' should equal: ", +"# unique gear types - 1 (i.e., q for reference ", +"type = 1).") +stop(errMsg) +} +} +} +} +source('../get_model_data.R') +source('../jointModel_helper.R') +cov = NULL +family = 'gamma' p10priors = c(1,20) -q = TRUE +q = FALSE phipriors = NULL multicore = FALSE initial_values = inits @@ -113,155 +316,36 @@ as.integer(n.iter.burn) + as.integer(n.iter.sample) init = inits, refresh = ifelse(verbose == TRUE,500,0) ) -rstan::summary(out,pars=c('mu','q'))$summary -model <- rstan::stan_model('inst/stan/traditional_count.stan') -# run model -out <- rstan::sampling( -model, -data = model_data, -cores = 1, -#' @srrstats {G2.4,G2.4a} explicit conversion to -#' integers for sampling arguments -chains = n.chain, -thin = as.integer(thin), -warmup = as.integer(n.iter.burn), -iter = ( -as.integer(n.iter.burn) + as.integer(n.iter.sample) -), -init = inits, -refresh = ifelse(verbose == TRUE,500,0) -) -rstan::summary(out,pars=c('mu','q'))$summary -test <- rstan::summary(out)$summary -View(test) -model <- rstan::stan_model('inst/stan/traditional_count.stan') -model <- rstan::stan_model('inst/stan/traditional_count.stan') -model <- rstan::stan_model('inst/stan/traditional_count.stan') -model <- rstan::stan_model('inst/stan/traditional_count.stan') -model <- rstan::stan_model('inst/stan/traditional_count.stan') -model <- rstan::stan_model('inst/stan/traditional_count.stan') -# run model -out <- rstan::sampling( -model, -data = model_data, -cores = 1, -#' @srrstats {G2.4,G2.4a} explicit conversion to -#' integers for sampling arguments -chains = n.chain, -thin = as.integer(thin), -warmup = as.integer(n.iter.burn), -iter = ( -as.integer(n.iter.burn) + as.integer(n.iter.sample) -), -init = inits, -refresh = ifelse(verbose == TRUE,500,0) -) -rstan::summary(out,pars=c('mu','q'))$summary -model_data$ctch -model_data$nparams -model <- rstan::stan_model('inst/stan/traditional_count.stan') -# run model -out <- rstan::sampling( -model, -data = model_data, -cores = 1, -#' @srrstats {G2.4,G2.4a} explicit conversion to -#' integers for sampling arguments -chains = n.chain, -thin = as.integer(thin), -warmup = as.integer(n.iter.burn), -iter = ( -as.integer(n.iter.burn) + as.integer(n.iter.sample) -), -init = inits, -refresh = ifelse(verbose == TRUE,500,0) -) -rstan::summary(out,pars=c('mu','q'))$summary +rstan::summary(out,pars=c('mu'))$summary model <- rstan::stan_model('inst/stan/traditional_count.stan') -# run model -out <- rstan::sampling( -model, -data = model_data, -cores = 1, -#' @srrstats {G2.4,G2.4a} explicit conversion to -#' integers for sampling arguments -chains = n.chain, -thin = as.integer(thin), -warmup = as.integer(n.iter.burn), -iter = ( -as.integer(n.iter.burn) + as.integer(n.iter.sample) -), -init = inits, -refresh = ifelse(verbose == TRUE,500,0) -) -rstan::summary(out,pars=c('mu','q'))$summary -model <- rstan::stan_model('inst/stan/traditional_count.stan') -# run model -out <- rstan::sampling( -model, -data = model_data, -cores = 1, -#' @srrstats {G2.4,G2.4a} explicit conversion to -#' integers for sampling arguments -chains = n.chain, -thin = as.integer(thin), -warmup = as.integer(n.iter.burn), -iter = ( -as.integer(n.iter.burn) + as.integer(n.iter.sample) -), -init = inits, -refresh = ifelse(verbose == TRUE,500,0) -) -rstan::summary(out,pars=c('mu','q'))$summary -model <- rstan::stan_model('inst/stan/traditional_count.stan') -model <- rstan::stan_model('inst/stan/traditional_count.stan') -# run model -out <- rstan::sampling( -model, -data = model_data, -cores = 1, -#' @srrstats {G2.4,G2.4a} explicit conversion to -#' integers for sampling arguments -chains = n.chain, -thin = as.integer(thin), -warmup = as.integer(n.iter.burn), -iter = ( -as.integer(n.iter.burn) + as.integer(n.iter.sample) -), -init = inits, -refresh = ifelse(verbose == TRUE,500,0) -) -test <- rstan::summary(out)$summary -View(test) -model <- rstan::stan_model('inst/stan/traditional_continuous.stan') +## 17. +# model, pois (traditional model) # constants nsite <- 20 nobs_count <- 100 # params mu <- rlnorm(nsite, meanlog = log(1), sdlog = 1) -q <- 2 -beta_gamma <- 1 -alpha_gamma <- mu * beta_gamma -# traditional type -count_type <- cbind(matrix(1, nrow = nsite, ncol = nobs_count/2), -matrix(2, nrow = nsite, ncol = nobs_count/2)) +# count +count <- matrix(NA, nrow = nsite, ncol = nobs_count) +for(i in 1:nsite){ +count[i,] <- rpois(n = nobs_count, mu[i]) +} # collect data data <- list( -count = count, -count.type = count_type +count = count ) # initial values inits <- list() inits[[1]] <- list( -alpha = mu, -beta = rep(1,length(mu)), -q = q +mu = mu ) -names(inits[[1]]) <- c('alpha','beta','q') -cov = NULL -family = 'gamma' +names(inits[[1]]) <- c('mu') +source('../get_model_data.R') +source('../jointModel_helper.R') +cov = NULL +family = 'poisson' p10priors = c(1,20) -q = TRUE +q = FALSE phipriors = NULL multicore = FALSE initial_values = inits @@ -297,12 +381,7 @@ as.integer(n.iter.burn) + as.integer(n.iter.sample) init = inits, refresh = ifelse(verbose == TRUE,500,0) ) -rstan::summary(out,pars=c('mu','q'))$summary -test <- rstan::summary(out)$summary -View(test) -model <- rstan::stan_model('inst/stan/joint_continuous.stan') -model <- rstan::stan_model('inst/stan/joint_continuous.stan') -model <- rstan::stan_model('inst/stan/joint_continuous.stan') +rstan::summary(out,pars=c('mu'))$summary # constants nsite <- 20 nobs_count <- 100 @@ -312,8 +391,7 @@ mu <- rlnorm(nsite, meanlog = log(1), sdlog = 1) beta <- 0.5 log_p10 <- -4.5 q <- 2 -beta_gamma <- 1 -alpha_gamma <- mu * beta_gamma +phi <- 1.2 # traditional type count_type <- cbind(matrix(1, nrow = nsite, ncol = nobs_count/2), matrix(2, nrow = nsite, ncol = nobs_count/2)) @@ -322,9 +400,9 @@ count <- matrix(NA, nrow = nsite, ncol = nobs_count) for(i in 1:nsite){ for(j in 1:nobs_count){ if(count_type[i,j]==1){ -count[i,j] <- rgamma(1,shape = alpha_gamma[i],rate = beta_gamma) +count[i,j] <- rnbinom(n = 1, mu = mu[i], size = phi) } else { -count[i,j] <- rgamma(1,shape = alpha_gamma[i]*q,rate = beta_gamma) +count[i,j] <- rnbinom(n = 1, mu = mu[i]*q, size = phi) } } } @@ -356,17 +434,16 @@ count.type = count_type # initial values inits <- list() inits[[1]] <- list( -alpha_gamma = mu, -beta_gamma = rep(1,length(mu)), +mu = mu, p10 = exp(log_p10), alpha = beta, +phi = phi, q = q ) -names(inits[[1]]) <- c('alpha_gamma','beta_gamma','p10','alpha','q') -source('../get_model_data.R') -source('../jointModel_helper.R') -cov = NULL -family = 'gamma' +names(inits[[1]]) <- c('mu','p10','alpha','phi','q') +source('R/utils-joint.R') +cov = NULL +family = 'negbin' p10priors = c(1,20) q = TRUE phipriors = NULL @@ -379,134 +456,57 @@ thin = 1 adapt_delta = 0.9 verbose = TRUE seed = NULL -moddat <- get_model_data(data,cov,family,p10priors,q,phipriors, -multicore,initial_values,n.chain,n.iter.burn, -n.iter.sample,thin,adapt_delta,verbose,seed) -# moddat <- get_model_data_trad(data, family,q,phipriors,multicore, -# initial_values,n.chain,n.iter.burn, -# n.iter.sample,thin,adapt_delta, -# verbose, seed) -model_data <- moddat[[1]] -inits <- moddat[[2]] -# run model -out <- rstan::sampling( -model, -data = model_data, -cores = 1, -#' @srrstats {G2.4,G2.4a} explicit conversion to -#' integers for sampling arguments -chains = n.chain, -thin = as.integer(thin), -warmup = as.integer(n.iter.burn), -iter = ( -as.integer(n.iter.burn) + as.integer(n.iter.sample) -), -init = inits, -refresh = ifelse(verbose == TRUE,500,0) -) -rstan::summary(out,pars=c('mu','q','p11_trad'))$summary -model <- rstan::stan_model('inst/stan/joint_count.stan') -## 2. -# model includes 'p10','beta','q' -# constants -nsite <- 20 -nobs_count <- 100 -nobs_pcr <- 8 -# params -mu <- rlnorm(nsite, meanlog = log(1), sdlog = 1) -beta <- 0.5 -log_p10 <- -4.5 -q <- 2 -# traditional type -count_type <- cbind(matrix(1, nrow = nsite, ncol = nobs_count/2), -matrix(2, nrow = nsite, ncol = nobs_count/2)) -# count -count <- matrix(NA, nrow = nsite, ncol = nobs_count) -for(i in 1:nsite){ -for(j in 1:nobs_count){ -if(count_type[i,j]==1){ -count[i,j] <- rpois(1, mu[i]) +# make character inputs case-insensitive +#' @srrstats {G2.3b} Allow case-insensitive character parameter values +family <- tolower(family) +# get phipriors +if(family != 'negbin'){ +phipriors <- NULL +} else if(family == 'negbin' && is.null(phipriors)){ +phipriors <- c(0.25,0.25) +} else if(family == 'negbin' && !is.null(phipriors)){ +phipriors <- phipriors +} +if (q == TRUE) { +# model with catchability coefficients +catchability_checks(data,cov) } else { -count[i,j] <- rpois(1, mu[i]*q) +# model without catchability coefficients +no_catchability_checks(data,cov) } +# model with covariates +if (all(!is.null(cov))) { +covariate_checks(data,cov) } +# all models +all_checks(data,cov,family,p10priors,phipriors,n.chain,n.iter.burn, +n.iter.sample,thin,adapt_delta,seed) +# initial value checks +if(all(!is.null(initial_values))){ +initial_values_checks(initial_values,data,cov,n.chain) } -# p11 (probability of true positive eDNA detection) and p (probability -# of eDNA detection) -p11 <- rep(NA,nsite) -p <- rep(NA,nsite) -for (i in 1:nsite){ -p11[i] <- mu[i] / (mu[i] + exp(beta)) -p[i] <- min(p11[i] + exp(log_p10),1) +if (!requireNamespace("rstan", quietly = TRUE)){ +stop ("The 'rstan' package is not installed.", call. = FALSE) } -# qPCR.N (# qPCR observations) -qPCR.N <- matrix(NA, nrow = nsite, ncol = nobs_pcr) -for(i in 1:nsite){ -qPCR.N[i,] <- rep(3,nobs_pcr) -} -# qPCR.K (# positive qPCR detections) -qPCR.K <- matrix(NA, nrow = nsite, ncol = nobs_pcr) -for (i in 1:nsite){ -qPCR.K[i,] <- rbinom(nobs_pcr, qPCR.N[i,], rep(p[i], nobs_pcr)) -} -# collect data -data <- list( -qPCR.N = qPCR.N, -qPCR.K = qPCR.K, -count = count, -count.type = count_type -) -# initial values -inits <- list() -inits[[1]] <- list( -mu = mu, -p10 = exp(log_p10), -alpha = beta, -q = q -) -names(inits[[1]]) <- c('mu','p10','alpha','q') -cov = NULL -family = 'poisson' -p10priors = c(1,20) -q = TRUE -phipriors = NULL -multicore = FALSE -initial_values = inits -n.chain = 1 -n.iter.burn = 100 -n.iter.sample = 75 -thin = 1 -adapt_delta = 0.9 -verbose = TRUE -seed = NULL -moddat <- get_model_data(data,cov,family,p10priors,q,phipriors, -multicore,initial_values,n.chain,n.iter.burn, -n.iter.sample,thin,adapt_delta,verbose,seed) -# moddat <- get_model_data_trad(data, family,q,phipriors,multicore, -# initial_values,n.chain,n.iter.burn, -# n.iter.sample,thin,adapt_delta, -# verbose, seed) -model_data <- moddat[[1]] -inits <- moddat[[2]] -# run model -out <- rstan::sampling( -model, -data = model_data, -cores = 1, -#' @srrstats {G2.4,G2.4a} explicit conversion to -#' integers for sampling arguments -chains = n.chain, -thin = as.integer(thin), -warmup = as.integer(n.iter.burn), -iter = ( -as.integer(n.iter.burn) + as.integer(n.iter.sample) -), -init = inits, -refresh = ifelse(verbose == TRUE,500,0) -) -rstan::summary(out,pars=c('mu','q','p11_trad'))$summary -model <- rstan::stan_model('inst/stan/joint_count.stan') -model <- rstan::stan_model('inst/stan/joint_count.stan') -model <- rstan::stan_model('inst/stan/joint_count.stan') -model <- rstan::stan_model('inst/stan/joint_count.stan') -model <- rstan::stan_model('inst/stan/joint_continuous.stan') +# 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') |> +#' @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() +# 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') |> +#' @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() +get_family_index('gamma') +roxygen2::roxygenise() +install.packages('roxygen2') +roxygen2::roxygenise() diff --git a/R/detectionPlot.R b/R/detectionPlot.R index d854778..84658c5 100644 --- 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 index 74e3eaa..b49ba00 100644 --- 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 index f199077..0dfd766 100644 --- 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,8 +172,8 @@ #' 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) #' } #' @@ -180,8 +181,8 @@ 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 index 071227d..48b7170 100644 --- 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,7 +121,7 @@ #' 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) #' } #' @@ -128,8 +129,8 @@ 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 index 8241bca..d382d07 100644 --- 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 index e563bb0..31bc23d 100644 --- 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 index 69fafd4..5442d74 100644 --- 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 index 541bad7..fb1a995 100644 --- a/inst/stan/joint_continuous.stan +++ b/inst/stan/joint_continuous.stan @@ -22,6 +22,9 @@ data{///////////////////////////////////////////////////////////////////// array[S_dna] int N_dna; // number of qPCR replicates per site of unpaired samples array[S_dna] int 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 nparams; // number of gear types array[n_C] int mat; // vector of gear type integers int 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 index 6952886..8cf2718 100644 --- a/inst/stan/joint_count.stan +++ b/inst/stan/joint_count.stan @@ -21,6 +21,7 @@ data{///////////////////////////////////////////////////////////////////// array[S_dna] int N_dna; // number of qPCR replicates per site of unpaired samples array[S_dna] int 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 nparams; // number of gear types array[n_C] int mat; // vector of gear type integers int 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 index 1a4904c..b15d8ac 100644 --- a/inst/stan/traditional_continuous.stan +++ b/inst/stan/traditional_continuous.stan @@ -11,6 +11,9 @@ data{///////////////////////////////////////////////////////////////////// int nparams; // number of gear types array[n_C] int mat; // vector of gear type integers int 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]); } diff --git a/tests/testthat/test-detectionCalculate.R b/tests/testthat/test-detectionCalculate.R index ac83c91..b429cbc 100644 --- a/tests/testthat/test-detectionCalculate.R +++ b/tests/testthat/test-detectionCalculate.R @@ -5,12 +5,13 @@ test_that("detectionCalculate input checks work", { # run joint model to do tests with model1 <- suppressWarnings({jointModel(data = gobyData, cov = c('Filter_time','Salinity'), - n.chain = 1,n.iter.burn = 25, - n.iter.sample = 75,multicore = FALSE)}) + n.chain = 1, n.warmup = 25, + n.iter = 100, multicore = FALSE)}) - model2 <- suppressWarnings({jointModel(data = greencrabData,family = 'negbin', - n.chain = 1,n.iter.burn = 25, - n.iter.sample = 75,multicore = FALSE)}) + model2 <- suppressWarnings({jointModel(data = greencrabData, + family = 'negbin', + n.chain = 1, n.warmup = 25, + n.iter = 100, multicore = FALSE)}) #1. make sure model fit is of class stanfit expect_error(detectionCalculate(as.matrix(model1$model), mu = c(0.1, 0.5), diff --git a/tests/testthat/test-detectionPlot.R b/tests/testthat/test-detectionPlot.R index 5789e5a..ffef0a6 100644 --- a/tests/testthat/test-detectionPlot.R +++ b/tests/testthat/test-detectionPlot.R @@ -5,14 +5,14 @@ test_that("detectionPlot input checks work", { # run joint model to do tests with model1 <- suppressWarnings({jointModel(data = gobyData, cov = c('Filter_time','Salinity'), - n.chain = 1,n.iter.burn = 25, - n.iter.sample = 75, + n.chain = 1,n.warmup = 25, + n.iter = 100, multicore = FALSE)}) model2 <- suppressWarnings({jointModel(data = greencrabData, family = 'negbin', - n.chain = 1,n.iter.burn = 25, - n.iter.sample = 75,multicore = FALSE)}) + n.chain = 1,n.warmup = 25, + n.iter = 100,multicore = FALSE)}) #1. make sure model fit is of class stanfit expect_error(detectionPlot(as.matrix(model1$model), mu.min = 0.1, diff --git a/tests/testthat/test-jointModel.R b/tests/testthat/test-jointModel.R index 566e5fb..a0541ef 100644 --- a/tests/testthat/test-jointModel.R +++ b/tests/testthat/test-jointModel.R @@ -740,36 +740,36 @@ test_that("jointModel input checks work", { n.chain = 0,multicore = FALSE), paste0("n.chain should be an integer > 0 and of length 1.")) - #48. check length and range of n.iter.sample + #48. check length and range of n.iter expect_error(jointModel(data = list(qPCR.K = rbind(c(1,1,1),c(1,1,NA)), qPCR.N = rbind(c(3,3,3),c(3,3,NA)), count = rbind(c(4,1,1),c(1,1,NA)), count.type = rbind(c(1,2,1),c(1,1,NA))), - n.iter.sample = c(1,1),multicore = FALSE), - paste0("n.iter.sample should be an integer > 0 and of length 1.") + n.iter = c(1,1), multicore = FALSE), + paste0("n.iter should be an integer > 0 and of length 1.") ) expect_error(jointModel(data = list(qPCR.K = rbind(c(1,1,1),c(1,1,NA)), qPCR.N = rbind(c(3,3,3),c(3,3,NA)), count = rbind(c(4,1,1),c(1,1,NA)), count.type = rbind(c(1,2,1),c(1,1,NA))), - n.iter.sample = 0,multicore = FALSE), - paste0("n.iter.sample should be an integer > 0 and of length 1.") + n.iter = 0,multicore = FALSE), + paste0("n.iter should be an integer > 0 and of length 1.") ) - #49. check length and range of n.iter.burn + #49. check length and range of n.warmup expect_error(jointModel(data = list(qPCR.K = rbind(c(1,1,1),c(1,1,NA)), qPCR.N = rbind(c(3,3,3),c(3,3,NA)), count = rbind(c(4,1,1),c(1,1,NA)), count.type = rbind(c(1,2,1),c(1,1,NA))), - n.iter.burn = c(1,1),multicore = FALSE), - paste0("n.iter.burn should be an integer > 0 and of length 1.") + n.warmup = c(1,1), multicore = FALSE), + paste0("n.warmup should be an integer > 0 and of length 1.") ) expect_error(jointModel(data = list(qPCR.K = rbind(c(1,1,1),c(1,1,NA)), qPCR.N = rbind(c(3,3,3),c(3,3,NA)), count = rbind(c(4,1,1),c(1,1,NA)), count.type = rbind(c(1,2,1),c(1,1,NA))), - n.iter.burn = 0,multicore = FALSE), - paste0("n.iter.burn should be an integer > 0 and of length 1.") + n.warmup = 0, multicore = FALSE), + paste0("n.warmup should be an integer > 0 and of length 1.") ) #50. check length and range of thin @@ -1287,8 +1287,8 @@ test_that("semi-paired model works", { # run model fit <- suppressWarnings({jointModel(data = data, q = TRUE, family = 'negbin', - initial_values = inits, n.iter.burn = 25, - n.iter.sample = 75, + initial_values = inits, n.warmup = 25, + n.iter = 100, n.chain = 1, multicore = FALSE, seed = 10 )}) @@ -1373,8 +1373,8 @@ test_that("semi-paired model works", { # run model fit <- suppressWarnings({jointModel(data = data, q = TRUE, n.chain = 1, multicore = FALSE, seed = 10, - initial_values = inits, n.iter.burn = 25, - n.iter.sample = 75)}) + initial_values = inits, n.warmup = 25, + n.iter = 100)}) # get output params output_params <- rownames(as.data.frame(jointSummarize(fit$model))) @@ -1459,8 +1459,8 @@ test_that("semi-paired model works", { # run model fit <- suppressWarnings({jointModel(data = data, q = TRUE, family = 'gamma', n.chain = 1, multicore = FALSE, seed = 10, - initial_values = inits, n.iter.burn = 25, - n.iter.sample = 75)}) + initial_values = inits, n.warmup = 25, + n.iter = 100)}) # get output params output_params <- rownames(as.data.frame(jointSummarize(fit$model))) @@ -1531,8 +1531,8 @@ test_that("semi-paired model works", { fit <- suppressWarnings({jointModel(data = data, family = 'negbin', initial_values = inits, n.chain = 1, multicore = FALSE, seed = 10, - n.iter.burn = 25, - n.iter.sample = 75)}) + n.warmup = 25, + n.iter = 100)}) # get output params output_params <- rownames(as.data.frame(jointSummarize(fit$model))) @@ -1602,8 +1602,8 @@ test_that("semi-paired model works", { # run model fit <- suppressWarnings({jointModel(data = data, initial_values = inits, n.chain = 1, multicore = FALSE, seed = 10, - n.iter.burn = 25, - n.iter.sample = 75)}) + n.warmup = 25, + n.iter = 100)}) # get output params output_params <- rownames(as.data.frame(jointSummarize(fit$model))) @@ -1677,8 +1677,8 @@ test_that("semi-paired model works", { fit <- suppressWarnings({jointModel(data = data, initial_values = inits, family = 'gamma', n.chain = 1, multicore = FALSE, seed = 10, - n.iter.burn = 25, - n.iter.sample = 75)}) + n.warmup = 25, + n.iter = 100)}) # get output params output_params <- rownames(as.data.frame(jointSummarize(fit$model))) @@ -1774,8 +1774,8 @@ test_that("semi-paired model works", { cov = c('var_a','var_b'), n.chain = 1, multicore = FALSE, seed = 10, initial_values = inits, adapt_delta = 0.99, - n.iter.burn = 25, - n.iter.sample = 75)}) + n.warmup = 25, + n.iter = 100)}) # get output params output_params <- rownames(as.data.frame(jointSummarize(fit$model))) @@ -1867,8 +1867,8 @@ test_that("semi-paired model works", { fit <- suppressWarnings({jointModel(data = data, q = TRUE, cov = c('var_a','var_b'), n.chain = 1, multicore = FALSE, seed = 10, - initial_values = inits, n.iter.burn = 25, - n.iter.sample = 75)}) + initial_values = inits, n.warmup = 25, + n.iter = 100)}) # get output params output_params <- rownames(as.data.frame(jointSummarize(fit$model))) @@ -1965,8 +1965,8 @@ test_that("semi-paired model works", { fit <- suppressWarnings({jointModel(data = data, q = TRUE, family = 'gamma', cov = c('var_a','var_b'), n.chain = 1, multicore = FALSE, seed = 10, - initial_values = inits, n.iter.burn = 25, - n.iter.sample = 75 + initial_values = inits, n.warmup = 25, + n.iter = 100 )}) # get output params @@ -2047,8 +2047,8 @@ test_that("semi-paired model works", { 'phi') # run model fit <- suppressWarnings({jointModel(data = data, family = 'negbin', - cov = c('var_a','var_b'),n.iter.burn = 25, - n.iter.sample = 75, + cov = c('var_a','var_b'), n.warmup = 25, + n.iter = 100, n.chain = 1, multicore = FALSE, seed = 10, initial_values = inits)}) @@ -2130,8 +2130,8 @@ test_that("semi-paired model works", { fit <- suppressWarnings({jointModel(data = data, cov = c('var_a','var_b'), n.chain = 1, multicore = FALSE, seed = 10, - initial_values = inits, n.iter.burn = 25, - n.iter.sample = 75)}) + initial_values = inits, n.warmup = 25, + n.iter = 100)}) # get output params output_params <- rownames(as.data.frame(jointSummarize(fit$model))) @@ -2216,8 +2216,8 @@ test_that("semi-paired model works", { fit <- suppressWarnings({jointModel(data = data,family = 'gamma', cov = c('var_a','var_b'), n.chain = 1, multicore = FALSE, seed = 10, - initial_values = inits, n.iter.burn = 25, - n.iter.sample = 75)}) + initial_values = inits, n.warmup = 25, + n.iter = 100)}) # get output params output_params <- rownames(as.data.frame(jointSummarize(fit$model))) diff --git a/tests/testthat/test-jointSelect.R b/tests/testthat/test-jointSelect.R index 7ba284f..795e49e 100644 --- a/tests/testthat/test-jointSelect.R +++ b/tests/testthat/test-jointSelect.R @@ -8,13 +8,14 @@ test_that("jointSelect input checks work", { model1 <- suppressWarnings({traditionalModel(data = greencrabData, family = 'negbin', - multicore = FALSE, n.chain = 1,n.iter.burn = 25, - n.iter.sample = 75)}) + multicore = FALSE, n.chain = 1, + n.warmup = 25, + n.iter = 100)}) model2 <- suppressWarnings({jointModel(data = greencrabData, family = 'negbin', - multicore = FALSE,n.chain = 1, - n.iter.burn = 25, - n.iter.sample = 75)}) + multicore = FALSE, n.chain = 1, + n.warmup = 25, + n.iter = 100)}) # 1. Check that model inputs are a list @@ -37,10 +38,10 @@ test_that("jointSelect output check", { # fit models fit.q1 <- jointModel(data = greencrabData, family = "negbin", p10priors = c(1,20), q = TRUE, multicore = FALSE, - n.chain = 1, n.iter.sample = 1000) + n.chain = 1, n.iter = 1000) fit.q2 <- jointModel(data = greencrabData, family = "negbin", p10priors = c(1,50), q = TRUE, multicore = FALSE, - n.chain = 1, n.iter.sample = 1000) + n.chain = 1, n.iter = 1000) # run select select <- jointSelect(modelfits = list(fit.q1$model, fit.q2$model)) diff --git a/tests/testthat/test-jointSummarize.R b/tests/testthat/test-jointSummarize.R index fb1fbd7..f20866a 100644 --- a/tests/testthat/test-jointSummarize.R +++ b/tests/testthat/test-jointSummarize.R @@ -8,7 +8,7 @@ test_that("jointSummarize input checks work", { out <- traditionalModel(data = greencrabData, family = 'negbin', multicore = FALSE, - n.chain = 1, n.iter.sample = 1000) + n.chain = 1, n.iter = 1000) #1. make sure model fit is of class stanfit data <- data.frame(y = c(1,2,3),x = c(1,2,3)) @@ -103,8 +103,8 @@ test_that("jointSummarize outputs work", { # run model fit <- suppressWarnings({jointModel(data = data, q = TRUE, family = 'negbin', - initial_values = inits, n.iter.burn = 25, - n.iter.sample = 75, + initial_values = inits, n.warmup = 25, + n.iter = 100, n.chain = 1, multicore = FALSE, seed = 10 )}) @@ -210,8 +210,8 @@ test_that("jointSummarize outputs work", { # run model fit <- suppressWarnings({jointModel(data = data, q = TRUE, n.chain = 1, multicore = FALSE, seed = 10, - initial_values = inits, n.iter.burn = 25, - n.iter.sample = 75)}) + initial_values = inits, n.warmup = 25, + n.iter = 100)}) # get output params output_params <- rownames(as.data.frame(jointSummarize(fit$model))) @@ -317,8 +317,8 @@ test_that("jointSummarize outputs work", { # run model fit <- suppressWarnings({jointModel(data = data, q = TRUE, family = 'gamma', n.chain = 1, multicore = FALSE, seed = 10, - initial_values = inits, n.iter.burn = 25, - n.iter.sample = 75)}) + initial_values = inits, n.warmup = 25, + n.iter = 100)}) # get output params output_params <- rownames(as.data.frame(jointSummarize(fit$model))) @@ -380,8 +380,8 @@ test_that("jointSummarize outputs work", { fit <- suppressWarnings({jointModel(data = data, family = 'negbin', initial_values = inits, n.chain = 1, multicore = FALSE, seed = 10, - n.iter.burn = 25, - n.iter.sample = 75)}) + n.warmup = 25, + n.iter = 100)}) # get output params output_params <- rownames(as.data.frame(jointSummarize(fit$model))) @@ -471,8 +471,8 @@ test_that("jointSummarize outputs work", { # run model fit <- suppressWarnings({jointModel(data = data, initial_values = inits, n.chain = 1, multicore = FALSE, seed = 10, - n.iter.burn = 25, - n.iter.sample = 75)}) + n.warmup = 25, + n.iter = 100)}) # get output params output_params <- rownames(as.data.frame(jointSummarize(fit$model))) @@ -566,8 +566,8 @@ test_that("jointSummarize outputs work", { fit <- suppressWarnings({jointModel(data = data, initial_values = inits, family = 'gamma', n.chain = 1, multicore = FALSE, seed = 10, - n.iter.burn = 25, - n.iter.sample = 75)}) + n.warmup = 25, + n.iter = 100)}) # get output params output_params <- rownames(as.data.frame(jointSummarize(fit$model))) @@ -653,8 +653,8 @@ test_that("jointSummarize outputs work", { cov = c('var_a','var_b'), n.chain = 1, multicore = FALSE, seed = 10, initial_values = inits, adapt_delta = 0.99, - n.iter.burn = 25, - n.iter.sample = 75)}) + n.warmup = 25, + n.iter = 100)}) # get output params output_params <- rownames(as.data.frame(jointSummarize(fit$model))) @@ -769,8 +769,8 @@ test_that("jointSummarize outputs work", { fit <- suppressWarnings({jointModel(data = data, q = TRUE, cov = c('var_a','var_b'), n.chain = 1, multicore = FALSE, seed = 10, - initial_values = inits, n.iter.burn = 25, - n.iter.sample = 75)}) + initial_values = inits, n.warmup = 25, + n.iter = 100)}) # get output params output_params <- rownames(as.data.frame(jointSummarize(fit$model))) @@ -889,8 +889,8 @@ test_that("jointSummarize outputs work", { fit <- suppressWarnings({jointModel(data = data, q = TRUE, family = 'gamma', cov = c('var_a','var_b'), n.chain = 1, multicore = FALSE, seed = 10, - initial_values = inits, n.iter.burn = 25, - n.iter.sample = 75 + initial_values = inits, n.warmup = 25, + n.iter = 100 )}) # get output params @@ -963,8 +963,8 @@ test_that("jointSummarize outputs work", { 'phi') # run model fit <- suppressWarnings({jointModel(data = data, family = 'negbin', - cov = c('var_a','var_b'),n.iter.burn = 25, - n.iter.sample = 75, + cov = c('var_a','var_b'),n.warmup = 25, + n.iter = 100, n.chain = 1, multicore = FALSE, seed = 10, initial_values = inits)}) @@ -1068,8 +1068,8 @@ test_that("jointSummarize outputs work", { fit <- suppressWarnings({jointModel(data = data, cov = c('var_a','var_b'), n.chain = 1, multicore = FALSE, seed = 10, - initial_values = inits, n.iter.burn = 25, - n.iter.sample = 75)}) + initial_values = inits, n.warmup = 25, + n.iter = 100)}) # get output params output_params <- rownames(as.data.frame(jointSummarize(fit$model))) @@ -1174,8 +1174,8 @@ test_that("jointSummarize outputs work", { fit <- suppressWarnings({jointModel(data = data,family = 'gamma', cov = c('var_a','var_b'), n.chain = 1, multicore = FALSE, seed = 10, - initial_values = inits, n.iter.burn = 25, - n.iter.sample = 75)}) + initial_values = inits, n.warmup = 25, + n.iter = 100)}) # get output params output_params <- rownames(as.data.frame(jointSummarize(fit$model))) @@ -1256,8 +1256,8 @@ test_that("jointSummarize outputs work", { n.chain = 1, multicore = FALSE, seed = 10, initial_values = inits, - n.iter.burn = 25, - n.iter.sample = 75)}) + n.warmup = 25, + n.iter = 100)}) # get output params output_params <- rownames(as.data.frame(jointSummarize(fit$model))) @@ -1335,8 +1335,8 @@ test_that("jointSummarize outputs work", { # run model fit <- suppressWarnings({traditionalModel(data = data, q = TRUE, n.chain = 1, multicore = FALSE, seed = 10, - initial_values = inits, n.iter.burn = 25, - n.iter.sample = 75)}) + initial_values = inits, n.warmup = 25, + n.iter = 100)}) # get output params output_params <- rownames(as.data.frame(jointSummarize(fit$model))) @@ -1408,8 +1408,8 @@ test_that("jointSummarize outputs work", { n.chain = 1, multicore = FALSE, seed = 10, initial_values = inits, - n.iter.burn = 25, - n.iter.sample = 75)}) + n.warmup = 25, + n.iter = 100)}) # get output params output_params <- rownames(as.data.frame(jointSummarize(fit$model))) @@ -1478,8 +1478,8 @@ test_that("jointSummarize outputs work", { n.chain = 1, multicore = FALSE, seed = 10, initial_values = inits, - n.iter.burn = 25, - n.iter.sample = 75)}) + n.warmup = 25, + n.iter = 100)}) # get output params output_params <- rownames(as.data.frame(jointSummarize(fit$model))) @@ -1546,8 +1546,8 @@ test_that("jointSummarize outputs work", { fit <- suppressWarnings({traditionalModel(data = data,n.chain = 1, multicore = FALSE, seed = 10, initial_values = inits, - n.iter.burn = 25, - n.iter.sample = 75)}) + n.warmup = 25, + n.iter = 100)}) # get output params output_params <- rownames(as.data.frame(jointSummarize(fit$model))) @@ -1615,8 +1615,8 @@ test_that("jointSummarize outputs work", { fit <- suppressWarnings({traditionalModel(data = data, n.chain = 1, family = 'gamma', multicore = FALSE, seed = 10, - n.iter.burn = 25, - n.iter.sample = 75, + n.warmup = 25, + n.iter = 100, initial_values = inits)}) # get output params diff --git a/tests/testthat/test-muCritical.R b/tests/testthat/test-muCritical.R index b05c7e5..3cb49e0 100644 --- a/tests/testthat/test-muCritical.R +++ b/tests/testthat/test-muCritical.R @@ -7,19 +7,19 @@ test_that("muCritical input checks work", { model1 <- suppressWarnings({jointModel(data = gobyData, cov = c('Filter_time','Salinity'), multicore = FALSE, n.chain = 1, - n.iter.burn = 25, - n.iter.sample = 75)}) + n.warmup = 25, + n.iter = 100)}) model2 <- suppressWarnings({traditionalModel(data = gobyData, multicore = FALSE, - n.chain = 1, n.iter.burn = 25, - n.iter.sample = 75)}) + n.chain = 1, n.warmup = 25, + n.iter = 100)}) model3 <- suppressWarnings({jointModel(data = greencrabData, family = 'negbin', multicore = FALSE, - n.chain = 1, n.iter.burn = 25, - n.iter.sample = 75)}) + n.chain = 1, n.warmup = 25, + n.iter = 100)}) #1. make sure model fit is of class stanfit expect_error(muCritical(as.matrix(model1$model), cov.val = c(0,0)), @@ -64,13 +64,13 @@ test_that("muCritical output check", { p10priors = c(1,20), q = FALSE, multicore = FALSE, n.chain = 1, - n.iter.burn = 25, - n.iter.sample = 75)}) + n.warmup = 25, + n.iter = 100)}) fit.q <- suppressWarnings({jointModel(data = greencrabData, cov = NULL, family = "negbin",p10priors = c(1,20), q = TRUE, multicore = FALSE, n.chain = 1, - n.iter.burn = 25, - n.iter.sample = 75)}) + n.warmup = 25, + n.iter = 100)}) # check lengths of muCritical output expect_equal(length(muCritical(fit.cov$model, cov.val = c(0,0), ci = 0.9)), diff --git a/tests/testthat/test-traditionalModel.R b/tests/testthat/test-traditionalModel.R index 5e72694..b777bb1 100644 --- a/tests/testthat/test-traditionalModel.R +++ b/tests/testthat/test-traditionalModel.R @@ -176,32 +176,32 @@ test_that("traditionalModel input checks work", { n.chain = 0, multicore = FALSE), paste0("n.chain should be an integer > 0 and of length 1.")) - #19. check length and range of n.iter.sample + #19. check length and range of n.iter expect_error(traditionalModel(data = list(count = rbind(c(4,1,1),c(1,1,NA)), count.type = rbind(c(1,2,1), c(1,1,NA))), - n.iter.sample = c(1,1), multicore = FALSE), - paste0("n.iter.sample should be an integer > 0 and of length 1.") + n.iter = c(1,1), multicore = FALSE), + paste0("n.iter should be an integer > 0 and of length 1.") ) expect_error(traditionalModel(data = list(count = rbind(c(4,1,1),c(1,1,NA)), count.type = rbind(c(1,2,1), c(1,1,NA))), - n.iter.sample = 0, multicore = FALSE), - paste0("n.iter.sample should be an integer > 0 and of length 1.") + n.iter = 0, multicore = FALSE), + paste0("n.iter should be an integer > 0 and of length 1.") ) - #20. check length and range of n.iter.burn + #20. check length and range of n.warmup expect_error(traditionalModel(data = list(count = rbind(c(4,1,1),c(1,1,NA)), count.type = rbind(c(1,2,1), c(1,1,NA))), - n.iter.burn = c(1,1), multicore = FALSE), - paste0("n.iter.burn should be an integer > 0 and of length 1.") + n.warmup = c(1,1), multicore = FALSE), + paste0("n.warmup should be an integer > 0 and of length 1.") ) expect_error(traditionalModel(data = list(count = rbind(c(4,1,1),c(1,1,NA)), count.type = rbind(c(1,2,1), c(1,1,NA))), - n.iter.burn = 0, multicore = FALSE), - paste0("n.iter.burn should be an integer > 0 and of length 1.") + n.warmup = 0, multicore = FALSE), + paste0("n.warmup should be an integer > 0 and of length 1.") ) #21. check length and range of thin