Skip to content

Commit

Permalink
collapsed count stan files
Browse files Browse the repository at this point in the history
  • Loading branch information
@abigailkeller
abigailkeller committed Dec 23, 2024
1 parent add447b commit 9d73549
Show file tree
Hide file tree
Showing 14 changed files with 169 additions and 379 deletions.
13 changes: 7 additions & 6 deletions R/detectionCalculate.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -169,8 +169,9 @@ isCatch <- function(pars){
out <- ifelse('q' %in% pars,TRUE,FALSE)
return(out)
}
isNegbin <- function(pars){
out <- ifelse('phi' %in% pars,TRUE,FALSE)
isNegbin <- function(modelfit){
out <- ifelse(modelfit@par_dims$phi == 1,
TRUE,FALSE)
return(out)
}

Expand All @@ -185,8 +186,8 @@ get_ntrad <- function(pars, modelfit, mu, probability){
ntrad <- seq(from = 0, to = 50000, by = 1)

# P(X = 0 | mu) in one traditional survey trial
if(isNegbin(pars)){
phi <- stats::median(unlist(rstan::extract(modelfit, pars = 'phi')))
if(isNegbin(modelfit)){
phi <- stats::median(unlist(rstan::extract(modelfit, pars = 'phi[1]')))
pr <- stats::pnbinom(q = 0, mu = mu[i], size = phi)
} else {
pr <- stats::ppois(q = 0, lambda = mu[i])
Expand Down Expand Up @@ -227,8 +228,8 @@ get_ntrad_q <- function(pars, modelfit, mu, probability){
ntrad <- seq(from = 0, to = 50000, by = 1)

# P(X = 0 | mu) in one traditional survey trial
if(isNegbin(pars)){
phi <- stats::median(unlist(rstan::extract(modelfit, pars = 'phi')))
if(isNegbin(modelfit)){
phi <- stats::median(unlist(rstan::extract(modelfit, pars = 'phi[1]')))
pr <- stats::pnbinom(q = 0, mu = q_list[j]*mu[i], size = phi)
} else {
pr <- stats::ppois(q = 0, lambda = q_list[j]*mu[i])
Expand Down
27 changes: 18 additions & 9 deletions R/jointModel.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -381,10 +381,17 @@ jointModel <- function(data, cov = NULL, family = 'poisson',
#)
#}
# append data if family == negbin
if(isNegbin_type(family)){
if(get_family_index(family)==2){
model_data <- rlist::list.append(
model_data,
phipriors = phipriors
phipriors = phipriors,
negbin = 1
)
} else if(get_family_index(family)==1){
model_data <- rlist::list.append(
model_data,
phipriors = c(1,1),
negbin = 0
)
}

Expand Down Expand Up @@ -412,11 +419,9 @@ jointModel <- function(data, cov = NULL, family = 'poisson',

# run model
out <- rstan::sampling(
c(stanmodels$joint_binary_cov_catchability_pois,
stanmodels$joint_binary_cov_catchability_negbin,
c(stanmodels$joint_binary_cov_catchability_count,
stanmodels$joint_binary_cov_catchability_gamma,
stanmodels$joint_binary_cov_pois,
stanmodels$joint_binary_cov_negbin,
stanmodels$joint_binary_cov_count,
stanmodels$joint_binary_cov_gamma)[model_index][[1]],
data = model_data,
cores = cores,
Expand Down Expand Up @@ -502,10 +507,14 @@ get_stan_model <- function(q, family){

index <- get_family_index(family)

if(isCatch_type(q)){
final_index <- index
if(isCatch_type(q) && index %in% c(1,2)){
final_index <- 1
} else if(isCatch_type(q) && index == 3){
final_index <- 2
} else if(!isCatch_type(q) && index %in% c(1,2)){
final_index <- 3
} else {
final_index <- index + 3
final_index <- 4
}
return(final_index)
}
Expand Down
12 changes: 6 additions & 6 deletions R/jointSummarize.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -84,7 +84,7 @@ jointSummarize <- function(modelfit, par = 'all', probs = c(0.025,0.975),
# get summary
if(all(par == 'all')){

params <- get_all_params(modelfit@model_pars)
params <- get_all_params(modelfit@model_pars,modelfit)
#' @srrstats {G2.4,G2.4a} explicit conversion to integers for sampling
#' arguments
out <- round(rstan::summary(modelfit, pars = params, probs = probs,
Expand All @@ -109,13 +109,13 @@ isCatch <- function(pars){
out <- ifelse('q' %in% pars,TRUE,FALSE)
return(out)
}
isNegbin <- function(pars){
out <- ifelse('phi' %in% pars,TRUE,FALSE)
isNegbin <- function(modelfit){
out <- ifelse(modelfit@par_dims$phi == 1,TRUE,FALSE)
return(out)
}

# function to get vector of all param names
get_all_params <- function(pars){
get_all_params <- function(pars,modelfit){
params <- c('mu')

# catchability
Expand All @@ -129,8 +129,8 @@ get_all_params <- function(pars){
}

# negbin
if(isNegbin(pars)){
params <- c(params,'phi')
if(isNegbin(modelfit)){
params <- c(params,'phi[1]')
}

return(params)
Expand Down
14 changes: 10 additions & 4 deletions R/traditionalModel.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -217,7 +217,7 @@ traditionalModel <- function(data, family = 'poisson',
inits <- init_trad_catchability(n.chain, count_all, q_names,
initial_values)
out <- rstan::sampling(c(
stanmodels$traditional_catchability_pois,
stanmodels$traditional_catchability_count,
stanmodels$traditional_catchability_gamma)[model_index][[1]],
data = list(
Nloc = length(unique(count_all$L)),
Expand All @@ -226,6 +226,8 @@ traditionalModel <- function(data, family = 'poisson',
E = count_all$count,
nparams = length(q_names),
mat = as.integer(count_all$count.type),
phipriors = c(1,1),
negbin = 0,
control = list(adapt_delta = adapt_delta)
),
cores = cores,
Expand All @@ -245,7 +247,7 @@ traditionalModel <- function(data, family = 'poisson',
##run model, catchability, negbin
inits <- init_trad_catchability(n.chain, count_all, q_names,
initial_values)
out <- rstan::sampling(stanmodels$traditional_catchability_negbin,
out <- rstan::sampling(stanmodels$traditional_catchability_count,
data = list(
Nloc = length(unique(count_all$L)),
C = nrow(count_all),
Expand All @@ -254,6 +256,7 @@ traditionalModel <- function(data, family = 'poisson',
nparams = length(q_names),
mat = as.integer(count_all$count.type),
phipriors = phipriors,
negbin = 1,
control = list(adapt_delta = adapt_delta)
),
cores = cores,
Expand All @@ -274,13 +277,15 @@ traditionalModel <- function(data, family = 'poisson',
model_index <- dplyr::case_when(family == 'poisson'~ 1,
family == 'gamma' ~ 2)
inits <- init_trad(n.chain, count_all, initial_values)
out <- rstan::sampling(c(stanmodels$traditional_pois,
out <- rstan::sampling(c(stanmodels$traditional_count,
stanmodels$traditional_gamma)[model_index][[1]],
data = list(
Nloc = length(unique(count_all$L)),
C = nrow(count_all),
R = count_all$L,
E = count_all$count,
phipriors = c(1,1),
negbin = 0,
control = list(adapt_delta = adapt_delta,
stepsize = 0.5)
),
Expand All @@ -300,13 +305,14 @@ traditionalModel <- function(data, family = 'poisson',
} else if(q == FALSE && family == 'negbin'){
##run model, no catchability, negbin
inits <- init_trad(n.chain, count_all, initial_values)
out <- rstan::sampling(stanmodels$traditional_negbin,
out <- rstan::sampling(stanmodels$traditional_count,
data = list(
Nloc = length(unique(count_all$L)),
C = nrow(count_all),
R = count_all$L,
E = count_all$count,
phipriors = phipriors,
negbin = 1,
control = list(adapt_delta = adapt_delta)
),
cores = cores,
Expand Down
28 changes: 20 additions & 8 deletions ...joint_binary_cov_catchability_negbin.stan → .../joint_binary_cov_catchability_count.stan
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,7 @@ data{/////////////////////////////////////////////////////////////////////
int<lower=0> nsitecov; // number of site-level covariates
matrix[Nloc_trad+Nloc_dna,nsitecov] mat_site; // matrix of site-level covariates
array[2] real phipriors; // priors for gamma distrib on phi
int<lower=0,upper=1> negbin; // binary indicator of negative binomial

}

Expand All @@ -29,7 +30,7 @@ parameters{/////////////////////////////////////////////////////////////////////
array[Nloc_dna] real<lower=0, upper = 1> p_dna; // total detection probability
vector<lower=-0.99999>[nparams] q_trans; // catchability coefficients
vector[nsitecov] alpha; // site-level beta covariates
real<lower=0> phi; // dispersion parameter
real<lower=0> phi[(negbin == 1) ? 1 : 0]; // dispersion parameter
}

transformed parameters{/////////////////////////////////////////////////////////////////////
Expand All @@ -46,9 +47,14 @@ transformed parameters{/////////////////////////////////////////////////////////
model{/////////////////////////////////////////////////////////////////////


for(j in 1:C){
E[j] ~ neg_binomial_2(coef[mat[j]]*mu_trad_1[R[j]], phi); // Eq. 1.1
}
if(negbin == 1)
for(j in 1:C){
E[j] ~ neg_binomial_2(coef[mat[j]]*mu_trad_1[R[j]], phi); // Eq. 1.1
}
if(negbin == 0)
for(j in 1:C){
E[j] ~ poisson(coef[mat[j]]*mu_trad_1[R[j]]); // Eq. 1.1
}

for (i in 1:S){
K[i] ~ binomial(N[i], p_trad[L[i]]); // Eq. 1.4
Expand All @@ -63,7 +69,8 @@ model{/////////////////////////////////////////////////////////////////////
//priors
log_p10 ~ normal(p10priors[1], p10priors[2]); // p10 prior
alpha ~ normal(0,10); // sitecov shrinkage priors
phi ~ gamma(phipriors[1], phipriors[2]); // phi prior
if(negbin == 1)
phi ~ gamma(phipriors[1], phipriors[2]); // phi prior

}

Expand Down Expand Up @@ -97,9 +104,14 @@ generated quantities{
////////////////////////////////
// get point-wise log likelihood

for(j in 1:C){
log_lik[j] = neg_binomial_2_lpmf(E[j] | coef[mat[j]]*mu_trad_1[R[j]], phi); //store log likelihood of traditional data given model
}
if(negbin == 1)
for(j in 1:C){
log_lik[j] = neg_binomial_2_lpmf(E[j] | coef[mat[j]]*mu_trad_1[R[j]], phi); //store log likelihood of traditional data given model
}
if(negbin == 0)
for(j in 1:C){
log_lik[j] = poisson_lpmf(E[j] | coef[mat[j]]*mu_trad_1[R[j]]); //store log likelihood of traditional data given model
}

for(i in 1:S){
log_lik[C+i] = binomial_lpmf(K[i] | N[i], p_trad[L[i]]); //store log likelihood of eDNA data given model
Expand Down
114 changes: 0 additions & 114 deletions inst/stan/joint_binary_cov_catchability_pois.stan
View file

This file was deleted.

Loading

0 comments on commit 9d73549

Please sign in to comment.