diff --git a/.Rhistory b/.Rhistory
index 2cd5e49..f71c811 100644
--- a/.Rhistory
+++ b/.Rhistory
@@ -1,15 +1,468 @@
-nobs_count <- 100
+q=TRUE,
+multicore=FALSE),
+paste0("The first gear type should be referenced as 1 in ",
+"count.type. Subsequent gear types should be ",
+"referenced 2, 3, 4, etc."))
+#18. count are integers
+expect_error(jointModel(data=list(qPCR.N=rbind(c(1,1,1),c(1,1,NA)),
+qPCR.K=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.N=rbind(c(0.99,1,1),c(1,1,NA)),
+qPCR.K=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),
+"All values in qPCR.N should be non-negative integers.")
+#20. qPCR.K are integers
+expect_error(jointModel(data=list(qPCR.N=rbind(c(1,1,1),c(1,1,NA)),
+qPCR.K=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),
+"All values in qPCR.K should be non-negative integers.")
+#21. count.type are integers
+expect_error(jointModel(data=list(qPCR.N=rbind(c(1,1,1),c(1,1,NA)),
+qPCR.K=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),
+"All values in count.type should be integers.")
+#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.N=rbind(c(1,1,1),c(1,1,NA)),
+qPCR.K=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),
+"site.cov should be numeric.")
+#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.N=rbind(c(1,1,1),c(1,1,NA)),
+qPCR.K=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),
+paste0("cov values should be listed in the column names of ",
+"site.cov in the data."))
+#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.N=rbind(c(1,1,1),c(1,1,NA)),
+qPCR.K=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),
+paste0("The number of rows in site.cov matrix should match the ",
+"number of rows in all other matrices."))
+#25. make sure count.type is not zero-length
+expect_error(jointModel(data=list(qPCR.N=rbind(c(1,1,1),c(1,1,NA)),
+qPCR.K=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),
+"count.type contains zero-length data.")
+#26. make sure no column is entirely NA in count.type
+expect_error(jointModel(data=list(qPCR.N=rbind(c(1,1,1),c(1,1,NA)),
+qPCR.K=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),
+"count.type contains a column with all NA.")
+#27. make sure no column is entirely NA in qPCR.N
+expect_error(jointModel(data=list(qPCR.N=rbind(c(1,1,NA),c(1,1,NA)),
+qPCR.K=rbind(c(3,3,3),c(3,3,NA)),
+count=rbind(c(4,1,1),c(1,1,NA))),
+multicore=FALSE),
+"qPCR.N contains a column with all NA.")
+#28. make sure no column is entirely NA in qPCR.K
+expect_error(jointModel(data=list(qPCR.N=rbind(c(1,1,1),c(1,1,NA)),
+qPCR.K=rbind(c(3,3,NA),c(3,3,NA)),
+count=rbind(c(4,1,1),c(1,1,NA))),
+multicore=FALSE),
+"qPCR.K contains a column with all NA.")
+#29. make sure no column is entirely NA in count
+expect_error(jointModel(data=list(qPCR.N=rbind(c(1,1,1),c(1,1,NA)),
+qPCR.K=rbind(c(3,3,3),c(3,3,NA)),
+count=rbind(c(4,1,NA),c(1,1,NA))),
+multicore=FALSE),
+"count contains a column with all NA.")
+#30. make sure no data are undefined
+expect_error(jointModel(data=list(qPCR.N=rbind(c(1,1,1),c(1,1,NA)),
+qPCR.K=rbind(c(3,3,3),c(3,3,NA)),
+count=rbind(c(4,1,Inf),c(1,1,NA))),
+multicore=FALSE),
+"count contains undefined values \\(i.e., Inf or -Inf\\)")
+#31. make sure no data are undefined
+expect_error(jointModel(data=list(qPCR.N=rbind(c(1,1,1),c(1,1,NA)),
+qPCR.K=rbind(c(3,3,Inf),c(3,3,NA)),
+count=rbind(c(4,1,1),c(1,1,NA))),
+multicore=FALSE),
+"qPCR.K contains undefined values \\(i.e., Inf or -Inf\\)")
+#32. make sure no data are undefined
+expect_error(jointModel(data=list(qPCR.N=rbind(c(1,1,Inf),c(1,1,NA)),
+qPCR.K=rbind(c(3,3,3),c(3,3,NA)),
+count=rbind(c(4,1,1),c(1,1,NA))),
+multicore=FALSE),
+"qPCR.N contains undefined values \\(i.e., Inf or -Inf\\)")
+#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.N=rbind(c(1,1,1),c(1,1,NA)),
+qPCR.K=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),
+"site.cov contains zero-length data.")
+#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.N=rbind(c(1,1,1),c(1,1,NA)),
+qPCR.K=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),
+"site.cov contains a column with all NA.")
+#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.N=rbind(c(1,1,1),c(1,1,NA)),
+qPCR.K=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),
+"site.cov contains undefined values \\(i.e., Inf or -Inf\\)")
+#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 <- 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.N=rbind(c(1,1,1),c(1,1,NA)),
+qPCR.K=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),
+paste0("The length of the list of initial values should equal ",
+"the number of chains \\(n.chain, default is 4\\)."))
+#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 <- 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.N=rbind(c(1,1,1),c(1,1,NA)),
+qPCR.K=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),
+"Initial values for 'mu' should be numeric values > 0.")
+#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)
+)
+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.N=rbind(c(1,1,1),c(1,1,NA)),
+qPCR.K=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),
+paste0("The length of initial values for 'mu' should ",
+"equal the number of sites."))
+#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.N=rbind(c(1,1,1),c(1,1,NA)),
+qPCR.K=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),
+"Initial values for 'p10' should be numeric.")
+#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 <- 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.N=rbind(c(1,1,1),c(1,1,NA)),
+qPCR.K=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),
+"The length of initial values for 'p10' should equal 1.")
+#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"
+)
+names(inits[[i]]) <- c('mu','p10','beta')
+}
+expect_error(jointModel(data=list(qPCR.N=rbind(c(1,1,1),c(1,1,NA)),
+qPCR.K=rbind(c(3,3,3),c(3,3,NA)),
+count=rbind(c(4,1,1),c(1,1,NA))),
+initial_values=inits,
+multicore=FALSE),
+"Initial values for 'beta' should be numeric.")
+#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)
+)
+names(inits[[i]]) <- c('mu','p10','beta')
+}
+expect_error(jointModel(data=list(qPCR.N=rbind(c(1,1,1),c(1,1,NA)),
+qPCR.K=rbind(c(3,3,3),c(3,3,NA)),
+count=rbind(c(4,1,1),c(1,1,NA))),
+initial_values=inits,
+multicore=FALSE),
+"The length of initial values for 'beta' should equal 1.")
+#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 <- 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.N=rbind(c(1,1,1),c(1,1,NA)),
+qPCR.K=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),
+paste0("Initial values for 'alpha' should be numeric."))
+#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)
+)
+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.N=rbind(c(1,1,1),c(1,1,NA)),
+qPCR.K=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),
+paste0("The length of initial values for 'alpha' should ",
+"equal\\: \\# covariates \\+ 1 \\(i.e., ",
+"including intercept\\)."))
+#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.N=rbind(c(1,1,1),c(1,1,NA)),
+qPCR.K=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),
+paste0("Initial values for 'q' should be numeric."))
+#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.N=rbind(c(1,1,1),c(1,1,NA)),
+qPCR.K=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),
+paste0("The length of initial values for 'q' should equal: ",
+"\\# unique gear types \\- 1 \\(i.e., q for reference ",
+"type = 1\\)."))
+#47. check length and range of n.chain
+expect_error(jointModel(data=list(qPCR.N=rbind(c(1,1,1),c(1,1,NA)),
+qPCR.K=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.N=rbind(c(1,1,1),c(1,1,NA)),
+qPCR.K=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.N=rbind(c(1,1,1),c(1,1,NA)),
+qPCR.K=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.N=rbind(c(1,1,1),c(1,1,NA)),
+qPCR.K=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.N=rbind(c(1,1,1),c(1,1,NA)),
+qPCR.K=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.N=rbind(c(1,1,1),c(1,1,NA)),
+qPCR.K=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.N=rbind(c(1,1,1),c(1,1,NA)),
+qPCR.K=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.N=rbind(c(1,1,1),c(1,1,NA)),
+qPCR.K=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.N=rbind(c(1,1,1),c(1,1,NA)),
+qPCR.K=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.N=rbind(c(1,1,1),c(1,1,NA)),
+qPCR.K=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.N=rbind(c(1,1,1),c(1,1,NA)),
+qPCR.K=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.")
+)
+})
+# 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
-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)
+count[i,] <- rpois(nobs_count,mu[i])
 }
 # site-level covariates
 mat_site <- matrix(NA,nrow=nsite,ncol=length(alpha))
