diff --git a/.Rhistory b/.Rhistory
index 3f074b6..0841f6e 100644
--- a/.Rhistory
+++ b/.Rhistory
@@ -1,180 +1,3 @@
-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 <- 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 <- 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: '),
@@ -510,3 +333,180 @@ ggplot2::geom_point(ggplot2::aes(x = c(1,2), y = c(1,2)))+
 ggplot2::ggtitle(paste0('# survey units necessary to detect',
 ' species presence, given mu'))+
 ggplot2::theme_minimal()
+library(eDNAjoint)
+library(eDNAjoint)
+data(gobyData)
+goby.fit <- jointModel(data = gobyData)
+library(eDNAjoint)
+data(gobyData)
+# fit model
+goby.fit <- jointModel(data = gobyData)
+# summarize false positive probability (p10) probability
+jointSummarize(goby.fit$model, par = 'p10')
+library(eDNAjoint)
+data(gobyData)
+# run the joint model with two covariates
+goby.fit <- jointModel(data = gobyData,
+cov = c('Filter_time','Salinity'))
+# summarize false positive probability (p10) probability
+jointSummarize(goby.fit$model, par = 'p10')
+# summarize false positive probability (p10) probability
+jointSummarize(goby.fit$model, par = 'p10')
+library(eDNAjoint)
+data(gobyData)
+# run the joint model with two covariates
+goby.fit <- jointModel(data = gobyData,
+cov = c('Filter_time','Salinity'))
+# summarize false positive probability probability
+jointSummarize(goby.fit$model, par = 'p10')
+detectionPlot(goby.fit$model,
+mu.min = 0.1, mu.max = 1,
+cov.val = c(0,0))
+# plot effort necessary to detect presence
+detectionPlot(goby.fit$model,
+mu.min = 0.1, mu.max = 1,
+cov.val = c(0,0))
+# plot effort necessary to detect presence
+# at turbid sites
+detectionPlot(goby.fit$model,
+mu.min = 0.1, mu.max = 1,
+cov.val = c(1,0))
+A
+modelfit <- goby.fit$model
+mu.min = 0.01
+mu.max <- 1
+# create mu vector
+mu <- seq(from = mu.min, to = mu.max, by = (mu.max-mu.min)/1000)
+cov.val <- c(0,0)
+probability <- 0.9
+qPCR.N <- 3
+# get n samples df
+out <- detectionCalculate(modelfit, mu, cov.val, probability, qPCR.N)
+'%>%' <- magrittr::`%>%`
+out_long <- as.data.frame(out) %>%
+tidyr::pivot_longer(cols =! mu, names_to = 'survey_type')
+length(unique(out_long$survey_type))
+unique(out_long$survey_type)
+sort(unique(out_long$survey_type))
+# get full names
+if(length(unique(out_long$survey_type))==1){
+names <- 'traditional'
+} else if(length(unique(out_long$survey_type))==2){
+names <- c('eDNA','traditional')
+} else{
+names <- 'eDNA'
+for(i in 1:(length(unique(out_long$survey_type))-1)){
+names <- c(names, paste0('traditional (gear type ',i,')'))
+}
+}
+ggplot2::ggplot(data = out_long)+
+ggplot2::geom_line(ggplot2::aes(x = mu, y = value, color = survey_type),
+linewidth = 1)+
+ggplot2::labs(x = 'mu (expected catch rate)',y = '# survey units',
+color = 'survey type')+
+ggplot2::scale_color_manual(values = sort(unique(out_long$survey_type)),
+names = names)+
+ggplot2::theme_minimal()
+ggplot2::ggplot(data = out_long)+
+ggplot2::geom_line(ggplot2::aes(x = mu, y = value, color = survey_type),
+linewidth = 1)+
+ggplot2::labs(x = 'mu (expected catch rate)',y = '# survey units',
+color = 'survey type')+
+ggplot2::scale_color_manual(names = sort(unique(out_long$survey_type)),
+values = names)+
+ggplot2::theme_minimal()
+ggplot2::ggplot(data = out_long)+
+ggplot2::geom_line(ggplot2::aes(x = mu, y = value, color = survey_type),
+linewidth = 1)+
+ggplot2::labs(x = 'mu (expected catch rate)',y = '# survey units',
+color = 'survey type')+
+ggplot2::scale_color_manual(values = sort(unique(out_long$survey_type)),
+labels = names)+
+ggplot2::theme_minimal()
+ggplot2::ggplot(data = out_long)+
+ggplot2::geom_line(ggplot2::aes(x = mu, y = value, color = survey_type),
+linewidth = 1)+
+ggplot2::labs(x = 'mu (expected catch rate)',y = '# survey units',
+color = 'survey type')+
+ggplot2::scale_color_manual(labels = names)+
+ggplot2::theme_minimal()
+names
+scales::hue_pal()
+scales::hue_pal()[1:2]
+scales::hue_pal(2)
+scales::hue_pal()(2)
+ggplot2::ggplot(data = out_long)+
+ggplot2::geom_line(ggplot2::aes(x = mu, y = value, color = survey_type),
+linewidth = 1)+
+ggplot2::labs(x = 'mu (expected catch rate)',y = '# survey units',
+color = 'survey type')+
+ggplot2::scale_color_manual(values = scales::hue_pal()(2),
+labels = names)+
+ggplot2::theme_minimal()
+out
+mu.min
+mu.min <- 0.1
+# create mu vector
+mu <- seq(from = mu.min, to = mu.max, by = (mu.max-mu.min)/1000)
+# get n samples df
+out <- detectionCalculate(modelfit, mu, cov.val, probability, qPCR.N)
+# convert to long df
+'%>%' <- magrittr::`%>%`
+out_long <- as.data.frame(out) %>%
+tidyr::pivot_longer(cols =! mu, names_to = 'survey_type')
+# get full names
+if(length(unique(out_long$survey_type))==1){
+names <- 'traditional'
+} else if(length(unique(out_long$survey_type))==2){
+names <- c('eDNA','traditional')
+} else{
+names <- 'eDNA'
+for(i in 1:(length(unique(out_long$survey_type))-1)){
+names <- c(names, paste0('traditional (gear type ',i,')'))
+}
+}
+ggplot2::ggplot(data = out_long)+
+ggplot2::geom_line(ggplot2::aes(x = mu, y = value, color = survey_type),
+linewidth = 1)+
+ggplot2::labs(x = 'mu (expected catch rate)',y = '# survey units',
+color = 'survey type')+
+ggplot2::scale_color_manual(values = scales::hue_pal()(length(names)),
+labels = names)+
+ggplot2::theme_minimal()
+fit.q <- jointModel(data = greencrabData,family = "negbin", q = TRUE,
+n.chain = 1, n.iter.burn = 500, n.iter.sample = 1000)
+cov.val = NULL
+modelfit <- fit.q$model
+# create mu vector
+mu <- seq(from = mu.min, to = mu.max, by = (mu.max-mu.min)/1000)
+# get n samples df
+out <- detectionCalculate(modelfit, mu, cov.val, probability, qPCR.N)
+# convert to long df
+'%>%' <- magrittr::`%>%`
+out_long <- as.data.frame(out) %>%
+tidyr::pivot_longer(cols =! mu, names_to = 'survey_type')
+# get full names
+if(length(unique(out_long$survey_type))==1){
+names <- 'traditional'
+} else if(length(unique(out_long$survey_type))==2){
+names <- c('eDNA','traditional')
+} else{
+names <- 'eDNA'
+for(i in 1:(length(unique(out_long$survey_type))-1)){
+names <- c(names, paste0('traditional (gear type ',i,')'))
+}
+}
+names
+ggplot2::ggplot(data = out_long)+
+ggplot2::geom_line(ggplot2::aes(x = mu, y = value, color = survey_type),
+linewidth = 1)+
+ggplot2::labs(x = 'mu (expected catch rate)',y = '# survey units',
+color = 'survey type')
+ggplot2::ggplot(data = out_long)+
+ggplot2::geom_line(ggplot2::aes(x = mu, y = value, color = survey_type),
+linewidth = 1)+
+ggplot2::labs(x = 'mu (expected catch rate)',y = '# survey units',
+color = 'survey type')+
+ggplot2::scale_color_manual(values = scales::hue_pal()(length(names)),
+labels = names)+
+ggplot2::theme_minimal()
diff --git a/DESCRIPTION b/DESCRIPTION
index 4978c8c..03ac491 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -29,7 +29,8 @@ Imports:
     rlist,
     rstan (>= 2.26.23),
     rstantools (>= 2.3.1.1),
-    tidyr
+    tidyr,
+    scales
 LinkingTo: 
     BH (>= 1.66.0),
     Rcpp (>= 0.12.0),