test_that("jointModel input checks work", {
testthat::skip_on_cran()
#' @srrstats {G5.2,G5.2b,BS2.15} Tests the assure function input checks are
#' behaving as expected.
#1. input tags are valid, q = FALSE, cov = FALSE
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))),
multicore = FALSE),
paste("Data should include 'qPCR.N', 'qPCR.K', and 'count'.",
'See the eDNAjoint guide for data formatting help: ',
paste0('https://bookdown.org/abigailkeller/eDNAjoint_',
'vignette/usecase1.html#prepare-the-data'),
sep='\n'))
#2. input tags are valid, q = FALSE, cov = TRUE
site.cov <- cbind(c(1,0),c(0.4,-0.4))
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(paste0("Data should include 'qPCR.N', 'qPCR.K', ",
"'count', and 'site.cov'."),
'See the eDNAjoint guide for data formatting help: ',
paste0('https://bookdown.org/abigailkeller/eDNAjoint_',
'vignette/usecase2.html#prepare-the-data'),
sep='\n'))
#3. input tags are valid, q = TRUE, cov = FALSE
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))),
q = TRUE,
multicore = FALSE),
paste(paste0("Data should include 'qPCR.N', 'qPCR.K', ",
"'count', and 'count.type'."),
'See the eDNAjoint guide for data formatting help: ',
paste0('https://bookdown.org/abigailkeller/eDNAjoint_',
'vignette/usecase3.html#prepare-the-data'),
sep='\n'))
#4. input tags are valid, q = TRUE, cov = TRUE
site.cov <- cbind(c(1,0),c(0.4,-0.4))
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("Data should include 'qPCR.N', 'qPCR.K', 'count', ",
"'count.type', and 'site.cov'."),
'See the eDNAjoint guide for data formatting help: ',
paste0('https://bookdown.org/abigailkeller/eDNAjoint_',
'vignette/usecase3.html#prepare-the-data'),
sep='\n'))
#5. make sure dimensions of qPCR.N and qPCR.K are equal
#' @srrstats {BS2.1a} Test to ensure pre-processing routines to ensure all
#' input data is dimensionally commensurate
expect_error(jointModel(data = list(qPCR.K = rbind(c(1,1,1,1),c(1,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("Dimensions of qPCR.N and qPCR.K do not match.",
'See the eDNAjoint guide for data formatting help: ',
paste0('https://bookdown.org/abigailkeller/eDNAjoint_',
'vignette/usecase1.html#prepare-the-data'),
sep='\n'))
#6. make sure dimensions of count and count.type are equal, if count.type is
# present
#' @srrstats {BS2.1a} Test to ensure pre-processing routines to ensure all
#' input data is dimensionally commensurate
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),c(1,2))),
q = TRUE,
multicore = FALSE),
paste("Dimensions of count and count.type do not match.",
'See the eDNAjoint guide for data formatting help: ',
paste0('https://bookdown.org/abigailkeller/eDNAjoint_',
'vignette/usecase3.html#prepare-the-data'),
sep='\n'))
#7. make sure number of rows in count = number of rows in qPCR.N and qPCR.K
#' @srrstats {BS2.1a} Test to ensure pre-processing routines to ensure all
#' input data is dimensionally commensurate
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),
c(4,1,1))),
multicore = FALSE),
paste(paste0("Number of sites \\(rows\\) in qPCR data and ",
"traditional survey count data do not match."),
'See the eDNAjoint guide for data formatting help: ',
paste0('https://bookdown.org/abigailkeller/eDNAjoint_',
'vignette/usecase1.html#prepare-the-data'),
sep='\n'))
#8. make sure all data is numeric -- if q == TRUE
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('NA',2,2),c(1,2,2))),
q = TRUE,
multicore = FALSE),
paste("Data should be numeric.",
'See the eDNAjoint guide for data formatting help: ',
paste0('https://bookdown.org/abigailkeller/eDNAjoint_',
'vignette/usecase3.html#prepare-the-data'),
sep='\n'))
#9. make sure all data is numeric -- if q == FALSE
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))),
multicore = FALSE),
paste("Data should be numeric.",
'See the eDNAjoint guide for data formatting help: ',
paste0('https://bookdown.org/abigailkeller/eDNAjoint_',
'vignette/usecase1.html#prepare-the-data'),
sep='\n'))
#10. make sure locations of NAs in count data match locations of NAs in
# count.type data
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(NA,2,2),c(1,2,2))),
q = TRUE,
multicore = FALSE),
paste(paste0("Empty data cells \\(NA\\) in count data should ",
"match ",
"empty data cells \\(NA\\) in count.type data."),
'See the eDNAjoint guide for data formatting help: ',
paste0('https://bookdown.org/abigailkeller/eDNAjoint_',
'vignette/usecase3.html#prepare-the-data'),
sep='\n'))
#11. make sure locations of NAs in qPCR.N data match locations of NAs in
# qPCR.K data
expect_error(jointModel(data = list(qPCR.K = rbind(c(1,1,1),c(1,NA,1)),
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("Empty data cells \\(NA\\) in qPCR.N data should ",
"match ",
"empty data cells \\(NA\\) in qPCR.K data."),
'See the eDNAjoint guide for data formatting help: ',
paste0('https://bookdown.org/abigailkeller/eDNAjoint_',
'vignette/usecase1.html#prepare-the-data'),
sep='\n'))
#12. make sure family is either 'poisson', 'negbin', or 'gamma'
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))),
family = 'normal',
multicore = FALSE),
paste0("Invalid family. Options include 'poisson', ",
"'negbin', and 'gamma'."))
#13. p10 priors is a vector of two numeric values
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))),
p10priors = c(1,1,2),
multicore = FALSE),
paste0("p10priors should be a vector of two positive numeric ",
"values. ex. c\\(1,20\\)"))
#14. p10 priors is a vector of two numeric values
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))),
multicore = FALSE,
p10priors = '1,20'),
paste0("p10priors should be a vector of two positive numeric ",
"values. ex. c\\(1,20\\)"))
#15. p10 priors is a vector of two numeric values
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))),
p10priors = c(0,20),
multicore = FALSE),
paste0("p10priors should be a vector of two positive numeric ",
"values. ex. c\\(1,20\\)"))
#16. phi priors is a vector of two numeric values
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))),
phipriors = c(0,1),family = 'negbin',
multicore = FALSE),
paste0("phipriors should be a vector of two positive numeric ",
"values. ex. c\\(0.25,0.25\\)"))
#17. the smallest count.type is 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(0,1,2),c(1,2,NA))),
q = TRUE,
multicore = FALSE),
paste(paste0("The first gear type should be referenced as 1 in ",
"count.type. Subsequent gear types should be ",
"referenced 2, 3, 4, etc."),
'See the eDNAjoint guide for data formatting help: ',
paste0('https://bookdown.org/abigailkeller/eDNAjoint_',
'vignette/usecase3.html#prepare-the-data'),
sep='\n'))
#18. count are integers
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,1),c(1,1,NA)),
count.type = rbind(c(1,1,2),c(1,2,NA))),
q = TRUE, family = 'negbin',
multicore = FALSE),
paste0("All values in count should be non-negative integers. ",
"Use family = 'gamma' if count is continuous."))
#19. qPCR.N are integers
expect_error(jointModel(data = list(qPCR.K = rbind(c(0.99,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,1,2),c(1,2,NA))),
q = TRUE,
multicore = FALSE),
paste("All values in qPCR.K should be non-negative integers.",
'See the eDNAjoint guide for data formatting help: ',
paste0('https://bookdown.org/abigailkeller/eDNAjoint_',
'vignette/usecase1.html#prepare-the-data'),
sep='\n'))
#20. qPCR.K are integers
expect_error(jointModel(data = list(qPCR.K = rbind(c(1,1,1),c(1,1,NA)),
qPCR.N = rbind(c(3.1,3,3),c(3,3,NA)),
count = rbind(c(4,1,1),c(1,1,NA)),
count.type = rbind(c(1,1,2),c(1,2,NA))),
q = TRUE,
multicore = FALSE),
paste("All values in qPCR.N should be non-negative integers.",
'See the eDNAjoint guide for data formatting help: ',
paste0('https://bookdown.org/abigailkeller/eDNAjoint_',
'vignette/usecase1.html#prepare-the-data'),
sep='\n'))
#21. count.type are integers
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.1,1,2),c(1,2,NA))),
q = TRUE,
multicore = FALSE),
paste("All values in count.type should be integers.",
'See the eDNAjoint guide for data formatting help: ',
paste0('https://bookdown.org/abigailkeller/eDNAjoint_',
'vignette/usecase3.html#prepare-the-data'),
sep='\n'))
#22. site.cov is numeric, if present
site.cov <- cbind(c('high','low'),c(0.4,-0.4))
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("site.cov should be numeric.",
'See the eDNAjoint guide for data formatting help: ',
paste0('https://bookdown.org/abigailkeller/eDNAjoint_',
'vignette/usecase2.html'),
sep='\n'))
#23. cov values match column names in site.cov
site.cov <- cbind(c(0,1),c(0.4,-0.4))
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("cov values should be listed in the column names of ",
"site.cov in the data."),
'See the eDNAjoint guide for data formatting help: ',
paste0('https://bookdown.org/abigailkeller/eDNAjoint_',
'vignette/usecase2.html'),
sep='\n'))
#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()
for(i in 1:n.chain){
inits[[i]] <- list(
mu <- stats::runif(3, 0.01, 5),
p10 <- exp(stats::runif(1,log(0.0001),log(0.08))),
alpha <- rep(0.1,3)
)
}
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()
for(i in 1:n.chain){
inits[[i]] <- list(
mu <- stats::runif(3, -1, 0),
p10 <- exp(stats::runif(1,log(0.0001),log(0.08))),
alpha <- rep(0.1,3)
)
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 <- exp(stats::runif(1,log(0.0001),log(0.08))),
alpha <- rep(0.1,3)
)
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)
)
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()
for(i in 1:n.chain){
inits[[i]] <- list(
mu <- stats::runif(2,0,1),
p10 <- exp(stats::runif(2,log(0.0001),log(0.08))),
alpha <- rep(0.1,3)
)
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 <- exp(stats::runif(1,log(0.0001),log(0.08))),
beta <- "1"
)
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 <- exp(stats::runif(1,log(0.0001),log(0.08))),
beta <- c(1,0)
)
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()
for(i in 1:n.chain){
inits[[i]] <- list(
mu <- stats::runif(2,0,1),
p10 <- exp(stats::runif(1,log(0.0001),log(0.08))),
alpha <- c("1","2")
)
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 <- exp(stats::runif(1,log(0.0001),log(0.08))),
alpha <- rep(0.1,2)
)
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"
)
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()
for(i in 1:n.chain){
inits[[i]] <- list(
q <- c(0.1,0.1)
)
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.")
)
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.")
)
#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.")
)
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.")
)
#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.")
)
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.")
)
#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')
)
})
# correctness and parameter recovery tests
#' @srrstats {G5.4, G5.6} Correctness/parameter recovery tests to test that
#' the implementation produces expected results given data with known
#' properties
#' @srrstats {PD4.0} These tests do not test for numeric equality of outputs,
#' but rather test for the recovery of parameter values. In these tests, I
#' simulate data with known parameter values, run the main statistical
#' function (jointModel()) with the data, and then test if the known parameter
#' values are within the 95% credibility interval of the parameters'
#' posteriors in the function output.
test_that("jointModel parameter recovery tests work",{
testthat::skip_on_cran()
################################
# model run 1: smaller dataset #
################################
# simulate data: seed 123
#' @srrstats {G5.5, G5.6} Running correctness test with a fixed random seed
set.seed(123)
# constants
nsite <- 20
nobs_count <- 200
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))
}
# collect data
data <- list(
qPCR.N = qPCR.N,
qPCR.K = qPCR.K,
count = count,
site.cov = mat_site
)
# run model
fit1 <- suppressWarnings({jointModel(data = data, cov = c('var_a','var_b'),
multicore = FALSE, seed = 10)})
# summary
summary1 <- as.data.frame(rstan::summary(fit1$model,
pars = c('mu','alpha','log_p10'),
probs = c(0.025, 0.975))$summary)
# set up empty vector to check if true values are in 95% interval of
# posterior estimates
check <- rep(NA,length(mu)+length(alpha)+length(log_p10))
# check mu
for(i in 1:nsite){
par <- paste0('mu[',i,']')
if(mu[i] > summary1[par,'2.5%'] && mu[i] < summary1[par,'97.5%']){
check[i] <- TRUE
} else {
check[i] <- FALSE
}
}
# check alpha
for(i in seq_along(alpha)){
par <- paste0('alpha[',i,']')
if(alpha[i] > summary1[par,'2.5%'] && alpha[i] < summary1[par,'97.5%']){
check[nsite+i] <- TRUE
} else {
check[nsite+i] <- FALSE
}
}
# check p10
if(log_p10 > summary1['log_p10','2.5%'] &&
log_p10 < summary1['log_p10','97.5%']){
check[nsite+length(alpha)+1] <- TRUE
} else {
check[nsite+length(alpha)+1] <- FALSE
}
#' @srrstats {G3.0, G5.6a} Instead of comparing floating point values for
#' equality, here the model is tested to determine if the true parameter
#' values are within the 95% quantiles of the posterior
# all should be equal to true
expect_equal(check,rep(TRUE,nsite+length(alpha)+length(log_p10)))
# test that output values are on the same scale as the data
mu_estimates <- rep(NA,nsite)
for(i in 1:nsite){
mu_estimates[i] <- summary1[paste0('mu[',i,']'),'mean']
}
# get mean of input count data at each site
data_mean <- rep(NA,nsite)
for(i in 1:nsite){
data_mean[i] <- mean(data$count[i,])
}
#' @srrstats {BS7.4, BS7.4a} Check to ensure that mean posterior estimates
#' are on the same scale as the mean of the input data (here checking
#' estimates of mu, i.e., expected catch rate at each site)
# check that estimates are on same scale as data
expect_equal(round(mu_estimates,0),round(data_mean,0))
##########################################
# model run 2: smaller dataset, new seed #
##########################################
#' @srrstats {G5.6b,G5.9b} Run test for multiple seeds with same data
# simulate data: seed 222
set.seed(222)
# constants
nsite <- 20
nobs_count <- 200
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))
}
# collect data
data <- list(
qPCR.N = qPCR.N,
qPCR.K = qPCR.K,
count = count,
site.cov = mat_site
)
# run model
fit2 <- suppressWarnings({jointModel(data = data, cov = c('var_a','var_b'),
multicore = FALSE, seed = 2)})
# summary
summary2 <- as.data.frame(rstan::summary(fit2$model,
pars = c('mu','alpha','log_p10'),
probs = c(0.025, 0.975))$summary)
# set up empty vector to check if true values are in 95% interval of
# posterior estimates
check <- rep(NA,length(mu)+length(alpha)+length(log_p10))
# check mu
for(i in 1:nsite){
par <- paste0('mu[',i,']')
if(mu[i] > summary2[par,'2.5%'] && mu[i] < summary2[par,'97.5%']){
check[i] <- TRUE
} else {
check[i] <- FALSE
}
}
# check alpha
for(i in seq_along(alpha)){
par <- paste0('alpha[',i,']')
if(alpha[i] > summary2[par,'2.5%'] && alpha[i] < summary2[par,'97.5%']){
check[nsite+i] <- TRUE
} else {
check[nsite+i] <- FALSE
}
}
# check p10
if(log_p10 > summary2['log_p10','2.5%'] &&
log_p10 < summary2['log_p10','97.5%']){
check[nsite+length(alpha)+1] <- TRUE
} else {
check[nsite+length(alpha)+1] <- FALSE
}
#' @srrstats {G3.0, G5.6a} Instead of comparing floating point values for
#' equality, here the model is tested to determine if the true parameter
#' values are within the 95% quantiles of the posterior
# all should be equal to true
expect_equal(check,rep(TRUE,nsite+length(alpha)+length(log_p10)))
################################
# model run 3: larger dataset #
################################
# simulate data: seed 123
set.seed(123)
# constants
nsite <- 20
nobs_count <- 200 # increased dataset size (doubled)
nobs_pcr <- 16 # increased dataset size (doubled)
# 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))
}
# collect data
data <- list(
qPCR.N = qPCR.N,
qPCR.K = qPCR.K,
count = count,
site.cov = mat_site
)
# run model
fit_large <- suppressWarnings({jointModel(data = data,
cov = c('var_a','var_b'),
multicore = FALSE, seed = 10)})
# summary
summary_large <- as.data.frame(rstan::summary(fit_large$model,
pars = 'alpha',
probs = c(0.025,
0.975))$summary)
# compare standard errors of estimates of alpha coefficients with
# implementation with a smaller dataset
se_large <- c(summary_large['alpha[2]','se_mean'],
summary_large['alpha[3]','se_mean']
)
# get values for smaller dataset
se_small <- c(summary1['alpha[2]','se_mean'],
summary1['alpha[3]','se_mean']
)
# set up empty vector to check if standard errors are smaller with
# larger dataset
check <- rep(NA,length(alpha)-1)
for(i in seq_along(check)){
check[i] <- se_large[i] < se_small[i]
}
#' @srrstats {G5.7} Check to see that implementation performs as expected
#' properties of data change (i.e., standard error of posteriors is smaller
#' if there are more data observations)
# all should be equal to true
expect_equal(check,rep(TRUE,length(alpha)-1))
})
test_that("jointModel probability distribution tests work",{
testthat::skip_on_cran()
#' @srrstats {PD4.2} This package fits models with fixed distributions to
#' data. Users can choose the distribution used to represent the
#' data generating process for traditional survey data through an input
#' argument to the function (family). This test assesses whether estimates
#' of mu, or the site-level expected catch rate, is different when a
#' negative binomial or poisson distribution is used to represent the data
#' generating process.
#######################################################
# generate data with a negative binomial distribution #
set.seed(222)
# constants
nsite <- 20
nobs_count <- 200
nobs_pcr <- 8
# params
mu <- rlnorm(nsite,meanlog = log(1),sdlog = 1)
beta <- 0.5
log_p10 <- -4.5
phi <- 0.1
# count
count <- matrix(NA,nrow = nsite,ncol = nobs_count)
for(i in 1:nsite){
count[i,] <- rnbinom(nobs_count,mu = mu[i],size = phi)
}
# 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)
}
# 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
)
##############################################################
# fit data with a negative binomial and poisson distribution #
# run model -- negative binomial distribution
negbin.fit <- suppressWarnings({jointModel(data = data, multicore = FALSE,
family = 'negbin',seed = 2)})
# run model -- poisson distribution
pois.fit <- suppressWarnings({jointModel(data = data, multicore = FALSE,
family = 'poisson',seed = 2)})
# summarize outputs
negbin_summary <- as.data.frame(rstan::summary(negbin.fit$model,
pars = 'mu',
probs = c(0.025,
0.975))$summary)
pois_summary <- as.data.frame(rstan::summary(pois.fit$model,
pars = 'mu',
probs = c(0.025,
0.975))$summary)
# summarize differences in mean estimates of mu
summary_table <- table(negbin_summary$mean > pois_summary$mean)
# calculate percentage of mean estimates of mu that are greater when
# a negative binomial distribution are used rather than a poisson distribution
percent <- summary_table[[2]]/(summary_table[[2]]+summary_table[[1]])
# test that estimates of mu when a negative binomial distribution is used is
# greater than estimates of mu when a poisson distribution is used for
# >80% of sites
expect_gt(percent,0.8)
})
test_that("semi-paired model works", {
testthat::skip_on_cran()
## 1.
# model includes 'p10','q','phi','beta'
# 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
phi <- 1.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)
} else {
count[i,j] <- rnbinom(n = 1, mu = mu[i]*q, size = phi)
}
}
}
# 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)
}
# 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
)
# initial values
inits <- list()
inits[[1]] <- list(
mu = mu,
p10 = exp(log_p10),
beta = beta,
phi = phi,
q = q
)
names(inits[[1]]) <- c('mu','p10','beta','phi','q')
# run model
fit <- suppressWarnings({jointModel(data = data, q = TRUE, family = 'negbin',
initial_values = inits, n.iter.burn = 25,
n.iter.sample = 75,
n.chain = 1, multicore = FALSE, seed = 10
)})
# get output params
output_params <- rownames(as.data.frame(jointSummarize(fit$model)))
# test expectation
expect_true(all(c('p10','q[1]','phi','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*2)
expect_true(is.numeric(mu_estimates))
## 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])
} else {
count[i,j] <- rpois(1, mu[i]*q)
}
}
}
# 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)
}
# 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
)
# initial values
inits <- list()
inits[[1]] <- list(
mu = mu,
p10 = exp(log_p10),
beta = beta,
q = q
)
names(inits[[1]]) <- c('mu','p10','beta','q')
# 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)})
# get output params
output_params <- rownames(as.data.frame(jointSummarize(fit$model)))
# test expectation
expect_true(all(c('p10','q[1]','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*2)
expect_true(is.numeric(mu_estimates))
## 3.
# model includes 'p10','beta','q', 'alpha_gamma','beta_gamma'
# 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
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)
}
}
}
# 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)
}
# 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
)
# initial values
inits <- list()
inits[[1]] <- list(
alpha_gamma = mu,
beta_gamma = rep(1,length(mu)),
p10 = exp(log_p10),
beta = beta,
q = q
)
names(inits[[1]]) <- c('alpha_gamma','beta_gamma','p10','beta','q')
# 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)})
# get output params
output_params <- rownames(as.data.frame(jointSummarize(fit$model)))
# test expectation
expect_true(all(c('p10','q[1]','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*2)
expect_true(is.numeric(mu_estimates))
## 4.
# model includes 'p10','phi','beta'
# constants
nsite <- 50
nobs_count <- 100
nobs_pcr <- 8
# params
mu <- rlnorm(nsite, meanlog = log(1), sdlog = 1)
beta <- 0.5
log_p10 <- -4.5
phi <- 1.2
# 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)
}
# 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)
}
# 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
)
# add NA to sites 2 and 4
count[c(2,4),] <- NA
# initial values
inits <- list()
inits[[1]] <- list(
mu = mu,
p10 = exp(log_p10),
beta = beta,
phi = phi
)
names(inits[[1]]) <- c('mu','p10','beta','phi')
# run model
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)})
# get output params
output_params <- rownames(as.data.frame(jointSummarize(fit$model)))
# test expectation
expect_true(all(c('p10','phi','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))
## 5.
# model includes 'p10','beta'
# constants
nsite <- 50
nobs_count <- 100
nobs_pcr <- 8
# params
mu <- rlnorm(nsite, meanlog = log(1), sdlog = 1)
beta <- 0.5
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])
}
# 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)
}
# 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
)
# initial values
inits <- list()
inits[[1]] <- list(
mu = mu,
p10 = exp(log_p10),
beta = beta
)
names(inits[[1]]) <- c('mu','p10','beta')
# 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)})
# 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))
## 6.
# model includes 'p10','beta','alpha_gamma','alpha_beta'
# 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
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)
}
# 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)
}
# 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
)
# initial values
inits <- list()
inits[[1]] <- list(
alpha_gamma = mu,
beta_gamma = rep(1,length(mu)),
p10 = exp(log_p10),
beta = beta
)
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
inits <- list()
inits[[1]] <- list(
mu = mu,
p10 = exp(log_p10),
alpha = alpha,
q = q,
phi = phi
)
names(inits[[1]]) <- c('mu','p10','alpha','q','phi'
)
# 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
)
# initial values
inits <- list()
inits[[1]] <- list(
mu = mu,
p10 = exp(log_p10),
alpha = alpha
)
names(inits[[1]]) <- c('mu','p10','alpha'
)
# 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
)
# initial values
inits <- list()
inits[[1]] <- list(
alpha_gamma = mu,
beta_gamma = rep(1,length(mu)),
p10 = exp(log_p10),
alpha = alpha,
q = q
)
names(inits[[1]]) <- c('alpha_gamma','beta_gamma','p10','alpha','q'
)
# 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
)
# initial values
inits <- list()
inits[[1]] <- list(
mu = mu,
p10 = exp(log_p10),
alpha = alpha,
phi = phi
)
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
)
# initial values
inits <- list()
inits[[1]] <- list(
mu = mu,
p10 = exp(log_p10),
alpha = alpha
)
names(inits[[1]]) <- c('mu','p10','alpha'
)
# 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
)
# initial values
inits <- list()
inits[[1]] <- list(
alpha_gamma = mu,
beta_gamma = rep(1,length(mu)),
p10 = exp(log_p10),
alpha = alpha
)
names(inits[[1]]) <- c('alpha_gamma','beta_gamma','p10','alpha'
)
# 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))
})