diff --git a/.Rhistory b/.Rhistory
index 1655775..623bcd9 100644
--- a/.Rhistory
+++ b/.Rhistory
@@ -1,12 +1,512 @@
+# 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 <- mu / (mu + exp(beta))
+p <- pmin(p11 + exp(log_p10), 1)
+# pcr_n (# qPCR observations)
+pcr_n <- matrix(3, nrow = nsite, ncol = nobs_pcr)
+# pcr_k (# positive qPCR detections)
+pcr_k <- matrix(NA, nrow = nsite, ncol = nobs_pcr)
+for (i in 1:nsite) {
+pcr_k[i, ] <- rbinom(nobs_pcr, pcr_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(
+pcr_n = pcr_n,
+pcr_k = pcr_k,
+count = count,
+count_type = count_type
+)
+# initial values
+inits <- list()
+inits[[1]] <- list(
+mu = mu,
+p10 = exp(log_p10),
+alpha = beta,
+q = 2
+)
+names(inits[[1]]) <- c("mu", "p10", "alpha", "q")
+# run model
+fit <- suppressWarnings({
+joint_model(data = data, q = TRUE, n_chain = 1, multicore = FALSE,
+seed = 10, initial_values = inits, n_warmup = 25, n_iter = 100)
+})
+devtools::load_all()
+## 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 <- mu / (mu + exp(beta))
+p <- pmin(p11 + exp(log_p10), 1)
+# pcr_n (# qPCR observations)
+pcr_n <- matrix(3, nrow = nsite, ncol = nobs_pcr)
+# pcr_k (# positive qPCR detections)
+pcr_k <- matrix(NA, nrow = nsite, ncol = nobs_pcr)
+for (i in 1:nsite) {
+pcr_k[i, ] <- rbinom(nobs_pcr, pcr_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(
+pcr_n = pcr_n,
+pcr_k = pcr_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),
+alpha = beta,
+q = q
+)
+names(inits[[1]]) <- c("alpha_gamma", "beta_gamma", "p10", "alpha", "q")
+# run model
+fit <- suppressWarnings({
+joint_model(data = data, q = TRUE, family = "gamma",
+n_chain = 1, multicore = FALSE, seed = 10,
+initial_values = inits, n_warmup = 25, n_iter = 100)
+})
+devtools::load_all()
+## 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 <- mu / (mu + exp(beta))
+p <- pmin(p11 + exp(log_p10), 1)
+# pcr_n (# qPCR observations)
+pcr_n <- matrix(3, nrow = nsite, ncol = nobs_pcr)
+# pcr_k (# positive qPCR detections)
+pcr_k <- matrix(NA, nrow = nsite, ncol = nobs_pcr)
+for (i in 1:nsite) {
+pcr_k[i, ] <- rbinom(nobs_pcr, pcr_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(
+pcr_n = pcr_n,
+pcr_k = pcr_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),
+alpha = beta,
+q = q
+)
+names(inits[[1]]) <- c("alpha_gamma", "beta_gamma", "p10", "alpha", "q")
+# run model
+fit <- suppressWarnings({
+joint_model(data = data, q = TRUE, family = "gamma",
+n_chain = 1, multicore = FALSE, seed = 10,
+initial_values = inits, n_warmup = 25, n_iter = 100)
+})
+devtools::load_all()
+qPCR.N <- matrix(3, nrow=12, ncol=12)
+colnames(qPCR.N) <- c("1","2","3","4","5","6","7","8","9","10","11","12")
+# this shows that for every one of 12 sites, there were 12 samples, and 3 replicates were collected for each sample
+RHCAqPCR.K <- matrix(nrow=12, ncol=12)
+colnames(RHCAqPCR.K) <- c("1","2","3","4","5","6","7","8","9","10","11","12")
+RHKEqPCR.K <- matrix(nrow=12, ncol=12)
+colnames(RHKEqPCR.K) <- c("1","2","3","4","5","6","7","8","9","10","11","12")
+# made empty 12 x 12 matrices, now let's fill them in with the data
+RHCAqPCR.K[1, ] <- c(0,0,0,0,0,0,0,0,0,0,0,0) #badger had zero positive replicates
+RHCAqPCR.K[2, ] <- c(3,0,0,2,0,0,0,0,0,0,0,0) #bear
+RHCAqPCR.K[3, ] <- c(0,0,1,1,2,0,1,1,1,1,0,0) #bolin
+RHCAqPCR.K[4, ] <- c(3,3,2,1,2,3,3,2,3,1,2,3) #bridal
+RHCAqPCR.K[5, ] <- c(2,1,2,1,3,1,3,0,1,0,1,1) #flat
+RHCAqPCR.K[6, ] <- c(3,3,3,3,3,3,3,3,3,3,3,3) #larch had all replicates positive
+RHCAqPCR.K[7, ] <- c(0,0,0,0,0,0,0,0,0,0,0,0) #salmon had zero positive replicates
+RHCAqPCR.K[8, ] <- c(1,0,3,1,1,0,1,1,1,0,2,1) #eagle
+RHCAqPCR.K[9, ] <- c(2,0,0,0,1,0,0,0,0,0,0,1) #santiam
+RHCAqPCR.K[10, ] <- c(0,1,0,1,3,1,1,0,0,0,2,1) #strawberry
+RHCAqPCR.K[11, ] <- c(3,3,0,0,0,0,0,3,3,1,3,3) #whiskey
+RHCAqPCR.K[12, ] <- c(0,0,3,2,1,2,0,0,0,0,0,1) #wyeth
+RHKEqPCR.K[1, ] <- c(1,3,3,3,3,3,3,3,2,3,3,3) #101
+RHKEqPCR.K[2, ] <- c(2,3,1,1,1,1,2,2,3,3,2,3) #clarence
+RHKEqPCR.K[3, ] <- c(2,2,1,0,1,1,1,0,1,3,1,1) #coffee
+RHKEqPCR.K[4, ] <- c(3,3,3,3,3,3,3,3,3,3,3,3) #crazy had all replicates positive
+RHKEqPCR.K[5, ] <- c(1,0,0,3,0,1,1,1,0,2,1,2) #edwards
+RHKEqPCR.K[6, ] <- c(3,3,3,3,3,3,3,3,3,3,3,3) #elochoman had all replicates positive
+RHKEqPCR.K[7, ] <- c(3,1,3,3,3,3,3,2,0,2,3,2) #falk
+RHKEqPCR.K[8, ] <- c(2,1,3,3,3,2,0,1,3,3,2,3) #fern
+RHKEqPCR.K[9, ] <- c(3,3,3,3,3,3,3,3,3,3,3,3) #holcombe had all replicates positive
+RHKEqPCR.K[10, ] <- c(3,3,3,3,3,3,3,3,3,3,3,3) #nehalem had all replicates positive
+RHKEqPCR.K[11, ] <- c(3,0,1,3,1,3,2,3,3,3,2,2) #rock
+RHKEqPCR.K[12, ] <- c(1,0,0,0,0,0,3,2,3,3,2,0) #wlest
+RHCAcount <- matrix(nrow=12, ncol=1)
+colnames(RHCAcount) <- c("1")
+RHKEcount <- matrix(nrow=12, ncol=1)
+colnames(RHKEcount) <- c("1")
+# 12 rows, one for each site, only 1 column because we only did physical samples once per site
+RHCAcount[1] <- c(1) #badger
+RHCAcount[2] <- c(2) #bear
+RHCAcount[3] <- c(8) #bolin
+RHCAcount[4] <- c(2) #bridal
+RHCAcount[5] <- c(8) #flat
+RHCAcount[6] <- c(15) #larch
+RHCAcount[7] <- c(13) #salmon
+RHCAcount[8] <- c(7) #eagle
+RHCAcount[9] <- c(11) #santiam
+RHCAcount[10] <- c(6) #strawberry
+RHCAcount[11] <- c(19) #whiskey
+RHCAcount[12] <- c(3) #wyeth
+RHKEcount[1] <- c(24) #101
+RHKEcount[2] <- c(4) #clarence
+RHKEcount[3] <- c(6) #coffee
+RHKEcount[4] <- c(17) #crazy
+RHKEcount[5] <- c(9) #edwards
+RHKEcount[6] <- c(16) #elochoman
+RHKEcount[7] <- c(13) #falk
+RHKEcount[8] <- c(29) #fern
+RHKEcount[9] <- c(53) #holcombe
+RHKEcount[10] <- c(11) #nehalem
+RHKEcount[11] <- c(11) #rock
+RHKEcount[12] <- c(26) #wlest
+RHCAcount2 <- matrix(nrow=12, ncol=2)
+colnames(RHCAcount2) <- c("1","2")
+RHKEcount2 <- matrix(nrow=12, ncol=2)
+colnames(RHKEcount2) <- c("1","2")
+# 12 rows, one for each site, only 1 column because we only did physical samples once per site
+RHCAcount2[1,] <- c(1,1) #badger
+RHCAcount2[2,] <- c(2,2) #bear
+RHCAcount2[3,] <- c(8,8) #bolin
+RHCAcount2[4,] <- c(2,2) #bridal
+RHCAcount2[5,] <- c(8,8) #flat
+RHCAcount2[6,] <- c(15,15) #larch
+RHCAcount2[7,] <- c(13,13) #salmon
+RHCAcount2[8,] <- c(7,7) #eagle
+RHCAcount2[9,] <- c(11,11) #santiam
+RHCAcount2[10,] <- c(6,6) #strawberry
+RHCAcount2[11,] <- c(19,19) #whiskey
+RHCAcount2[12,] <- c(3,3) #wyeth
+RHKEcount2[1,] <- c(24,24) #101
+RHKEcount2[2,] <- c(4,4) #clarence
+RHKEcount2[3,] <- c(6,6) #coffee
+RHKEcount2[4,] <- c(17,17) #crazy
+RHKEcount2[5,] <- c(9,9) #edwards
+RHKEcount2[6,] <- c(16,16) #elochoman
+RHKEcount2[7,] <- c(13,13) #falk
+RHKEcount2[8,] <- c(29,29) #fern
+RHKEcount2[9,] <- c(53,53) #holcombe
+RHKEcount2[10,] <- c(11,11) #nehalem
+RHKEcount2[11,] <- c(11,11) #rock
+RHKEcount2[12,] <- c(26,26) #wlest
+# AK changed
+RHCA_list <- list(
+pcr_n = qPCR.N,  # Number of eDNA qPCR replicates per site
+pcr_k = RHCAqPCR.K, # Number of positive eDNA qPCR detections per site
+count = RHCAcount  # Traditional survey count data
+)
+# AK changed
+RHKE_list <- list(
+pcr_n = qPCR.N,  # Number of eDNA qPCR replicates per site
+pcr_k = RHKEqPCR.K, # Number of positive eDNA qPCR detections per site
+count = RHKEcount  # Traditional survey count data
+)
+# AK changed
+RHCA_list2 <- list(
+pcr_n = qPCR.N,  # Number of eDNA qPCR replicates per site
+pcr_k = RHCAqPCR.K, # Number of positive eDNA qPCR detections per site
+count = RHCAcount2  # Traditional survey count data
+)
+RHKE_list2 <- list(
+pcr_n = qPCR.N,  # Number of eDNA qPCR replicates per site
+pcr_k = RHKEqPCR.K, # Number of positive eDNA qPCR detections per site
+count = RHKEcount2  # Traditional survey count data
+)
+RHCA.fit1 <- joint_model(data = RHCA_list, family = 'poisson',
+p10_priors = c(0.1,99.99), q=FALSE,
+n_warmup = 1000, n_iter = 8000,
+adapt_delta = 0.9999, multicore = TRUE)
+RHKE.fit1 <- joint_model(data = RHKE_list, family = 'poisson',
+p10_priors = c(0.1,99.9), q=FALSE,
+n_warmup = 1000, n_iter = 8000,
+adapt_delta = 0.9999, multicore = TRUE)
+joint_summarize(RHKE.fit1)
+joint_summarize(RHKE.fit1$model)
+RHKE_list
+ggplot()+
+geom_point(aes(x=rowMeans(RHCA_list$count, na.rm = TRUE),
+y=rowMeans(RHCA_list$pcr_k, na.rm = TRUE)),
+size = 3)+
+labs(x='mean count',y='mean positive eDNA detections')+
+theme_minimal()
+library(tidyverse)
+ggplot()+
+geom_point(aes(x=rowMeans(RHCA_list$count, na.rm = TRUE),
+y=rowMeans(RHCA_list$pcr_k, na.rm = TRUE)),
+size = 3)+
+labs(x='mean count',y='mean positive eDNA detections')+
+theme_minimal()
+ggplot()+
+geom_point(aes(x=rowMeans(RHKE_list$count, na.rm = TRUE),
+y=rowMeans(RHKE_list$pcr_k, na.rm = TRUE)),
+size = 3)+
+labs(x='mean count',y='mean positive eDNA detections')+
+theme_minimal()
+jointSummarize(RHKE.fit1$model)
+joint_summarize(RHKE.fit1$model)
+hist(rbeta(10000, 0.1, 99.9))
+mean(rbeta(10000, 0.1, 99.9))
+var(rbeta(10000, 0.1, 99.9))
+mean(rbeta(10000, 1, 100))
+mean(rbeta(10000, 1, 1000))
+var(rbeta(10000, 1, 1000))
+RHKE.fit1 <- joint_model(data = RHKE_list, family = 'poisson',
+p10_priors = c(1, 1000), q=FALSE,
+n_warmup = 1000, n_iter = 8000,
+adapt_delta = 0.9999, multicore = TRUE)
+mean1 <- mean(rbeta(10000, 0.1, 99.9))
+var1 <- var(rbeta(10000, 0.1, 99.9))
+print('mean: ', mean1)
+mean1 <- mean(rbeta(10000, 0.1, 99.9))
+var1 <- var(rbeta(10000, 0.1, 99.9))
+print('mean: ', mean1)
+mean1 <- mean(rbeta(10000, 0.1, 99.9))
+var1 <- var(rbeta(10000, 0.1, 99.9))
+print(paste0('mean: ', mean1))
+print(paste0('var: ', var1))
+mean2 <- mean(rbeta(10000, 1, 1000))
+var2 <- var(rbeta(10000, 0.1, 1000))
+print(paste0('mean: ', mean2))
+print(paste0('var: ', var2))
+library(tidyverse)
+ggplot() +
+geom_histogram(x = rbeta(10000, 0.1, 99.9), color = "blue") +
+geom_histogram(x = rbeta(10000, 1, 1000), color = "purple") +
+theme_minimal()
+ggplot() +
+geom_histogram(aes(x = rbeta(10000, 0.1, 99.9), color = "blue")) +
+geom_histogram(aes(x = rbeta(10000, 1, 1000), color = "purple")) +
+theme_minimal()
+library(tidyverse)
+ggplot() +
+geom_histogram(aes(x = rbeta(10000, 0.1, 99.9),
+fill = "blue", alpha = 0.4)) +
+geom_histogram(aes(x = rbeta(10000, 1, 1000),
+fill = "purple", alpha = 0.4)) +
+theme_minimal()
+library(tidyverse)
+ggplot() +
+geom_histogram(aes(x = rbeta(10000, 0.1, 99.9),
+fill = "blue", alpha = 0.4)) +
+geom_histogram(aes(x = rbeta(10000, 1, 1000),
+fill = "purple", alpha = 0.4)) +
+theme_minimal() +
+theme(legend.position = "None")
+ggplot() +
+geom_histogram(aes(x = rbeta(10000, 0.1, 99.9),
+fill = "blue", alpha = 0.4)) +
+geom_histogram(aes(x = rbeta(10000, 1, 1000),
+fill = "purple", alpha = 0.4)) +
+labs(x = "p10 value", y = "frequency") +
+theme_minimal() +
+theme(legend.position = "None")
+remove.packages('eDNAjoint')
+install.packages("eDNAjoint", repos = "https://ropensci.r-universe.dev")
+install.packages("eDNAjoint", repos = "https://ropensci.r-universe.dev")
+install.packages("eDNAjoint", repos = "https://ropensci.r-universe.dev")
+library(eDNAjoint)
+#install.packages("eDNAjoint")
+library(eDNAjoint)
+library(bayesplot)
+qPCR.N <- matrix(3, nrow=12, ncol=12)
+colnames(qPCR.N) <- c("1","2","3","4","5","6","7","8","9","10","11","12")
+# this shows that for every one of 12 sites, there were 12 samples, and 3 replicates were collected for each sample
+RHCAqPCR.K <- matrix(nrow=12, ncol=12)
+colnames(RHCAqPCR.K) <- c("1","2","3","4","5","6","7","8","9","10","11","12")
+RHKEqPCR.K <- matrix(nrow=12, ncol=12)
+colnames(RHKEqPCR.K) <- c("1","2","3","4","5","6","7","8","9","10","11","12")
+# made empty 12 x 12 matrices, now let's fill them in with the data
+RHCAqPCR.K[1, ] <- c(0,0,0,0,0,0,0,0,0,0,0,0) #badger had zero positive replicates
+RHCAqPCR.K[2, ] <- c(3,0,0,2,0,0,0,0,0,0,0,0) #bear
+RHCAqPCR.K[3, ] <- c(0,0,1,1,2,0,1,1,1,1,0,0) #bolin
+RHCAqPCR.K[4, ] <- c(3,3,2,1,2,3,3,2,3,1,2,3) #bridal
+RHCAqPCR.K[5, ] <- c(2,1,2,1,3,1,3,0,1,0,1,1) #flat
+RHCAqPCR.K[6, ] <- c(3,3,3,3,3,3,3,3,3,3,3,3) #larch had all replicates positive
+RHCAqPCR.K[7, ] <- c(0,0,0,0,0,0,0,0,0,0,0,0) #salmon had zero positive replicates
+RHCAqPCR.K[8, ] <- c(1,0,3,1,1,0,1,1,1,0,2,1) #eagle
+RHCAqPCR.K[9, ] <- c(2,0,0,0,1,0,0,0,0,0,0,1) #santiam
+RHCAqPCR.K[10, ] <- c(0,1,0,1,3,1,1,0,0,0,2,1) #strawberry
+RHCAqPCR.K[11, ] <- c(3,3,0,0,0,0,0,3,3,1,3,3) #whiskey
+RHCAqPCR.K[12, ] <- c(0,0,3,2,1,2,0,0,0,0,0,1) #wyeth
+RHKEqPCR.K[1, ] <- c(1,3,3,3,3,3,3,3,2,3,3,3) #101
+RHKEqPCR.K[2, ] <- c(2,3,1,1,1,1,2,2,3,3,2,3) #clarence
+RHKEqPCR.K[3, ] <- c(2,2,1,0,1,1,1,0,1,3,1,1) #coffee
+RHKEqPCR.K[4, ] <- c(3,3,3,3,3,3,3,3,3,3,3,3) #crazy had all replicates positive
+RHKEqPCR.K[5, ] <- c(1,0,0,3,0,1,1,1,0,2,1,2) #edwards
+RHKEqPCR.K[6, ] <- c(3,3,3,3,3,3,3,3,3,3,3,3) #elochoman had all replicates positive
+RHKEqPCR.K[7, ] <- c(3,1,3,3,3,3,3,2,0,2,3,2) #falk
+RHKEqPCR.K[8, ] <- c(2,1,3,3,3,2,0,1,3,3,2,3) #fern
+RHKEqPCR.K[9, ] <- c(3,3,3,3,3,3,3,3,3,3,3,3) #holcombe had all replicates positive
+RHKEqPCR.K[10, ] <- c(3,3,3,3,3,3,3,3,3,3,3,3) #nehalem had all replicates positive
+RHKEqPCR.K[11, ] <- c(3,0,1,3,1,3,2,3,3,3,2,2) #rock
+RHKEqPCR.K[12, ] <- c(1,0,0,0,0,0,3,2,3,3,2,0) #wlest
+RHCAcount <- matrix(nrow=12, ncol=1)
+colnames(RHCAcount) <- c("1")
+RHKEcount <- matrix(nrow=12, ncol=1)
+colnames(RHKEcount) <- c("1")
+# 12 rows, one for each site, only 1 column because we only did physical samples once per site
+RHCAcount[1] <- c(1) #badger
+RHCAcount[2] <- c(2) #bear
+RHCAcount[3] <- c(8) #bolin
+RHCAcount[4] <- c(2) #bridal
+RHCAcount[5] <- c(8) #flat
+RHCAcount[6] <- c(15) #larch
+RHCAcount[7] <- c(13) #salmon
+RHCAcount[8] <- c(7) #eagle
+RHCAcount[9] <- c(11) #santiam
+RHCAcount[10] <- c(6) #strawberry
+RHCAcount[11] <- c(19) #whiskey
+RHCAcount[12] <- c(3) #wyeth
+RHKEcount[1] <- c(24) #101
+RHKEcount[2] <- c(4) #clarence
+RHKEcount[3] <- c(6) #coffee
+RHKEcount[4] <- c(17) #crazy
+RHKEcount[5] <- c(9) #edwards
+RHKEcount[6] <- c(16) #elochoman
+RHKEcount[7] <- c(13) #falk
+RHKEcount[8] <- c(29) #fern
+RHKEcount[9] <- c(53) #holcombe
+RHKEcount[10] <- c(11) #nehalem
+RHKEcount[11] <- c(11) #rock
+RHKEcount[12] <- c(26) #wlest
+RHCAcount2 <- matrix(nrow=12, ncol=2)
+colnames(RHCAcount2) <- c("1","2")
+RHKEcount2 <- matrix(nrow=12, ncol=2)
+colnames(RHKEcount2) <- c("1","2")
+# 12 rows, one for each site, only 1 column because we only did physical samples once per site
+RHCAcount2[1,] <- c(1,1) #badger
+RHCAcount2[2,] <- c(2,2) #bear
+RHCAcount2[3,] <- c(8,8) #bolin
+RHCAcount2[4,] <- c(2,2) #bridal
+RHCAcount2[5,] <- c(8,8) #flat
+RHCAcount2[6,] <- c(15,15) #larch
+RHCAcount2[7,] <- c(13,13) #salmon
+RHCAcount2[8,] <- c(7,7) #eagle
+RHCAcount2[9,] <- c(11,11) #santiam
+RHCAcount2[10,] <- c(6,6) #strawberry
+RHCAcount2[11,] <- c(19,19) #whiskey
+RHCAcount2[12,] <- c(3,3) #wyeth
+RHKEcount2[1,] <- c(24,24) #101
+RHKEcount2[2,] <- c(4,4) #clarence
+RHKEcount2[3,] <- c(6,6) #coffee
+RHKEcount2[4,] <- c(17,17) #crazy
+RHKEcount2[5,] <- c(9,9) #edwards
+RHKEcount2[6,] <- c(16,16) #elochoman
+RHKEcount2[7,] <- c(13,13) #falk
+RHKEcount2[8,] <- c(29,29) #fern
+RHKEcount2[9,] <- c(53,53) #holcombe
+RHKEcount2[10,] <- c(11,11) #nehalem
+RHKEcount2[11,] <- c(11,11) #rock
+RHKEcount2[12,] <- c(26,26) #wlest
+# AK changed
+RHCA_list <- list(
+pcr_n = qPCR.N,  # Number of eDNA qPCR replicates per site
+pcr_k = RHCAqPCR.K, # Number of positive eDNA qPCR detections per site
+count = RHCAcount  # Traditional survey count data
+)
+# AK changed
+RHKE_list <- list(
+pcr_n = qPCR.N,  # Number of eDNA qPCR replicates per site
+pcr_k = RHKEqPCR.K, # Number of positive eDNA qPCR detections per site
+count = RHKEcount  # Traditional survey count data
+)
+RHKE.fit1 <- joint_model(data = RHKE_list, family = 'poisson',
+p10_priors = c(1, 1000), # new prior
+q=FALSE,
+n_warmup = 1000, n_iter = 8000,
+adapt_delta = 0.9999,
+multicore = TRUE # this argument runs the chains in parallel (i.e., faster)
+)
+roxygen2::roxygenise()
+lintr::lint_package()
+lintr::lint_package()
+install.packages('bayesplot')
diff --git a/inst/stan/joint_continuous.stan b/inst/stan/joint_continuous.stan
index bcaeca5..84e79ee 100644
--- a/inst/stan/joint_continuous.stan
+++ b/inst/stan/joint_continuous.stan
@@ -96,7 +96,6 @@ transformed parameters {
 }
 
 model {
-
   // get lambda
   real lambda[n_C];
   lambda = get_lambda_continuous(ctch, coef, mat, alpha_gamma, R_ind, n_C);
@@ -115,13 +114,11 @@ model {
     }
   }
 
-
   //priors
   log_p10 ~ normal(p10_priors[1], p10_priors[2]); // p10 prior
   alpha ~ normal(alphapriors[1], alphapriors[2]); // sitecov shrinkage priors
   beta_gamma ~ gamma(bgammapriors[1], bgammapriors[2]);
   alpha_gamma ~ gamma(agammapriors[1], agammapriors[2]);
-
 }
 
 generated quantities {
diff --git a/inst/stan/joint_count.stan b/inst/stan/joint_count.stan
index 8f6e9d4..8911204 100644
--- a/inst/stan/joint_count.stan
+++ b/inst/stan/joint_count.stan
@@ -87,7 +87,6 @@ transformed parameters {
 }
 
 model {
-
   // get lambda
   real lambda[n_C];
   lambda = get_lambda_count(ctch, coef, mat, mu_trad, R_ind, n_C);
@@ -118,7 +117,6 @@ model {
   if (negbin == 1) {
     phi ~ gamma(phi_priors[1], phi_priors[2]); // phi prior
   }
-
 }
 
 generated quantities {
@@ -144,7 +142,6 @@ generated quantities {
     Nloc_dna, Nloc_trad, p_dna, p10, mat_site, alpha
   );
 
-
   ////////////////////////////////
   // get point-wise log likelihood
 
@@ -153,5 +150,4 @@ generated quantities {
     n_E, n_K, n_N, p_trad, L_ind, n_C, n_S, S_dna,
     Nloc_dna, K_dna, N_dna, p_dna, L_dna
   );
-
 }
diff --git a/inst/stan/traditional_continuous.stan b/inst/stan/traditional_continuous.stan
index 018776e..66945ff 100644
--- a/inst/stan/traditional_continuous.stan
+++ b/inst/stan/traditional_continuous.stan
@@ -49,7 +49,6 @@ transformed parameters {
 }
 
 model {
-
   // get lambda
   real lambda[n_C];
   lambda = get_lambda_continuous(ctch, coef, mat, alpha, R_ind, n_C);
@@ -60,7 +59,6 @@ model {
 
   alpha ~ gamma(alphapriors[1], alphapriors[2]);
   beta ~ gamma(betapriors[1], betapriors[2]);
-
 }
 
 generated quantities{
@@ -83,5 +81,4 @@ generated quantities{
   log_lik = calc_loglik_tradmod_continuous(
     beta, E_trans, R_ind, n_C, ctch, coef, mat, alpha
   );
-
 }
diff --git a/inst/stan/traditional_count.stan b/inst/stan/traditional_count.stan
index 8fe6b6c..3696d0f 100644
--- a/inst/stan/traditional_count.stan
+++ b/inst/stan/traditional_count.stan
@@ -43,7 +43,6 @@ transformed parameters {
 }
 
 model {
-
   // get lambda
   real lambda[n_C];
   lambda = get_lambda_count(ctch, coef, mat, mu_1, R_ind, n_C);
@@ -61,7 +60,6 @@ model {
   if (negbin == 1) {
     phi ~ gamma(phi_priors[1], phi_priors[2]); // phi prior
   }
-
 }
 
 generated quantities{
@@ -83,5 +81,4 @@ generated quantities{
   log_lik = calc_loglik_tradmod_count(
     negbin, phi, n_E, n_C, ctch, coef, mat, mu_1, R_ind
   );
-
 }