finished adding functions

inst/stan/functions/calc_loglik.stan

@@ -1,5 +1,6 @@
 /* functions for calculating log likelihood */
 
+// calculate log likelihood of eDNA data
 vector calc_loglik_dna(
   int S,
   int S_dna,
@@ -28,6 +29,7 @@ vector calc_loglik_dna(
     return log_lik;
   }
 
+// calculate log likelihood of traditional count data
 vector calc_loglik_trad_count(
   real[] lambda,
   int negbin,
@@ -50,6 +52,7 @@ vector calc_loglik_trad_count(
     return log_lik;
   }
 
+// calculate log likelihood of traditional continuous data
 vector calc_loglik_trad_continuous(
   real[] lambda,
   vector beta_gamma,
@@ -67,7 +70,7 @@ vector calc_loglik_trad_continuous(
     return log_lik;
   }
 
-
+// calculate log likelihood of data in joint count model
 vector calc_loglik_count(
   int ctch,
   real[] coef,
@@ -108,6 +111,7 @@ vector calc_loglik_count(
 
   }
 
+// calculate log likelihood of data in joint continuous model
 vector calc_loglik_continuous(
   int ctch,
   real[] coef,
@@ -147,6 +151,7 @@ vector calc_loglik_continuous(
 
   }
 
+// calculate lambda for count data
 real[] get_lambda_count(
   int ctch,
   real[] coef,
@@ -164,6 +169,7 @@ real[] get_lambda_count(
     return lambda;
 }
 
+// calculate lambda for continuous data
 real[] get_lambda_continuous(
   int ctch,
   real[] coef,
@@ -181,3 +187,50 @@ real[] get_lambda_continuous(
     return lambda;
 }
 
+// calculate log likelihood of data in traditional count model
+vector calc_loglik_tradmod_count(
+  int negbin,
+  real[] phi,
+  int[] E,
+  int C,
+  int ctch,
+  real[] coef,
+  int[] mat,
+  vector mu_1,
+  int[] R){
+
+    vector[C] log_lik;
+
+    //get lambda
+    real lambda[C];
+    lambda = get_lambda_count(ctch, coef, mat, mu_1, R, C);
+
+    //store log likelihood of traditional data given model
+    log_lik = calc_loglik_trad_count(lambda, negbin, phi, E, C);
+
+    return log_lik;
+  }
+
+// calculate log likelihood of data in traditional continuous model
+vector calc_loglik_tradmod_continuous(
+  vector beta,
+  real[] E_trans,
+  int[] R,
+  int C,
+  int ctch,
+  real[] coef,
+  int[] mat,
+  vector alpha){
+
+    vector[C] log_lik;
+
+    //get lambda
+    real lambda[C];
+    lambda = get_lambda_continuous(ctch, coef, mat, alpha, R, C);
+
+    //store log likelihood of traditional data given model
+    log_lik = calc_loglik_trad_continuous(lambda, beta, E_trans, R, C);
+
+    return log_lik;
+  }
+

inst/stan/functions/calc_mu.stan

@@ -0,0 +1,80 @@
+/* functions for creating mu as a generated quantity */
+
+// calculate mu for joint models
+matrix calc_mu(
+  int[] trad_ind,
+  int[] dna_ind,
+  vector mu_trad,
+  int ctch,
+  int nparams,
+  vector q,
+  int Nloc_dna,
+  int Nloc_trad,
+  real[] p_dna,
+  real p10,
+  matrix mat_site,
+  vector alpha){
+
+    matrix[Nloc_dna+Nloc_trad,nparams+1] mu;
+
+    mu[trad_ind, 1] = mu_trad;
+    if(ctch == 1)
+      mu[trad_ind, 2:(nparams + 1)] = mu_trad * q';
+
+    if(Nloc_dna > 0)
+     for (i in 1:Nloc_dna){
+       real p11_dna[Nloc_dna];
+       p11_dna[i] = p_dna[i] - p10;
+       mu[dna_ind[i],1] = p11_dna[i]*exp(dot_product(to_vector(mat_site[dna_ind[i]]),alpha))/(1-p11_dna[i]);
+     }
+     if(ctch == 1)
+       for (i in 1:Nloc_dna){
+         mu[dna_ind[i], 2:(nparams + 1)] = mu[dna_ind[i], 1] * q';
+       }
+
+    return mu;
+  }
+
+
+// calculate mu for traditional count model
+matrix calc_mu_trad_count(
+  int Nloc,
+  int nparams,
+  vector mu_1,
+  vector q,
+  int ctch){
+
+    matrix[Nloc,nparams+1] mu;
+
+    mu[, 1] = mu_1;
+
+    if(ctch == 1){
+      mu[, 2:(nparams + 1)] = mu_1 * q';
+    }
+
+    return mu;
+  }
+
+// calculate mu for traditional continuous model
+matrix calc_mu_trad_continuous(
+  int Nloc,
+  int nparams,
+  vector alpha,
+  vector beta,
+  vector q,
+  int ctch){
+
+    matrix[Nloc,nparams+1] mu;
+
+    for(j in 1:Nloc){
+      mu[j,1] = alpha[j]/beta[j];
+    }
+
+    if(ctch == 1)
+      for(i in 1:nparams){
+        mu[,i+1] = to_vector(mu[,1])*q[i];
+      }
+
+    return mu;
+  }
+

inst/stan/functions/calc_p11.stan

@@ -0,0 +1,18 @@
+/* function for calculating p11 */
+
+// calculate p11 for joint models
+vector calc_p11(
+  int Nloc_trad,
+  vector mu_trad,
+  matrix mat_site,
+  int[] trad_ind,
+  vector alpha){
+
+    vector[Nloc_trad] p11_trad;
+
+    p11_trad = mu_trad ./ (mu_trad + exp(mat_site[trad_ind, ] * alpha)); // Eq. 1.2
+
+    return p11_trad;
+  }
+
+

inst/stan/joint_continuous.stan

@@ -1,5 +1,7 @@
 functions {
   #include /functions/calc_loglik.stan
+  #include /functions/calc_mu.stan
+  #include /functions/calc_p11.stan
 }
 
 
@@ -45,7 +47,7 @@ transformed parameters{/////////////////////////////////////////////////////////
   array[C] real<lower=0> E_trans;
 
   mu_trad = alpha_gamma ./ beta_gamma;
-  p11_trad = mu_trad ./ (mu_trad + exp(mat_site[trad_ind, ] * alpha)); // Eq. 1.2
+  p11_trad = calc_p11(Nloc_trad, mu_trad, mat_site, trad_ind, alpha); // Eq. 1.2
   p_trad = p11_trad + exp(log_p10); // Eq. 1.3
 
   if(ctch == 1)
@@ -90,30 +92,20 @@ generated quantities{
   vector[C+S+S_dna] log_lik;
   real p10;
   matrix[Nloc_dna+Nloc_trad,nparams+1] mu;  // matrix of catch rates
-  array[Nloc_dna] real<lower=0, upper = 1> p11_dna; // true-positive detection probability
   vector[Nloc_trad] beta;
 
   ////////////////////////////////////
   // transform to interpretable params
   p10 = exp(log_p10);
 
-  mu[trad_ind, 1] = mu_trad;
   if(ctch == 1)
     q = q_trans + 1;
-    mu[trad_ind, 2:(nparams + 1)] = mu_trad * q';
 
   beta = mat_site[trad_ind] * alpha;
 
-
-  if(Nloc_dna > 0)
-     for (i in 1:Nloc_dna){
-       p11_dna[i] = p_dna[i] - p10;
-       mu[dna_ind[i],1] = p11_dna[i]*exp(dot_product(mat_site[dna_ind[i]],alpha))/(1-p11_dna[i]);
-     }
-     if(ctch == 1)
-       for (i in 1:Nloc_dna){
-         mu[dna_ind[i], 2:(nparams + 1)] = mu[dna_ind[i], 1] * q';
-       }
+  mu = calc_mu(trad_ind, dna_ind, mu_trad, ctch, nparams, q,
+               Nloc_dna, Nloc_trad, p_dna, p10,
+               mat_site, alpha);
 
   ////////////////////////////////
   // get point-wise log likelihood

inst/stan/joint_count.stan

@@ -1,5 +1,7 @@
 functions {
   #include /functions/calc_loglik.stan
+  #include /functions/calc_mu.stan
+  #include /functions/calc_p11.stan
 }
 
 data{/////////////////////////////////////////////////////////////////////
@@ -43,7 +45,7 @@ transformed parameters{/////////////////////////////////////////////////////////
   vector<lower=0, upper = 1>[Nloc_trad] p_trad;   // total detection probability
   real<lower=0> coef[(ctch == 1) ? nparams+1 :  0];
 
-  p11_trad = mu_trad ./ (mu_trad + exp(mat_site[trad_ind, ] * alpha)); // Eq. 1.2
+  p11_trad = calc_p11(Nloc_trad, mu_trad, mat_site, trad_ind, alpha); // Eq. 1.2
   p_trad = p11_trad + exp(log_p10); // Eq. 1.3
 
   if(ctch == 1)
@@ -91,30 +93,21 @@ generated quantities{
   vector[C+S+S_dna] log_lik;
   real p10;
   matrix[Nloc_dna+Nloc_trad,nparams+1] mu;  // matrix of expected catch rates
-  array[Nloc_dna] real<lower=0, upper = 1> p11_dna; // true-positive detection probability
   vector[Nloc_trad] beta;
 
   ////////////////////////////////////
   // transform to interpretable params
   p10 = exp(log_p10);
 
-  mu[trad_ind, 1] = mu_trad;
   if(ctch == 1)
     q = q_trans + 1;
-    mu[trad_ind, 2:(nparams + 1)] = mu_trad * q';
 
   beta = mat_site[trad_ind] * alpha;
 
+  mu = calc_mu(trad_ind, dna_ind, mu_trad, ctch, nparams, q,
+               Nloc_dna, Nloc_trad, p_dna, p10,
+               mat_site, alpha);
 
-  if(Nloc_dna > 0)
-     for (i in 1:Nloc_dna){
-       p11_dna[i] = p_dna[i] - p10;
-       mu[dna_ind[i],1] = p11_dna[i]*exp(dot_product(mat_site[dna_ind[i]],alpha))/(1-p11_dna[i]);
-     }
-     if(ctch == 1)
-       for (i in 1:Nloc_dna){
-         mu[dna_ind[i], 2:(nparams + 1)] = mu[dna_ind[i], 1] * q';
-       }
 
   ////////////////////////////////
   // get point-wise log likelihood

inst/stan/traditional_continuous.stan

@@ -1,5 +1,6 @@
 functions {
   #include /functions/calc_loglik.stan
+  #include /functions/calc_mu.stan
 }
 
 data{/////////////////////////////////////////////////////////////////////
@@ -54,23 +55,17 @@ generated quantities{
 
   ////////////////////////////////////
   // transform to interpretable params
-  for(j in 1:Nloc){
-    mu[j,1] = alpha[j]/beta[j];
-  }
-
   if(ctch == 1)
     q = q_trans + 1;
-    for(i in 1:nparams){
-      mu[,i+1] = to_vector(mu[,1])*q[i];
-    }
+
+  mu = calc_mu_trad_continuous(Nloc, nparams, alpha, beta, q, ctch);
 
   ////////////////////////////////
   // get point-wise log likelihood
 
   // get lambda
-  real lambda[C];
-  lambda = get_lambda_continuous(ctch, coef, mat, alpha, R, C);
-  log_lik = calc_loglik_trad_continuous(lambda, beta, E_trans, R, C);
+  log_lik = calc_loglik_tradmod_continuous(beta, E_trans, R, C,
+                                           ctch, coef, mat, alpha);
 
 }
 

inst/stan/traditional_count.stan

@@ -1,5 +1,6 @@
 functions {
   #include /functions/calc_loglik.stan
+  #include /functions/calc_mu.stan
 }
 
 data{/////////////////////////////////////////////////////////////////////
@@ -58,18 +59,16 @@ generated quantities{
 
   ////////////////////////////////////
   // transform to interpretable params
-  mu[, 1] = to_vector(mu_1);
 
   if(ctch == 1)
     q = q_trans + 1;
-    mu[, 2:(nparams + 1)] = to_vector(mu_1) * q';
+
+  mu = calc_mu_trad_count(Nloc, nparams, mu_1, q, ctch);
 
   ////////////////////////////////
   // get point-wise log likelihood
-
-  real lambda[C];
-  lambda = get_lambda_count(ctch, coef, mat, mu_1, R, C);
-  log_lik = calc_loglik_trad_count(lambda, negbin, phi, E, C);
+  log_lik = calc_loglik_tradmod_count(negbin, phi, E, C, ctch,
+                                      coef, mat, mu_1, R);
 
 }