From 12ddbd97c59e528d24d1aa32328a4b0209c4791a Mon Sep 17 00:00:00 2001 From: Abby Date: Tue, 13 Aug 2024 09:37:09 -0400 Subject: [PATCH] updated lifecycle badge --- .Rhistory | 976 ++++++++++++++++++++++++++--------------------------- README.Rmd | 2 +- README.md | 6 +- 3 files changed, 492 insertions(+), 492 deletions(-) diff --git a/.Rhistory b/.Rhistory index 147e4b3..791ef6f 100644 --- a/.Rhistory +++ b/.Rhistory @@ -1,512 +1,512 @@ -# collect data -data <- list( -qPCR.N = qPCR.N, -qPCR.K = qPCR.K, -count = count -) -# initial values +#24. site.cov has same number of rows as qPCR.N and count, if present +site.cov <- cbind(c(0,1,1),c(0.4,-0.4,1)) +colnames(site.cov) <- c('var_a','var_b') +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,2,NA)), +site.cov = site.cov), +cov = c('var_a','var_b'),q = TRUE, +multicore = FALSE), +paste(paste0("The number of rows in site.cov matrix should match ", +"the number of rows in all other matrices."), +'See the eDNAjoint guide for data formatting help: ', +paste0('https://bookdown.org/abigailkeller/eDNAjoint_', +'vignette/usecase2.html'), +sep='\n')) +#25. make sure count.type is not zero-length +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 = matrix(NA,ncol = 3, +nrow = 0)), +q = TRUE, +multicore = FALSE), +paste("count.type contains zero-length data.", +'See the eDNAjoint guide for data formatting help: ', +paste0('https://bookdown.org/abigailkeller/eDNAjoint_', +'vignette/usecase3.html#prepare-the-data'), +sep='\n')) +#26. make sure no column is entirely NA in count.type +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(4,1,NA),c(1,1,NA))), +q = TRUE, +multicore = FALSE), +paste("count.type contains a column with all NA.", +'See the eDNAjoint guide for data formatting help: ', +paste0('https://bookdown.org/abigailkeller/eDNAjoint_', +'vignette/usecase3.html#prepare-the-data'), +sep='\n')) +#27. make sure no column is entirely NA in qPCR.K +expect_error(jointModel(data = list(qPCR.K = rbind(c(1,1,NA),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))), +multicore = FALSE), +paste("qPCR.K contains a column with all NA.", +'See the eDNAjoint guide for data formatting help: ', +paste0('https://bookdown.org/abigailkeller/eDNAjoint_', +'vignette/usecase1.html#prepare-the-data'), +sep='\n')) +#28. make sure no column is entirely NA in qPCR.N +expect_error(jointModel(data = list(qPCR.K = rbind(c(1,1,1),c(1,1,NA)), +qPCR.N = rbind(c(3,3,NA),c(3,3,NA)), +count = rbind(c(4,1,1),c(1,1,NA))), +multicore = FALSE), +paste("qPCR.N contains a column with all NA.", +'See the eDNAjoint guide for data formatting help: ', +paste0('https://bookdown.org/abigailkeller/eDNAjoint_', +'vignette/usecase1.html#prepare-the-data'), +sep='\n')) +#29. make sure no column is entirely NA in count +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,NA),c(1,1,NA))), +multicore = FALSE), +paste("count contains a column with all NA.", +'See the eDNAjoint guide for data formatting help: ', +paste0('https://bookdown.org/abigailkeller/eDNAjoint_', +'vignette/usecase1.html#prepare-the-data'), +sep='\n')) +#30. make sure no data are undefined +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,Inf),c(1,1,NA))), +multicore = FALSE), +paste("count contains undefined values \\(i.e., Inf or -Inf\\)", +'See the eDNAjoint guide for data formatting help: ', +paste0('https://bookdown.org/abigailkeller/eDNAjoint_', +'vignette/usecase1.html#prepare-the-data'), +sep='\n')) +#31. make sure no data are undefined +expect_error(jointModel(data = list(qPCR.K = rbind(c(1,1,1),c(1,1,NA)), +qPCR.N = rbind(c(3,3,Inf),c(3,3,NA)), +count = rbind(c(4,1,1),c(1,1,NA))), +multicore = FALSE), +paste("qPCR.N contains undefined values \\(i.e., Inf or -Inf\\)", +'See the eDNAjoint guide for data formatting help: ', +paste0('https://bookdown.org/abigailkeller/eDNAjoint_', +'vignette/usecase1.html#prepare-the-data'), +sep='\n')) +#32. make sure no data are undefined +expect_error(jointModel(data = list(qPCR.K = rbind(c(1,1,Inf),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))), +multicore = FALSE), +paste("qPCR.K contains undefined values \\(i.e., Inf or -Inf\\)", +'See the eDNAjoint guide for data formatting help: ', +paste0('https://bookdown.org/abigailkeller/eDNAjoint_', +'vignette/usecase1.html#prepare-the-data'), +sep='\n')) +#33. make sure site.cov is not zero-length +site.cov <- matrix(NA,ncol = 2,nrow = 0) +colnames(site.cov) <- c('var_a','var_b') +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)), +site.cov = site.cov), +cov = c('var_a','var_b'), +multicore = FALSE), +paste("site.cov contains zero-length data.", +'See the eDNAjoint guide for data formatting help: ', +paste0('https://bookdown.org/abigailkeller/eDNAjoint_', +'vignette/usecase2.html'), +sep='\n')) +#34. make sure no column is entirely NA in site.cov +site.cov <- rbind(c(4,1,NA),c(1,1,NA)) +colnames(site.cov) <- c('var_a','var_b','var_c') +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)), +site.cov = site.cov), +cov = c('var_a','var_b'), +multicore = FALSE), +paste("site.cov contains a column with all NA.", +'See the eDNAjoint guide for data formatting help: ', +paste0('https://bookdown.org/abigailkeller/eDNAjoint_', +'vignette/usecase2.html'), +sep='\n')) +#35. make sure no data are undefined +site.cov <- rbind(c(4,1,Inf),c(1,1,NA)) +colnames(site.cov) <- c('var_a','var_b','var_c') +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)), +site.cov = site.cov), +cov = c('var_a','var_b'), +multicore = FALSE), +paste(paste0("site.cov contains undefined values \\(i.e., ", +"Inf or -Inf\\)"), +'See the eDNAjoint guide for data formatting help: ', +paste0('https://bookdown.org/abigailkeller/eDNAjoint_', +'vignette/usecase2.html'), +sep='\n')) +#36. length of initial values is equal to the number of chains +n.chain <- 4 inits <- list() -inits[[1]] <- list( -alpha_gamma = mu, -beta_gamma = rep(1,length(mu)), -p10 = log_p10, -beta = beta +for(i in 1:n.chain){ +inits[[i]] <- list( +mu <- stats::runif(3, 0.01, 5), +p10 <- stats::runif(1,log(0.0001),log(0.08)), +alpha <- rep(0.1,3) ) -names(inits[[1]]) <- c('alpha_gamma','beta_gamma','p10','beta') -# run model -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)}) -# get output params -output_params <- rownames(as.data.frame(jointSummarize(fit$model))) -# test expectation -expect_true(all(c('p10','beta') %in% output_params)) -# check length of of mu -mu_estimates <- as.vector(rstan::summary(fit$model,par='mu')$summary[,1]) -expect_equal(length(mu_estimates),nsite) -expect_true(is.numeric(mu_estimates)) -## 7. -# model includes 'p10','q','phi','alpha' -# constants -nsite <- 20 -nobs_count <- 100 -nobs_pcr <- 8 -# params -mu <- rlnorm(nsite, meanlog = log(1), sdlog = 1) -alpha <- c(0.5, 0.1, -0.4) -log_p10 <- -4.5 -q <- 2 -phi <- 10 -# 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) -} else { -count[i,j] <- rnbinom(n = 1, mu = mu[i]*q, size = phi) -} -} -} -# site-level covariates -mat_site <- matrix(NA, nrow = nsite, ncol = length(alpha)) -mat_site[,1] <- 1 # intercept -for(i in 2:length(alpha)){ -mat_site[,i] <- rnorm(nsite,0,1) -} -colnames(mat_site) <- c('int','var_a','var_b') -# 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(sum(mat_site[i,]*alpha))) -p[i] <- min(p11[i] + exp(log_p10),1) } -# 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)) -} -# add NA to sites 2 and 4 -count[c(2,4),] <- NA -count_type[c(2,4),] <- NA -# collect data -data <- list( -qPCR.N = qPCR.N, -qPCR.K = qPCR.K, -count = count, -count.type = count_type, -site.cov = mat_site -) -# initial values +site.cov <- rbind(c(4,1),c(1,1)) +colnames(site.cov) <- c('var_a','var_b') +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)), +site.cov = site.cov), +cov = c('var_a','var_b'), initial_values = inits, +n.chain = 5, +multicore = FALSE), +paste(paste0("The length of the list of initial values ", +"should equal the number of chains \\(n.chain, ", +"default is 4\\)."), +paste0('See the eDNAjoint guide for help formatting ', +'initial values: '), +paste0('https://bookdown.org/abigailkeller/eDNAjoint_', +'vignette/usecase1.html#initialvalues'), +sep='\n')) +#37. initial values check: if mu is numeric +n.chain <- 4 inits <- list() -inits[[1]] <- list( -mu = mu, -p10 = log_p10, -alpha = alpha, -q = q, -phi = phi +for(i in 1:n.chain){ +inits[[i]] <- list( +mu <- stats::runif(3, -1, 0), +p10 <- stats::runif(1,log(0.0001),log(0.08)), +alpha <- rep(0.1,3) ) -names(inits[[1]]) <- c('mu','p10','alpha','q','phi' +names(inits[[i]]) <- c('mu','p10','alpha') +} +site.cov <- rbind(c(4,1),c(1,1)) +colnames(site.cov) <- c('var_a','var_b') +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)), +site.cov = site.cov), +cov = c('var_a','var_b'), initial_values = inits, +multicore = FALSE), +paste("Initial values for 'mu' should be numeric values > 0.", +paste0('See the eDNAjoint guide for help formatting ', +'initial values: '), +paste0('https://bookdown.org/abigailkeller/eDNAjoint_', +'vignette/usecase1.html#initialvalues'), +sep='\n')) +#38. initial values check: mu length +n.chain <- 4 +inits <- list() +for(i in 1:n.chain){ +inits[[i]] <- list( +mu <- stats::runif(4, 0.1, 1), +p10 <- stats::runif(1,log(0.0001),log(0.08)), +alpha <- rep(0.1,3) ) -# run model -fit <- suppressWarnings({jointModel(data = data, family = 'negbin', q = TRUE, -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)}) -# get output params -output_params <- rownames(as.data.frame(jointSummarize(fit$model))) -# test expectation -expect_true(all(c('p10','q[1]','phi','alpha[1]') %in% output_params)) -# check length of of mu -mu_estimates <- as.vector(rstan::summary(fit$model,par='mu')$summary[,1]) -expect_equal(length(mu_estimates),nsite*2) -expect_true(is.numeric(mu_estimates)) -## 8. -# model includes 'p10','alpha','q' -# constants -nsite <- 20 -nobs_count <- 100 -nobs_pcr <- 8 -# params -mu <- rlnorm(nsite, meanlog = log(1), sdlog = 1) -alpha <- c(0.5, 0.1, -0.4) -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(n = 1, mu[i]) -} else { -count[i,j] <- rpois(n = 1, mu[i]*q) -} -} -} -# site-level covariates -mat_site <- matrix(NA, nrow = nsite, ncol = length(alpha)) -mat_site[,1] <- 1 # intercept -for(i in 2:length(alpha)){ -mat_site[,i] <- rnorm(nsite,0,1) -} -colnames(mat_site) <- c('int','var_a','var_b') -# 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(sum(mat_site[i,]*alpha))) -p[i] <- min(p11[i] + exp(log_p10),1) -} -# 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)) -} -# add NA to sites 2 and 4 -count[c(2,4),] <- NA -count_type[c(2,4),] <- NA -# collect data -data <- list( -qPCR.N = qPCR.N, -qPCR.K = qPCR.K, -count = count, -count.type = count_type, -site.cov = mat_site +names(inits[[i]]) <- c('mu','p10','alpha') +} +site.cov <- rbind(c(4,1),c(1,1)) +colnames(site.cov) <- c('var_a','var_b') +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)), +site.cov = site.cov), +cov = c('var_a','var_b'), initial_values = inits, +multicore = FALSE), +paste(paste0("The length of initial values for 'mu' should ", +"equal the number of sites."), +paste0('See the eDNAjoint guide for help formatting ', +'initial values: '), +paste0('https://bookdown.org/abigailkeller/eDNAjoint_', +'vignette/usecase1.html#initialvalues'), +sep='\n')) +#39. initial values check: p10 numeric +n.chain <- 4 +inits <- list() +for(i in 1:n.chain){ +inits[[i]] <- list( +mu <- stats::runif(2,0,1), +p10 <- "-1", +alpha <- rep(0.1,3) ) -# initial values +names(inits[[i]]) <- c('mu','p10','alpha') +} +site.cov <- rbind(c(4,1),c(1,1)) +colnames(site.cov) <- c('var_a','var_b') +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)), +site.cov = site.cov), +cov = c('var_a','var_b'), initial_values = inits, +multicore = FALSE), +paste("Initial values for 'p10' should be numeric.", +paste0('See the eDNAjoint guide for help formatting ', +'initial values: '), +paste0('https://bookdown.org/abigailkeller/eDNAjoint_', +'vignette/usecase1.html#initialvalues'), +sep='\n')) +#40. initial values check: p10 length +n.chain <- 4 inits <- list() -inits[[1]] <- list( -mu = mu, -p10 = log_p10, -alpha = alpha +for(i in 1:n.chain){ +inits[[i]] <- list( +mu <- stats::runif(2,0,1), +p10 <- stats::runif(2,log(0.0001),log(0.08)), +alpha <- rep(0.1,3) ) -names(inits[[1]]) <- c('mu','p10','alpha' +names(inits[[i]]) <- c('mu','p10','alpha') +} +site.cov <- rbind(c(4,1),c(1,1)) +colnames(site.cov) <- c('var_a','var_b') +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)), +site.cov = site.cov), +cov = c('var_a','var_b'), initial_values = inits, +multicore = FALSE), +paste("The length of initial values for 'p10' should equal 1.", +paste0('See the eDNAjoint guide for help formatting ', +'initial values: '), +paste0('https://bookdown.org/abigailkeller/eDNAjoint_', +'vignette/usecase1.html#initialvalues'), +sep='\n')) +#41. initial values check: beta numeric +n.chain <- 4 +inits <- list() +for(i in 1:n.chain){ +inits[[i]] <- list( +mu <- stats::runif(2,0,1), +p10 <- stats::runif(1,log(0.0001),log(0.08)), +beta <- "1" ) -# run model -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)}) -# get output params -output_params <- rownames(as.data.frame(jointSummarize(fit$model))) -# test expectation -expect_true(all(c('p10','q[1]','alpha[1]') %in% output_params)) -# check length of of mu -mu_estimates <- as.vector(rstan::summary(fit$model,par='mu')$summary[,1]) -expect_equal(length(mu_estimates),nsite*2) -expect_true(is.numeric(mu_estimates)) -## 9. -# model includes 'p10','alpha','q', 'alpha_gamma', 'alpha_beta' -# constants -nsite <- 20 -nobs_count <- 100 -nobs_pcr <- 8 -# params -mu <- rlnorm(nsite, meanlog = log(1), sdlog = 1) -alpha <- c(0.5, 0.1, -0.4) -log_p10 <- -4.5 -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){ -for(j in 1:nobs_count){ -if(count_type[i,j]==1){ -count[i,j] <- rgamma(1,shape = alpha_gamma[i],rate = beta_gamma) -} else { -count[i,j] <- rgamma(1,shape = alpha_gamma[i]*q,rate = beta_gamma) -} -} -} -# site-level covariates -mat_site <- matrix(NA, nrow = nsite, ncol = length(alpha)) -mat_site[,1] <- 1 # intercept -for(i in 2:length(alpha)){ -mat_site[,i] <- rnorm(nsite,0,1) -} -colnames(mat_site) <- c('int','var_a','var_b') -# 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(sum(mat_site[i,]*alpha))) -p[i] <- min(p11[i] + exp(log_p10),1) -} -# 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)) -} -# add NA to sites 2 and 4 -count[c(2,4),] <- NA -count_type[c(2,4),] <- NA -# collect data -data <- list( -qPCR.N = qPCR.N, -qPCR.K = qPCR.K, -count = count, -count.type = count_type, -site.cov = mat_site +names(inits[[i]]) <- c('mu','p10','beta') +} +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))), +initial_values = inits, +multicore = FALSE), +paste("Initial values for 'beta' should be numeric.", +paste0('See the eDNAjoint guide for help formatting ', +'initial values: '), +paste0('https://bookdown.org/abigailkeller/eDNAjoint_', +'vignette/usecase1.html#initialvalues'), +sep='\n')) +#42. initial values check: beta length +n.chain <- 4 +inits <- list() +for(i in 1:n.chain){ +inits[[i]] <- list( +mu <- stats::runif(2,0,1), +p10 <- stats::runif(1,log(0.0001),log(0.08)), +beta <- c(1,0) ) -# initial values +names(inits[[i]]) <- c('mu','p10','beta') +} +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))), +initial_values = inits, +multicore = FALSE), +paste("The length of initial values for 'beta' should equal 1.", +paste0('See the eDNAjoint guide for help formatting ', +'initial values: '), +paste0('https://bookdown.org/abigailkeller/eDNAjoint_', +'vignette/usecase1.html#initialvalues'), +sep='\n')) +#43. initial values check: alpha numeric +n.chain <- 4 inits <- list() -inits[[1]] <- list( -alpha_gamma = mu, -beta_gamma = rep(1,length(mu)), -p10 = log_p10, -alpha = alpha, -q = q +for(i in 1:n.chain){ +inits[[i]] <- list( +mu <- stats::runif(2,0,1), +p10 <- stats::runif(1,log(0.0001),log(0.08)), +alpha <- c("1","2") ) -names(inits[[1]]) <- c('alpha_gamma','beta_gamma','p10','alpha','q' +names(inits[[i]]) <- c('mu','p10','alpha') +} +site.cov <- rbind(c(4,1),c(1,1)) +colnames(site.cov) <- c('var_a','var_b') +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)), +site.cov = site.cov), +cov = c('var_a','var_b'), initial_values = inits, +multicore = FALSE), +paste("Initial values for 'alpha' should be numeric.", +paste0('See the eDNAjoint guide for help formatting ', +'initial values: '), +paste0('https://bookdown.org/abigailkeller/eDNAjoint_', +'vignette/usecase2.html#initialvalues'), +sep='\n')) +#44. initial values check: alpha length +n.chain <- 4 +inits <- list() +for(i in 1:n.chain){ +inits[[i]] <- list( +mu <- stats::runif(2,0,1), +p10 <- stats::runif(1,log(0.0001),log(0.08)), +alpha <- rep(0.1,2) ) -# run model -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 -)}) -# get output params -output_params <- rownames(as.data.frame(jointSummarize(fit$model))) -# test expectation -expect_true(all(c('p10','q[1]','alpha[1]') %in% output_params)) -# check length of of mu -mu_estimates <- as.vector(rstan::summary(fit$model,par='mu')$summary[,1]) -expect_equal(length(mu_estimates),nsite*2) -expect_true(is.numeric(mu_estimates)) -## 10. -# model includes 'p10','phi','alpha' -# constants -nsite <- 20 -nobs_count <- 100 -nobs_pcr <- 8 -# params -mu <- rlnorm(nsite, meanlog = log(1), sdlog = 1) -alpha <- c(0.5, 0.1, -0.4) -log_p10 <- -4.5 -phi <- 10 -# count -count <- matrix(NA, nrow = nsite, ncol = nobs_count) -for(i in 1:nsite){ -count[i,] <- rnbinom(n = nobs_count, mu = mu[i], size = phi) -} -# site-level covariates -mat_site <- matrix(NA, nrow = nsite, ncol = length(alpha)) -mat_site[,1] <- 1 # intercept -for(i in 2:length(alpha)){ -mat_site[,i] <- rnorm(nsite,0,1) -} -colnames(mat_site) <- c('int','var_a','var_b') -# 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(sum(mat_site[i,]*alpha))) -p[i] <- min(p11[i] + exp(log_p10),1) -} -# 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)) -} -# add NA to sites 2 and 4 -count[c(2,4),] <- NA -# collect data -data <- list( -qPCR.N = qPCR.N, -qPCR.K = qPCR.K, -count = count, -site.cov = mat_site +names(inits[[i]]) <- c('mu','p10','alpha') +} +site.cov <- rbind(c(4,1),c(1,1)) +colnames(site.cov) <- c('var_a','var_b') +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)), +site.cov = site.cov), +cov = c('var_a','var_b'), initial_values = inits, +multicore = FALSE), +paste(paste0("The length of initial values for 'alpha' should ", +"equal\\: \\# covariates \\+ 1 \\(i.e., ", +"including intercept\\)."), +paste0('See the eDNAjoint guide for help formatting ', +'initial values: '), +paste0('https://bookdown.org/abigailkeller/eDNAjoint_', +'vignette/usecase2.html#initialvalues'), +sep='\n')) +#45. initial values check: q numeric +n.chain <- 4 +inits <- list() +for(i in 1:n.chain){ +inits[[i]] <- list( +q <- "0.1" ) -# initial values +names(inits[[i]]) <- c('q') +} +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))), +initial_values = inits, +multicore = FALSE), +paste("Initial values for 'q' should be numeric.", +paste0('See the eDNAjoint guide for help formatting ', +'initial values: '), +paste0('https://bookdown.org/abigailkeller/eDNAjoint_', +'vignette/usecase1.html#initialvalues'), +sep='\n')) +#46. initial values check: q length +n.chain <- 4 inits <- list() -inits[[1]] <- list( -mu = mu, -p10 = log_p10, -alpha = alpha, -phi = phi +for(i in 1:n.chain){ +inits[[i]] <- list( +q <- c(0.1,0.1) ) -names(inits[[1]]) <- c('mu','p10','alpha', -'phi') -# run model -fit <- suppressWarnings({jointModel(data = data, family = 'negbin', -cov = c('var_a','var_b'),n.iter.burn = 25, -n.iter.sample = 75, -n.chain = 1, multicore = FALSE, seed = 10, -initial_values = inits)}) -# get output params -output_params <- rownames(as.data.frame(jointSummarize(fit$model))) -# test expectation -expect_true(all(c('p10','phi','alpha[1]') %in% output_params)) -# check length of of mu -mu_estimates <- as.vector(rstan::summary(fit$model,par='mu')$summary[,1]) -expect_equal(length(mu_estimates),nsite) -expect_true(is.numeric(mu_estimates)) -## 11. -# model includes 'p10','alpha' -# constants -nsite <- 20 -nobs_count <- 100 -nobs_pcr <- 8 -# params -mu <- rlnorm(nsite, meanlog = log(1), sdlog = 1) -alpha <- c(0.5, 0.1, -0.4) -log_p10 <- -4.5 -# count -count <- matrix(NA, nrow = nsite, ncol = nobs_count) -for(i in 1:nsite){ -count[i,] <- rpois(nobs_count, mu[i]) -} -# site-level covariates -mat_site <- matrix(NA, nrow = nsite, ncol = length(alpha)) -mat_site[,1] <- 1 # intercept -for(i in 2:length(alpha)){ -mat_site[,i] <- rnorm(nsite,0,1) -} -colnames(mat_site) <- c('int','var_a','var_b') -# 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(sum(mat_site[i,]*alpha))) -p[i] <- min(p11[i] + exp(log_p10),1) -} -# 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)) -} -# add NA to sites 2 and 4 -count[c(2,4),] <- NA -# collect data -data <- list( -qPCR.N = qPCR.N, -qPCR.K = qPCR.K, -count = count, -site.cov = mat_site +names(inits[[i]]) <- c('q') +} +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))), +initial_values = inits, +multicore = FALSE), +paste(paste0("The length of initial values for 'q' should ", +"equal: \\# unique gear types \\- 1 \\(i.e., q ", +"for reference type = 1\\)."), +paste0('See the eDNAjoint guide for help formatting ', +'initial values: '), +paste0('https://bookdown.org/abigailkeller/eDNAjoint_', +'vignette/usecase1.html#initialvalues'), +sep='\n')) +#47. check length and range of n.chain +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.chain = c(1,1),multicore = FALSE), +paste0("n.chain 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.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 +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.") ) -# initial values -inits <- list() -inits[[1]] <- list( -mu = mu, -p10 = log_p10, -alpha = alpha +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.") ) -names(inits[[1]]) <- c('mu','p10','alpha' +#49. check length and range of n.iter.burn +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.") ) -# run model -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)}) -# get output params -output_params <- rownames(as.data.frame(jointSummarize(fit$model))) -# test expectation -expect_true(all(c('p10','alpha[1]') %in% output_params)) -# check length of of mu -mu_estimates <- as.vector(rstan::summary(fit$model,par='mu')$summary[,1]) -expect_equal(length(mu_estimates),nsite) -expect_true(is.numeric(mu_estimates)) -## 12. -# model includes 'p10','alpha','alpha_gamma','beta_gamma' -# constants -nsite <- 20 -nobs_count <- 100 -nobs_pcr <- 8 -# params -mu <- rlnorm(nsite, meanlog = log(1), sdlog = 1) -alpha <- c(0.5, 0.1, -0.4) -log_p10 <- -4.5 -beta_gamma <- 1 -alpha_gamma <- mu * beta_gamma -# count -count <- matrix(NA, nrow = nsite, ncol = nobs_count) -for(i in 1:nsite){ -count[i,] <- rgamma(nobs_count,shape = alpha_gamma[i],rate = beta_gamma) -} -# site-level covariates -mat_site <- matrix(NA, nrow = nsite, ncol = length(alpha)) -mat_site[,1] <- 1 # intercept -for(i in 2:length(alpha)){ -mat_site[,i] <- rnorm(nsite,0,1) -} -colnames(mat_site) <- c('int','var_a','var_b') -# 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(sum(mat_site[i,]*alpha))) -p[i] <- min(p11[i] + exp(log_p10),1) -} -# 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)) -} -# add NA to sites 2 and 4 -count[c(2,4),] <- NA -# collect data -data <- list( -qPCR.N = qPCR.N, -qPCR.K = qPCR.K, -count = count, -site.cov = mat_site +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.") ) -# initial values -inits <- list() -inits[[1]] <- list( -alpha_gamma = mu, -beta_gamma = rep(1,length(mu)), -p10 = log_p10, -alpha = alpha +#50. check length and range of thin +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))), +thin = c(1,1),multicore = FALSE), +paste0("thin should be an integer > 0 and of length 1.") ) -names(inits[[1]]) <- c('alpha_gamma','beta_gamma','p10','alpha' +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))), +thin = 0,multicore = FALSE), +paste0("thin should be an integer > 0 and of length 1.") ) -# run model -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)}) -# get output params -output_params <- rownames(as.data.frame(jointSummarize(fit$model))) -# test expectation -expect_true(all(c('p10','alpha[1]') %in% output_params)) -# check length of of mu -mu_estimates <- as.vector(rstan::summary(fit$model,par='mu')$summary[,1]) -expect_equal(length(mu_estimates),nsite) -expect_true(is.numeric(mu_estimates)) +#51. check length and range of adapt_delta +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))), +adapt_delta = c(0.9,0.9),multicore = FALSE), +paste0("adapt_delta should be a numeric value > 0 and < 1 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))), +adapt_delta = 1.2,multicore = FALSE), +paste0("adapt_delta should be a numeric value > 0 and < 1 and ", +"of length 1.") +) +#52. check length of seed +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))), +seed = c(1,2),multicore = FALSE), +paste0("seed should be an integer of length 1.") +) +#53. check K <= N +expect_error(jointModel(data = list(qPCR.K = rbind(c(4,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)), +multicore = FALSE)), +paste(paste0("N should be >= K in qPCR data. N is the number ", +"of qPCR replicates per sample, and K is the ", +"number of positive detections among replicates."), +'See the eDNAjoint guide for data formatting help: ', +paste0('https://bookdown.org/abigailkeller/eDNAjoint_', +'vignette/usecase1.html#prepare-the-data'), +sep='\n') +) +}) devtools::install() +library(devtools) +install_github('abigailkeller/eDNAjoint') diff --git a/README.Rmd b/README.Rmd index 0854565..dfdc343 100644 --- a/README.Rmd +++ b/README.Rmd @@ -27,7 +27,7 @@ knitr::opts_chunk$set( [![Status at rOpenSci Software Peer Review](https://badges.ropensci.org/642_status.svg)](https://github.com/ropensci/software-review/issues/642) [![Lifecycle: -experimental](https://img.shields.io/badge/lifecycle-experimental-orange.svg)](https://lifecycle.r-lib.org/articles/stages.html#experimental) +stable](https://img.shields.io/badge/lifecycle-stable-green.svg)](https://lifecycle.r-lib.org/articles/stages.html#stable) diff --git a/README.md b/README.md index 83420c6..ec851c3 100644 --- a/README.md +++ b/README.md @@ -12,7 +12,7 @@ [![Status at rOpenSci Software Peer Review](https://badges.ropensci.org/642_status.svg)](https://github.com/ropensci/software-review/issues/642) [![Lifecycle: -experimental](https://img.shields.io/badge/lifecycle-experimental-orange.svg)](https://lifecycle.r-lib.org/articles/stages.html#experimental) +stable](https://img.shields.io/badge/lifecycle-stable-green.svg)](https://lifecycle.r-lib.org/articles/stages.html#stable) @@ -73,8 +73,8 @@ detection, $p_{10}$: ``` r # summarize p10 posterior jointSummarize(goby.fit$model, par = 'p10') -#> mean se_mean sd 2.5% 97.5% n_eff Rhat -#> p10 0.003 0 0.001 0.001 0.007 17646.07 1 +#> mean se_mean sd 2.5% 97.5% n_eff Rhat +#> p10 0.003 0 0.001 0.001 0.007 15361.8 1 ``` Or to find the number of eDNA samples and traditional survey samples