diff --git a/DESCRIPTION b/DESCRIPTION index 485830b..cf357a8 100644 --- a/DESCRIPTION +++ b/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( diff --git a/NEWS.md b/NEWS.md index 946da7e..66b1d81 100644 --- a/NEWS.md +++ b/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. diff --git a/inst/stan/COA_Standard_gaussian.stan b/inst/stan/COA_Standard_gaussian.stan index adcf85a..d18ebe1 100644 --- a/inst/stan/COA_Standard_gaussian.stan +++ b/inst/stan/COA_Standard_gaussian.stan @@ -1,41 +1,22 @@ // Declare data data { - int nind; // number of individuals - int nrec; // number of receivers - int ntime; // number of time steps - int ntrans; // number of trials/expected number of transmissions per time step - array[nind, ntime, nrec] int 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 alpha0; // detection probability intercept on the logit scale - bounds are to ensure only searching reasonable parameter space - real alpha1; // coef. for decline in detection probability with distance - - // latent variables - matrix [nind, ntime] sx; // E-W center of activity coordinate - bounds reflect spatial extent - matrix [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 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 { diff --git a/inst/stan/COA_Standard_logistic.stan b/inst/stan/COA_Standard_logistic.stan index 6ec29be..357d9f3 100644 --- a/inst/stan/COA_Standard_logistic.stan +++ b/inst/stan/COA_Standard_logistic.stan @@ -1,41 +1,19 @@ // Declare data data { - int nind; // number of individuals - int nrec; // number of receivers - int ntime; // number of time steps - int ntrans; // number of trials/expected number of transmissions per time step - array[nind, ntime, nrec] int 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 alpha0; // detection probability intercept on the logit scale - bounds are to ensure only searching reasonable parameter space - real alpha1; // coef. for decline in detection probability with distance - - // latent variables - matrix [nind, ntime] sx; // E-W center of activity coordinate - bounds reflect spatial extent - matrix [nind, ntime] sy; // N-S center of activity coordinate - bounds reflect spatial extent + real 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 { diff --git a/inst/stan/COA_Tag_Integrated_gaussian.stan b/inst/stan/COA_Tag_Integrated_gaussian.stan index ef709fc..b192839 100644 --- a/inst/stan/COA_Tag_Integrated_gaussian.stan +++ b/inst/stan/COA_Tag_Integrated_gaussian.stan @@ -1,26 +1,16 @@ // Declare data data { - int nind; // number of individuals - int nrec; // number of receivers - int ntime; // number of time steps + #include include/shared_data.stan int ntest; // number of test tags - - // number of trials/expected number of transmissions per time step - int ntrans; - // number of detections for each individual at each receiver in each time step - array[nind, ntime, nrec] int y; // number of detections from each test tag at each receiver in each time step array[ntest, ntime, nrec] int 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[ntime, nrec] alpha0; // time effect - real alpha1; // coef. for decline in detection probability with distance - - // latent variables - // E-W center of activity coordinate - bounds reflect spatial extent - matrix[nind, ntime] sx; - // N-S center of activity coordinate - bounds reflect spatial extent - matrix[nind, ntime] sy; + matrix[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 { diff --git a/inst/stan/COA_Tag_Integrated_logistic.stan b/inst/stan/COA_Tag_Integrated_logistic.stan index d339fe8..3d8d968 100644 --- a/inst/stan/COA_Tag_Integrated_logistic.stan +++ b/inst/stan/COA_Tag_Integrated_logistic.stan @@ -1,26 +1,15 @@ // Declare data data { - int nind; // number of individuals - int nrec; // number of receivers - int ntime; // number of time steps + #include include/shared_data.stan int ntest; // number of test tags - - // number of trials/expected number of transmissions per time step - int ntrans; - // number of detections for each individual at each receiver in each time step - array[nind, ntime, nrec] int y; // number of detections from each test tag at each receiver in each time step array[ntest, ntime, nrec] int 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[ntime, nrec] alpha0; // time effect - real alpha1; // coef. for decline in detection probability with distance - - // latent variables - // E-W center of activity coordinate - bounds reflect spatial extent - matrix[nind, ntime] sx; - // N-S center of activity coordinate - bounds reflect spatial extent - matrix[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 { diff --git a/inst/stan/COA_TimeVarying_gaussian.stan b/inst/stan/COA_TimeVarying_gaussian.stan index d138cf3..3313124 100644 --- a/inst/stan/COA_TimeVarying_gaussian.stan +++ b/inst/stan/COA_TimeVarying_gaussian.stan @@ -1,46 +1,21 @@ // Declare data data { - int nind; // number of individuals - int nrec; // number of receivers - int ntime; // number of time steps - int ntrans; // number of trials/expected number of transmissions per time step - array[nind, ntime, nrec] int 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[ntime, nrec] alpha0; // time effect - real alpha1; // coef. for decline in detection probability with distance - - // latent variables - // E-W center of activity coordinate - bounds reflect spatial extent - matrix[nind, ntime] sx; - // N-S center of activity coordinate - bounds reflect spatial extent - matrix[nind, ntime] sy; + matrix[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 { diff --git a/inst/stan/COA_TimeVarying_logistic.stan b/inst/stan/COA_TimeVarying_logistic.stan index d893b9c..855b659 100644 --- a/inst/stan/COA_TimeVarying_logistic.stan +++ b/inst/stan/COA_TimeVarying_logistic.stan @@ -1,49 +1,20 @@ -// Declare data data { - int nind; // number of individuals - int nrec; // number of receivers - int ntime; // number of time steps - int ntrans; // number of trials/expected number of transmissions per time step - array[nind, ntime, nrec] int 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 = 1; } // Declare parameters parameters { - // fixed effects - matrix[ntime, nrec] alpha0; // time effect - real alpha1; // coef. for decline in detection probability with distance - - // latent variables - // E-W center of activity coordinate - bounds reflect spatial extent - matrix[nind, ntime] sx; - // N-S center of activity coordinate - bounds reflect spatial extent - matrix[nind, ntime] sy; + matrix[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 (i in 1:nind) { - for (t in 1:ntime) { - // Calculate distance - vector[nrec] d = sqrt(square(recX - sx[i, t]) + square(recY - sy[i, t])); - - // row(alpha0, t) is a row_vector; ' converts to column vector - 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 { // Detection probability at a distance of 0 - time-varying diff --git a/inst/stan/include/likelihood_coa_static.stan b/inst/stan/include/likelihood_coa_static.stan new file mode 100644 index 0000000..7f96a92 --- /dev/null +++ b/inst/stan/include/likelihood_coa_static.stan @@ -0,0 +1,13 @@ +alpha0 ~ cauchy(0, 2.5); +alpha1 ~ cauchy(0, 2.5); + +// Standard Model (alpha0 is static) +for (t in 1:ntime) { + for (i in 1:nind) { + vector[nrec] d2 = square(recX - sx[i, t]) + square(recY - sy[i, t]); + vector[nrec] d = logistic ? sqrt(d2) : d2; + + y[i, t] ~ binomial_logit(ntrans, alpha0 - (alpha1 * d)); + } +} + diff --git a/inst/stan/include/likelihood_coa_time_varying.stan b/inst/stan/include/likelihood_coa_time_varying.stan new file mode 100644 index 0000000..d656893 --- /dev/null +++ b/inst/stan/include/likelihood_coa_time_varying.stan @@ -0,0 +1,15 @@ +to_vector(alpha0) ~ cauchy(0, 2.5); +alpha1 ~ cauchy(0, 2.5); + +// Time-varying Model (alpha0 changes over time and receiver) +for (t in 1:ntime) { + row_vector[nrec] alpha0_t = row(alpha0, t); + + for (i in 1:nind) { + vector[nrec] d2 = square(recX - sx[i, t]) + square(recY - sy[i, t]); + vector[nrec] d = logistic ? sqrt(d2) : d2; + + y[i, t] ~ binomial_logit(ntrans, alpha0_t' - (alpha1 * d)); + } +} + diff --git a/inst/stan/include/shared_data.stan b/inst/stan/include/shared_data.stan new file mode 100644 index 0000000..b6d31e8 --- /dev/null +++ b/inst/stan/include/shared_data.stan @@ -0,0 +1,11 @@ + int nind; // number of individuals + int nrec; // number of receivers + int ntime; // number of time steps + int ntrans; // number of trials/expected number of transmissions per time step + array[nind, ntime, nrec] int 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 + + \ No newline at end of file diff --git a/inst/stan/include/shared_parameters.stan b/inst/stan/include/shared_parameters.stan new file mode 100644 index 0000000..884a263 --- /dev/null +++ b/inst/stan/include/shared_parameters.stan @@ -0,0 +1,10 @@ +// fixed effect +// detection probability intercept - max of ~1. [-7,7] covers 99.9% prob. space +real alpha1; // coef. for decline in detection probability with distance + +// latent variables +// E-W center of activity coordinate - bounds reflect spatial extent +matrix[nind, ntime] sx; +// N-S center of activity coordinate - bounds reflect spatial extent +matrix[nind, ntime] sy; +