Merge pull request #13 from trackyverse/includes

Shift redundant code to #includes files

Author: mhpob

DESCRIPTION

@@ -1,5 +1,5 @@
 Package: TelemetrySpace
-Version: 1.3.0
+Version: 1.3.1
 Title: Spatial point process and random field models for electronic tagging data
 Description: A collection of tools to fit spatial models to several types of telemetry data. Models are provided to account for the detection process when estimating individual centers of activity from acoustic telemetry data and to incorporate data from stationary test transmitters when available. Bayesian versions of models are fitted using Stan (http://mc-stan.org/). Maximum likelihood versions are fitted using Template Model Builder (https://kaskr.github.io/adcomp/index.html).
 Authors@R: c(

NEWS.md

@@ -1,5 +1,9 @@
-# TelemetrySpace 1.3.0
-  - Code refactor, **BREAKING CHANGE**: The structure is now "per individual" (loop) > "per time" (loop) > _vectorized over_ receivers. Because of this, input data must be an array with dimensions 1: individual; 2: time; 3: receiver. Previously it was individual, receiver, time.
+# TelemetrySpace 1.3
+## 1.3.1
+  - Code refactor ([PR #13](https://github.com/trackyverse/TelemetrySpace/pull/13)): leverage [Stan's "Includes" framework](https://mc-stan.org/docs/reference-manual/includes.html) to reduce redundant code across files.
+
+## 1.3.0
+  - Code refactor ([PR #12](https://github.com/trackyverse/TelemetrySpace/pull/12)), **BREAKING CHANGE**: The structure is now "per individual" (loop) > "per time" (loop) > _vectorized over_ receivers. Because of this, input data must be an array with dimensions 1: individual; 2: time; 3: receiver. Previously it was individual, receiver, time.
   - `p0` and `sigma` are generated quantities rather than transformed parameters. This should not affect anything on the user side, but the code runs a bit faster.
   - Calculated distances between receivers and COAs have been moved from the transformed parameters to the model block
   - The model now operates on squared distances rather than Euclidean distance.

inst/stan/COA_Standard_gaussian.stan

@@ -1,41 +1,22 @@
 // Declare data
 data {
-  int<lower = 0> nind;               // number of individuals
-  int<lower = 0> nrec;               // number of receivers
-  int<lower = 0> ntime;              // number of time steps
-  int<lower = 0> ntrans;             // number of trials/expected number of transmissions per time step
-  array[nind, ntime, nrec] int<lower = 0> y; // number of detections for each individual at each receiver in each time step
-  vector[nrec] recX;                // receiver locations in east-west direction
-  vector[nrec] recY;              // receiver locations in north-south direction
-  vector[2] xlim;                    // area bounds east-west
-  vector[2] ylim;                    // area boundes north-south
+  #include include/shared_data.stan
+}
+
+transformed data {
+  int logistic = 0; // Tells include file to use squared distance
 }
 
 // Declare parameters
 parameters {
-  // fixed effects
-  real<lower = -7, upper = 7> alpha0;  // detection probability intercept on the logit scale - bounds are to ensure only searching reasonable parameter space
-  real<lower = 0> alpha1;  // coef. for decline in detection probability with distance
-
-  // latent variables
-  matrix<lower = xlim[1], upper = xlim[2]> [nind, ntime] sx; // E-W center of activity coordinate - bounds reflect spatial extent
-  matrix<lower = ylim[1], upper = ylim[2]> [nind, ntime] sy; // N-S center of activity coordinate - bounds reflect spatial extent
+  // detection probability intercept on the logit scale
+  // bounds are to ensure only searching reasonable parameter space
+  real<lower = -5, upper = 5> alpha0;  
+  #include include/shared_parameters.stan
 }
 
 model {
-  // priors
-  alpha0 ~ cauchy(0, 2.5);
-  alpha1 ~ cauchy(0, 2.5);
-
-  for (i in 1:nind) {
-    for (t in 1:ntime) {
-      // Calculate squared distance (d2)
-      vector[nrec] d2 = square(recX - sx[i, t]) + square(recY - sy[i, t]);
-      
-      // Run binomial on logit scale
-      y[i, t] ~ binomial_logit(ntrans, alpha0 - (alpha1 * d2));
-    }
-  }
+  #include include/likelihood_coa_static.stan
 }  
 
 generated quantities {

inst/stan/COA_Standard_logistic.stan

@@ -1,41 +1,19 @@
 // Declare data
 data {
-  int<lower = 0> nind;               // number of individuals
-  int<lower = 0> nrec;               // number of receivers
-  int<lower = 0> ntime;              // number of time steps
-  int<lower = 0> ntrans;             // number of trials/expected number of transmissions per time step
-  array[nind, ntime, nrec] int<lower = 0> y; // number of detections for each individual at each receiver in each time step
-  vector[nrec] recX;                // receiver locations in east-west direction
-  vector[nrec] recY;              // receiver locations in north-south direction
-  vector[2] xlim;                    // area bounds east-west
-  vector[2] ylim;                    // area boundes north-south
+  #include include/shared_data.stan
+}
+transformed data {
+  int logistic = 1; // Tells include file to use linear distance (sqrt)
 }
-
 // Declare parameters
 parameters {
   // fixed effects
-  real<lower = -7, upper = 7> alpha0;  // detection probability intercept on the logit scale - bounds are to ensure only searching reasonable parameter space
-  real<lower = 0> alpha1;  // coef. for decline in detection probability with distance
-
-  // latent variables
-  matrix<lower = xlim[1], upper = xlim[2]> [nind, ntime] sx; // E-W center of activity coordinate - bounds reflect spatial extent
-  matrix<lower = ylim[1], upper = ylim[2]> [nind, ntime] sy; // N-S center of activity coordinate - bounds reflect spatial extent
+  real<lower = -5, upper = 5> alpha0;  // detection probability intercept on the logit scale - bounds are to ensure only searching reasonable parameter space
+  #include include/shared_parameters.stan
 }
 
 model {
-  // priors
-  alpha0 ~ cauchy(0, 2.5);
-  alpha1 ~ cauchy(0, 2.5);
-
-  for (i in 1:nind) {
-    for (t in 1:ntime) {
-      // Calculate euclidean distance (d)
-      vector[nrec] d = sqrt(square(recX - sx[i, t]) + square(recY - sy[i, t]));
-      
-      // Run binomial on logit scale
-      y[i, t] ~ binomial_logit(ntrans, alpha0 - (alpha1 * d));
-    }
-  }
+  #include include/likelihood_coa_static.stan
 }  
 
 generated quantities {

inst/stan/COA_Tag_Integrated_gaussian.stan

@@ -1,26 +1,16 @@
 // Declare data
 data {
-  int<lower = 0> nind;               // number of individuals
-  int<lower = 0> nrec;               // number of receivers
-  int<lower = 0> ntime;              // number of time steps
+  #include include/shared_data.stan
   int<lower = 0> ntest;              // number of test tags
-
-  // number of trials/expected number of transmissions per time step
-  int<lower = 0> ntrans;
-  // number of detections for each individual at each receiver in each time step
-  array[nind, ntime, nrec] int<lower = 0> y;
   // number of detections from each test tag at each receiver in each time step
   array[ntest, ntime, nrec] int<lower = 0> test;
-
-  vector[nrec] recX; // trap locations in east-west direction
-  vector[nrec] recY; // trap locations in north-south direction
-  vector[2] xlim;  // area bounds east-west
-  vector[2] ylim;                    // area boundes north-south
   vector[ntest] testX;               // test tag locations east-west
   vector[ntest] testY;               // test tag locations north-south
 }
 
 transformed data {
+  int logistic = 0;
+
   // Pre-calculate squared distances from receiver for fixed test tags
   matrix[ntest, nrec] td2;
   for (s in 1:ntest) {
@@ -30,48 +20,26 @@ transformed data {
 
 // Declare parameters
 parameters {
-  // fixed effects
-  // detection probability intercept - max of ~1
-  matrix<lower = -5, upper = 5>[ntime, nrec] alpha0; // time effect
-  real<lower = 0> alpha1;  // coef. for decline in detection probability with distance
-
-  // latent variables
-  // E-W center of activity coordinate - bounds reflect spatial extent
-  matrix<lower = xlim[1], upper = xlim[2]>[nind, ntime] sx;
-  // N-S center of activity coordinate - bounds reflect spatial extent
-  matrix<lower = ylim[1], upper = ylim[2]>[nind, ntime] sy;
+  matrix<lower = -5, upper = 5>[ntime, nrec] alpha0;
+  #include include/shared_parameters.stan
 }
 
 // Model specification
 model {
-  // priors
-  to_vector(alpha0) ~ cauchy(0, 2.5);
-  alpha1 ~ cauchy(0, 2.5);
-
   // Likelihood for test tags (fixed locations)
-  for (s in 1:ntest) { // For each test tag
-    // decay over distance portion of the binomial model
-    vector[nrec] dist_decay = alpha1 * td2[s, ]';
-
-    for (t in 1:ntime) { // Run binomial on logit scale
-      // row(p0, t) pulls the alpha0 vector for all receivers at time t
+  for (t in 1:ntime) {
+    row_vector[nrec] alpha0_t = row(alpha0, t);
+    
+    for (s in 1:ntest) {
+      // decay over distance portion of the binomial model
+      vector[nrec] dist_decay = alpha1 * td2[s, ]';
+      
       test[s, t] ~ binomial_logit(ntrans, row(alpha0, t)' - dist_decay);
     }
   }
 
   // Likelihood for individual COAs (estimated locations)
-  for (i in 1:nind) {
-    for (t in 1:ntime) {
-      // Distance squared from each COA to each receiver at each time
-      vector[nrec] d2 = square(recX - sx[i, t]) + square(recY - sy[i, t]);
-      
-      // Calculate linear predictor
-      vector[nrec] lp = row(alpha0, t)' - (alpha1 * d2);
-      
-      // Run binomial on logit scale
-      y[i, t] ~ binomial_logit(ntrans, lp);
-    }
-  }
+  #include include/likelihood_coa_time_varying.stan
 }  
 
 generated quantities {

inst/stan/COA_Tag_Integrated_logistic.stan

@@ -1,26 +1,15 @@
 // Declare data
 data {
-  int<lower = 0> nind;               // number of individuals
-  int<lower = 0> nrec;               // number of receivers
-  int<lower = 0> ntime;              // number of time steps
+  #include include/shared_data.stan
   int<lower = 0> ntest;              // number of test tags
-
-  // number of trials/expected number of transmissions per time step
-  int<lower = 0> ntrans;
-  // number of detections for each individual at each receiver in each time step
-  array[nind, ntime, nrec] int<lower = 0> y;
   // number of detections from each test tag at each receiver in each time step
   array[ntest, ntime, nrec] int<lower = 0> test;
-
-  vector[nrec] recX; // trap locations in east-west direction
-  vector[nrec] recY; // trap locations in north-south direction
-  vector[2] xlim;  // area bounds east-west
-  vector[2] ylim;                    // area boundes north-south
   vector[ntest] testX;               // test tag locations east-west
   vector[ntest] testY;               // test tag locations north-south
 }
 
 transformed data {
+  int logistic = 1;
   // Pre-calculate squared distances from receiver for fixed test tags
   matrix[ntest, nrec] td;
   for (s in 1:ntest) {
@@ -30,48 +19,27 @@ transformed data {
 
 // Declare parameters
 parameters {
-  // fixed effects
   // detection probability intercept - max of ~1
   matrix<lower = -5, upper = 5>[ntime, nrec] alpha0; // time effect
-  real<lower = 0> alpha1;  // coef. for decline in detection probability with distance
-
-  // latent variables
-  // E-W center of activity coordinate - bounds reflect spatial extent
-  matrix<lower = xlim[1], upper = xlim[2]>[nind, ntime] sx;
-  // N-S center of activity coordinate - bounds reflect spatial extent
-  matrix<lower = ylim[1], upper = ylim[2]>[nind, ntime] sy;
+  #include include/shared_parameters.stan
 }
 
 // Model specification
 model {
-  // priors
-  to_vector(alpha0) ~ cauchy(0, 2.5);
-  alpha1 ~ cauchy(0, 2.5);
-
   // Likelihood for test tags (fixed locations)
-  for (s in 1:ntest) { // For each test tag
-    // decay over distance portion of the binomial model
-    vector[nrec] dist_decay = alpha1 * td[s, ]';
+  for (t in 1:ntime) {
+    row_vector[nrec] alpha0_t = row(alpha0, t);
+    
 
-    for (t in 1:ntime) { // Run binomial on logit scale
-      // row(p0, t) pulls the alpha0 vector for all receivers at time t
+    for (s in 1:ntest) {
+      // decay over distance portion of the binomial model
+      vector[nrec] dist_decay = alpha1 * td[s, ]';
+      
       test[s, t] ~ binomial_logit(ntrans, row(alpha0, t)' - dist_decay);
     }
   }
 
-  // Likelihood for individual COAs (estimated locations)
-  for (i in 1:nind) {
-    for (t in 1:ntime) {
-      // Distance squared from each COA to each receiver at each time
-      vector[nrec] d = sqrt(square(recX - sx[i, t]) + square(recY - sy[i, t]));
-      
-      // Calculate linear predictor
-      vector[nrec] lp = row(alpha0, t)' - (alpha1 * d);
-      
-      // Run binomial on logit scale
-      y[i, t] ~ binomial_logit(ntrans, lp);
-    }
-  }
+  #include include/likelihood_coa_time_varying.stan
 }  
 
 generated quantities {

inst/stan/COA_TimeVarying_gaussian.stan

@@ -1,46 +1,21 @@
 // Declare data
 data {
-  int<lower = 0> nind;               // number of individuals
-  int<lower = 0> nrec;               // number of receivers
-  int<lower = 0> ntime;              // number of time steps
-  int<lower = 0> ntrans;             // number of trials/expected number of transmissions per time step
-  array[nind, ntime, nrec] int<lower = 0> y; // number of detections for each individual at each receiver in each time step
-  vector[nrec] recX;                 // trap locations in east-west direction
-  vector[nrec] recY;                 // trap locations in north-south direction
-  vector[2] xlim;                    // area bounds east-west
-  vector[2] ylim;                    // area boundes north-south
+  #include include/shared_data.stan
+}
+
+transformed data {
+  int logistic = 0;
 }
 
 // Declare parameters
 parameters {
-  // fixed effects
-  matrix<lower = -7, upper = 7>[ntime, nrec] alpha0; // time effect
-  real<lower = 0> alpha1;  // coef. for decline in detection probability with distance
-
-  // latent variables
-  // E-W center of activity coordinate - bounds reflect spatial extent
-  matrix<lower = xlim[1], upper = xlim[2]>[nind, ntime] sx;
-  // N-S center of activity coordinate - bounds reflect spatial extent
-  matrix<lower = ylim[1], upper = ylim[2]>[nind, ntime] sy;
+  matrix<lower = -5, upper = 5>[ntime, nrec] alpha0; // time effect
+  #include include/shared_parameters.stan
 }
 
 // Model specification
 model {
-  // priors
-  to_vector(alpha0) ~ cauchy(0, 2.5);
-  alpha1 ~ cauchy(0, 2.5);
-
-  // likelihood
-  for (i in 1:nind) {
-    for (t in 1:ntime) {
-      // Calculate squared distance
-      vector[nrec] d2 = square(recX - sx[i, t]) + square(recY - sy[i, t]);
-
-      // Run binomial on logit scale
-      // row(alpha0, t) is a row_vector; ' converts to column vector
-      y[i, t] ~ binomial_logit(ntrans, row(alpha0, t)' - (alpha1 * d2));
-    }
-  }
+  #include include/likelihood_coa_time_varying.stan
 } 
 
 generated quantities {