@@ -18,8 +471,8 @@ 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 (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){
@@ -43,470 +496,17 @@ 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 = log_p10,
-alpha = alpha
-)
-names(inits[[1]]) <- c('alpha_gamma','beta_gamma','p10','alpha'
-)
 # run model
-fit <- jointModel(data=data,family='gamma',
-cov=c('var_a','var_b'),
-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','alpha[1]') %in% output_params))
-# constants
-nsite <- 20
-nobs_count <- 100
-# params
-mu <- rlnorm(nsite,meanlog=log(1),sdlog=1)
-phi <- 1.2
-q <- 2
-# traditional type
-count_type <- cbind(matrix(1,nrow=nsite,ncol=nobs_count/2),
-matrix(2,nrow=nsite,ncol=nobs_count/2))
-# count
-count <- matrix(NA,nrow=nsite,ncol=nobs_count)
-for(i in 1:nsite){
-for(j in 1:nobs_count){
-if(count_type[i,j]==1){
-count[i,j] <- rnbinom(n=1,mu=mu[i],size=phi)
-} else {
-count[i,j] <- rnbinom(n=1,mu=mu[i]*q,size=phi)
-}
-}
-}
-# collect data
-data <- list(
-count = count,
-count.type = count_type
-)
-# initial values
-inits <- list()
-inits[[1]] <- list(
-mu = mu,
-phi = phi
-)
-names(inits[[1]]) <- c('mu','phi')
-# run model
-fit <- traditionalModel(data=data, q=TRUE, family='negbin',
-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('q[1]','phi') %in% output_params))
-# constants
-nsite <- 20
-nobs_count <- 100
-# params
-mu <- rlnorm(nsite,meanlog=log(1),sdlog=1)
-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)
-}
-}
-}
-# collect data
-data <- list(
-count = count,
-count.type = count_type
-)
-# initial values
-inits <- list()
-inits[[1]] <- list(
-mu = mu
-)
-names(inits[[1]]) <- c('mu')
-# run model
-fit <- traditionalModel(data=data, q=TRUE,
-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('q[1]') %in% output_params))
-# constants
-nsite <- 20
-nobs_count <- 100
-# params
-mu <- rlnorm(nsite,meanlog=log(1),sdlog=1)
-q <- 2
-beta_gamma <- 1
-alpha_gamma <- mu * beta_gamma
-# traditional type
-count_type <- cbind(matrix(1,nrow=nsite,ncol=nobs_count/2),
-matrix(2,nrow=nsite,ncol=nobs_count/2))
-# collect data
-data <- list(
-count = count,
-count.type = count_type
-)
-# initial values
-inits <- list()
-inits[[1]] <- list(
-alpha_gamma = mu,
-beta_gamma = rep(1,length(mu))
-)
-names(inits[[1]]) <- c('alpha_gamma','beta_gamma')
-# run model
-fit <- traditionalModel(data=data, q=TRUE,family='gamma',
-n.chain=1, multicore=FALSE, seed = 10,
-initial_values=inits)
-# get output params
-output_params <- rownames(as.data.frame(jointSummarize(fit$model)))
-fit$model
-# constants
-nsite <- 20
-nobs_count <- 100
-# params
-mu <- rlnorm(nsite,meanlog=log(1),sdlog=1)
-q <- 2
-beta_gamma <- 1
-alpha_gamma <- mu * beta_gamma
-# traditional type
-count_type <- cbind(matrix(1,nrow=nsite,ncol=nobs_count/2),
-matrix(2,nrow=nsite,ncol=nobs_count/2))
-# collect data
-data <- list(
-count = count,
-count.type = count_type
-)
-# initial values
-inits <- list()
-inits[[1]] <- list(
-alpha = mu,
-beta = rep(1,length(mu)),
-q=q
-)
-names(inits[[1]]) <- c('alpha','beta','q')
-# run model
-fit <- traditionalModel(data=data, q=TRUE,family='gamma',
-n.chain=1, multicore=FALSE, seed = 10,
-initial_values=inits)
-devtools::load_all()
-# get output params
-output_params <- rownames(as.data.frame(jointSummarize(fit$model)))
-fit$model
-modelfit <- fit$model
-all(c('q','alpha','beta') %in% modelfit@model_pars)==TRUE &&
-all(!c('p10','beta','phi') %in% modelfit@model_pars==TRUE) &&
-all(par == 'all')
-modelfit@model_pars
-par <- 'all'
-all(c('q','alpha','beta') %in% modelfit@model_pars)==TRUE &&
-all(!c('p10','beta','phi') %in% modelfit@model_pars==TRUE) &&
-all(par == 'all')
-all(c('q','alpha','beta') %in% modelfit@model_pars)==TRUE
-all(!c('p10','beta','phi') %in% modelfit@model_pars==TRUE)
-all(c('q','alpha','beta') %in% modelfit@model_pars)==TRUE &&
-all(!c('p10','phi') %in% modelfit@model_pars==TRUE) &&
-all(par == 'all')
-devtools::load_all()
-# get output params
-output_params <- rownames(as.data.frame(jointSummarize(fit$model)))
-# test expectation
-expect_true(all(c('q[1]') %in% output_params))
-# constants
-nsite <- 20
-nobs_count <- 100
-# params
-mu <- rlnorm(nsite,meanlog=log(1),sdlog=1)
-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)
-}
-# collect data
-data <- list(
-count = count
-)
-# initial values
-inits <- list()
-inits[[1]] <- list(
-mu = mu,
-phi = phi
-)
-names(inits[[1]]) <- c('mu','phi')
-# run model
-fit <- traditionalModel(data=data,family='negbin',
-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('phi') %in% output_params))
-## 17.
-# model, pois (traditional model)
-# constants
-nsite <- 20
-nobs_count <- 100
-# params
-mu <- rlnorm(nsite,meanlog=log(1),sdlog=1)
-# count
-count <- matrix(NA,nrow=nsite,ncol=nobs_count)
-for(i in 1:nsite){
-count[i,] <- rpois(n=nobs_count,mu[i])
-}
-# collect data
-data <- list(
-count = count
-)
-# initial values
-inits <- list()
-inits[[1]] <- list(
-mu = mu
-)
-names(inits[[1]]) <- c('mu')
-# run model
-fit <- traditionalModel(data=data,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','beta','q','phi') %in% output_params))
-devtools::load_all()
-## 18.
-# model, gamma (traditional model)
-# constants
-nsite <- 20
-nobs_count <- 100
-# params
-mu <- rlnorm(nsite,meanlog=log(1),sdlog=1)
-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)
-}
-# collect data
-data <- list(
-count = count
-)
-# initial values
-inits <- list()
-inits[[1]] <- list(
-alpha = mu,
-beta = rep(1,length(mu))
-)
-names(inits[[1]]) <- c('alpha','beta')
-# run model
-fit <- traditionalModel(data=data,n.chain=1, family='gamma',
-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','q','phi') %in% output_params))
+fit1 <- jointModel(data=data, cov=c('var_a','var_b'),
+multicore=FALSE, seed = 10)
 roxygen2::roxygenise()
-library(covr)
-package_coverage()
-getwd()
-cov <- package_coverage("/home/abby/eDNAjoint")
-posterior <- read.csv('../eDNA_presentation/Joint_posterior.csv')
-library(tidyverse)
-ggplot()+
-geom_density(aes(x=posterior$p10))
-ggplot()+
-geom_density(aes(x=exp(posterior$p10)))
-rounded <- round(posterior$p10,7)
-head(roudned)
-head(rounded)
-table(rounded)
-rounded <- round(posterior$p10,3)
-table(rounded)
-names(table(rounded))
-names <- as.numeric(names(table(rounded)))
-names
-values <- as.vector(tables(rounded))
-values <- as.vector(table(rounded))
-ggplot()+
-geom_line(aes(x=names,y=values))
-ggplot()+
-geom_line(aes(x=names,y=exp(values)))
-rounded <- round(posterior$p10,7)
-names <- as.numeric(names(table(rounded)))
-values <- as.vector(table(rounded))
-ggplot()+
-geom_line(aes(x=names,y=exp(values)))
-ggplot()+
-geom_line(aes(x=exp(names),y=values))
-rounded <- round(posterior$p10,3)
-names <- as.numeric(names(table(rounded)))
-values <- as.vector(table(rounded))
-ggplot()+
-geom_line(aes(x=exp(names),y=values))
-ggplot()+
-geom_smooth(aes(x=exp(names),y=values))
-ggplot()+
-geom_density(aes(x=exp(posterior$p10)))
-ggplot()+
-geom_density(aes(x=exp(posterior$p10)))+
-labs(x='probability of false positive eDNA detection',
-y='density')+
-theme_minimal()
-ggplot()+
-geom_density(aes(x=exp(posterior$p10)),linewidth=1)+
-labs(x='probability of false positive eDNA detection',
-y='density')+
-theme_minimal()
-ggplot()+
-geom_density(aes(x=exp(posterior$p10)),linewidth=1)+
-labs(x='probability of false positive eDNA detection',
-y='density')+
-theme_minimal()+
-theme(text=element_text(size=12))
-p10_plot <- ggplot()+
-geom_density(aes(x=exp(posterior$p10)),linewidth=1)+
-labs(x='probability of false positive eDNA detection',
-y='density')+
-theme_minimal()+
-theme(text=element_text(size=12))
-ggsave('../eDNA_presentation/p10_plot.png',p10_plot,dpi=400)
-raa_plot <- ggplot()+
-geom_density(aes(x=posterior$mu_crab.15),linewidth=1)+
-labs(x='EGC density (crabs/trap)',
-y='density')+
-theme_minimal()+
-theme(text=element_text(size=12))
-raa_plot
-raa_plot <- ggplot()+
-geom_density(aes(x=posterior$mu.15),linewidth=1)+
-labs(x='EGC density (crabs/trap)',
-y='density')+
-theme_minimal()+
-theme(text=element_text(size=12))
-raa_plot
-raa_plot <- ggplot()+
-geom_density(aes(x=posterior$mu.15),linewidth=1)+
-scale_x_continuous(limits=c(0,1.5))+
-labs(x='EGC density (crabs/trap)',
-y='density')+
-theme_minimal()+
-theme(text=element_text(size=12))
-raa_plot
-ggsave('../eDNA_presentation/raa_plot.png',p10_plot,dpi=400)
-ggsave('../eDNA_presentation/raa_plot.png',raa_plot,dpi=400)
-cov <- package_coverage("/home/abby/eDNAjoint")
-cov$`detectionCalculate.R:429:7:429:70:7:70:431:431`
-dim(cov)
-str(cov)
-cov$`detectionCalculate.R:429:7:429:70:7:70:431:431`
-length(cov)
-test <- as.data.frame(cov)
-View(test)
-table(test$value)
-zero_coverage(cov)
-getwd()
-saveRDS(test,'../eDNAjoint_coverage.rds')
+options(repos = c(
+ropenscireviewtools = "https://ropensci-review-tools.r-universe.dev",
+CRAN = "https://cloud.r-project.org"))
+install.packages("srr")
 roxygen2::roxygenise()
-devtools::load_all()
-library(eDNAjoint)
-#run the joint model with two covariates
-#' @srrstats {G1.5} Example code that reproduces results in the publication
-#' (Keller et al., 2022) where the model/algorithm was first developed
-#' @srrstats {BS1.2c} documentation of jointModel() function
-goby.fit1 <- jointModel(data = gobyData, cov=c('Filter_time','Salinity'),
-family = 'poisson', p10priors = c(1,20), q=FALSE,
-multicore=FALSE)
-# fit a new model with one site-level covariate
-goby.fit2 <- jointModel(data = gobyData, cov='Other_fishes',
-family = 'poisson', p10priors = c(1,20), q=FALSE,
-multicore=FALSE)
-#' @srrstats {BS2.7,BS2.11} Example of providing initial values
-# set number of chains
-n.chain <- 4
-# initial values should be a list of named lists
-inits <- list()
-for(i in 1:n.chain){
-inits[[i]] <- list(
-# length should equal the number of sites (dim(gobyData$count)[1])
-# for each chain
-mu = stats::runif(dim(gobyData$count)[1], 0.01, 5),
-# length should equal 1 for each chain
-p10 = stats::runif(1,log(0.0001),log(0.08)),
-# length should equal the number of covariates plus 1
-# (to account for intercept in regression)
-alpha = rep(0.1,length(c('Filter_time','Salinity'))+1)
-)
-names(inits[[i]]) <- c('mu','p10','alpha')
-}
-# now fit the model
-fit.w.inits <- jointModel(data = gobyData, cov=c('Filter_time','Salinity'),
-initial_values = inits, multicore=FALSE)
-#run the joint model with catchability coefficients
-greencrab.fit.q.negbin <- jointModel(data = greencrabData, cov=NULL,
-family = 'negbin',
-p10priors = c(1,20), q=TRUE,
-multicore=FALSE)
-help(jointModel)
-library(eDNAjoint)
-devtools::install()
-library(eDNAjoint)
-devtools::install()
-devtools::install()
-remove.packages('eDNAjoint')
-.libPaths()
-devtools::install()
 devtools::install()
-library(eDNAjoint)
-help(jointModel)
 library(autotest)
-x <- autotest_package(package='eDNAjoint',functions='jointModel',test=FALSE)
-View(x)
-roxygen2::roxygenise()
-devtools::install()
-library(eDNAjoint)
-help(jointModel)
-help(traditionalModel)
-help(jointModel)
 library(eDNAjoint)
-help(jointModel)
-library(autotest)
 x <- autotest_package(package='eDNAjoint',functions='jointModel',test=FALSE)
-View(x)
-getwd()
-saveRDS(x,'../autotest_eDNAjoint/tests_jointModel.rds')
-unique(x$test_name)
-subset <- which(x$test_name=='par_is_documented')
-subset
-help(autotest_package)
-install.package('MASS')
-install.packages('MASS')
-test <- x[x$test_name=='par_is_documented',]
-View(test)
-help(jointModel)
-hist(rbeta(1000,1.1,20.1))
-subset
-x_mod <- x
-x_mod$test <- FALSE
-x_mod[subset,'test'] <- TRUE
-View(x_mod)
-y <- autotest_package(package='eDNAjoint',functions='jointModel',test=TRUE,test_data=x_mod)
-y
-summary(y)
-View(y)
-family <- 'Negbin'
-lower(family)
-roxygen2::roxygenise()
-devtools::install()
-library(eDNAjoint)
-help(jointModel)
+warnings()
diff --git a/tests/testthat/test-jointSummarize.R b/tests/testthat/test-jointSummarize.R
index 7d56b1f..a207694 100644
--- a/tests/testthat/test-jointSummarize.R
+++ b/tests/testthat/test-jointSummarize.R
@@ -369,7 +369,8 @@ test_that("jointSummarize outputs work", {
   )
   names(inits[[1]]) <- c('mu','p10','beta','phi')
   # run model
-  fit <- suppressWarnings({jointModel(data=data, family = 'negbin', initial_values=inits,
+  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)})
@@ -1142,7 +1143,7 @@ test_that("jointSummarize outputs work", {
   fit <- suppressWarnings({jointModel(data=data,family='gamma',
                     cov=c('var_a','var_b'),
                     n.chain=1, multicore=FALSE, seed = 10,
-                    initial_values=initsn.iter.burn = 25,
+                    initial_values=inits, n.iter.burn = 25,
                     n.iter.sample = 75)})
 
   # get output params