diff --git a/vignettes/_ps_estimate_coa_vignette.Rmd b/vignettes/_estimate_ps_coa_vignette.Rmd similarity index 90% rename from vignettes/_ps_estimate_coa_vignette.Rmd rename to vignettes/_estimate_ps_coa_vignette.Rmd index 2d2d6b4..e297a90 100644 --- a/vignettes/_ps_estimate_coa_vignette.Rmd +++ b/vignettes/_estimate_ps_coa_vignette.Rmd @@ -19,7 +19,7 @@ knitr::opts_chunk$set( ``` ## Introduction -This vignette will walk you through the analyses presented in [Winton et al. 2018](https://besjournals.onlinelibrary.wiley.com/doi/full/10.1111/2041-210X.13080), who describes the use of spatial point process models to estimate individual centers of activity (COA) from passive acoustic telemetry data. This vignette walks through how to prepare the data, using the models, and interpretating the results. We will be using the simplest case, which assumes that detection probabilities/receiver detection ranges remain constant over time, to a more complex application of a test-tag integrated model, that incorporates detection data from one or more stationary test transmitters to estimate time-varying detection ranges. The models are fitted in a Bayesian framework using the Stan software ([Carpenter et al. 2017](https://www.jstatsoft.org/article/view/v076i01/0)); code was modified from models found in [Royle et al. 2013](https://www.sciencedirect.com/book/monograph/9780124059399/spatial-capture-recapture) for fitting spatial point process models to data from camera traps. We prefer the Bayesian approach for COA estimation due to the treatment of uncertainty, but realize the longer computational time required may be prohibitive for some applications. We'd also like to note that the models described can support varying degrees of complexity - not all applications will require (or have the data to support) the most complex version of the model. The simpler the model, the shorter the run-time. +This vignette will walk you through the analyses presented in [Winton et al. 2018](https://besjournals.onlinelibrary.wiley.com/doi/full/10.1111/2041-210X.13080), who describes the use of spatial point process models to estimate individual centers of activity (COA) from passive acoustic telemetry data. We will walk through how to prepare the data, using the models, and interpretating the results. We will be using the simplest case, which assumes that detection probabilities/receiver detection ranges remain constant over time, to a more complex application of a test-tag integrated model, that incorporates detection data from one or more stationary test transmitters to estimate time-varying detection ranges. The models are fitted in a Bayesian framework using the Stan software ([Carpenter et al. 2017](https://www.jstatsoft.org/article/view/v076i01/0)); code was modified from models found in [Royle et al. 2013](https://www.sciencedirect.com/book/monograph/9780124059399/spatial-capture-recapture) for fitting spatial point process models to data from camera traps. We prefer the Bayesian approach for COA estimation due to the treatment of uncertainty, but realize the longer computational time required may be prohibitive for some applications. We'd also like to note that the models described can support varying degrees of complexity - not all applications will require (or have the data to support) the most complex version of the model. The simpler the model, the shorter the run-time. We have tried to make the instructions outlined in this vignette user-friendly since we are a group of applied biologists with varying degrees of statistical experience. If some of the statistical notation outlined here or in [Winton et al. 2018](https://besjournals.onlinelibrary.wiley.com/doi/full/10.1111/2041-210X.13080) remains unclear, feel free to contact us with questions for clarification. This is a new package, so if you find bugs, places where code efficiency could be improved, or instances where the documentation could be made more user-friendly, please let us know through the [issues]( https://github.com/trackyverse/TelemetrySpace) on the GitHub repository! @@ -39,6 +39,7 @@ To run spatial point process/detection probability models in Stan we need to do We have several example datasets along with several functions to assist and streamline the data preparation. First, we need to evaluate the positions of the receivers and create a [Azimuthal Equal Distance projection (aeqd)](https://en.wikipedia.org/wiki/Azimuthal_equidistant_projection). We use an aeqd project set in kilometers (km) for several reasons, the first being that in Stan it is very efficient to calculate distances among receivers, however, the distance between receivers needs to be on the same x and y plane and ranges need to be between 0.2 - 20 km apart. The second is because....(Mike I need you to add stuff here as you know more about aeqd). We will be working with example data for a single Lake Trout (*Salvelinus namaycush*) that was implanted with an acoustic transmitter in Parry Sound which is a large embayment of Georgian Bay, Lake Huron. First, we load all the packages needed to carry out the analysis. + ```{r setup, message=FALSE} { library(bayesplot) @@ -110,7 +111,7 @@ We can notice that this detection data consists of 5 columns with 577 rows. To b For our detection data, we have a few things we need to do, the first is we need to build a time bin that we will create COAs. For this data we are going to use 1 hour but this time bin can range from 30 mins - 1 day or more and depends on the questions you are asking and the species you are working with. -Let's build our time bins +Let's build our time bins. ```{r build time bin} ps_det_example <- build_time_bin(ps_det_example, unit = "1 hour") @@ -168,7 +169,7 @@ There are a few things to know about running a Bayesian analysis, we suggest rea ### Priors -Bayesian analyses rely on supplying uninformed or informed prior distributions for each parameter (coefficient; predictor) in the model. For the puproses of the deteciton probablity model the following priors are assumed and are not adjustable. To learn more about the structure of the priors and likelihoods see [Standard detection probability](https://telemetryspace.trackyverse.org/articles/coa_standard_gaussian_model.html). +Bayesian analyses rely on supplying uninformed or informed prior distributions for each parameter (coefficient; predictor) in the model. For the puproses of the deteciton probablity model the following priors are assumed and are not adjustable. To learn more about the structure of the priors and likelihoods see [standard detection probability model](https://telemetryspace.trackyverse.org/articles/coa_standard_gaussian_model.html). $$ \begin{aligned} @@ -183,7 +184,7 @@ The Cauchy priors on $\alpha_0$ and $\alpha_1$ are weakly informative and regularize the intercept and decay-rate parameters toward zero while allowing heavy tails. Because both parameters carry explicit `lower`/`upper` bounds in the `parameters` block, the priors are implicitly truncated to -those bounds — Stan automatically renormalizes the density over the +those bounds. Stan automatically renormalizes the density over the constrained support, so no separate normalizing constant needs to be added by hand. The activity-center coordinates $s_{x,i,t}$ and $s_{y,i,t}$ have no explicit sampling statement, so they receive Stan's implicit flat (uniform) @@ -224,6 +225,7 @@ m <- COA_Standard( ) ``` +### Convergance and model performance We can inspect the object created with the first object containing the Stan model. To view and work with the model itself call `m$model`. However, as you will see there @@ -233,7 +235,7 @@ Are a bunch of other objects in the objected created. These objects have pulled ```{r summary of model object} summary(m) ``` -Let's look at our trace plots for the model parameters and posterior distributions to ensure the model converged properly. Remember the trace plots should look grassy or caterpillar like. A trace plot is the posterior draw for a given iteration plotted with the iteration number on the x axis and the posterior value for a given parameter or latent variable on the y. We evaluate this for both chains and want to see that both chains are converging on a smiler posterior draw for a given parameter or latent variable. We will first look at parameters of the model. +Let's look at our trace plots for the model parameters and posterior distributions to ensure the model converged properly. Remember the trace plots should look grassy or caterpillar like. A trace plot is the posterior draw for a given iteration plotted with the iteration number on the x-axis and the posterior value for a given parameter or latent variable on the y-axis. We evaluate this for both chains and want to see that both chains are converging on a similar posterior draw for a given parameter or latent variable. We will first look at parameters of the model. ``` {r trace plots and post dist, fig.cap=""} stan_trace(m$model, pars = c("alpha0", "alpha1", "p0", "sigma")) @@ -246,12 +248,11 @@ stan_dens( ) ``` -We can see our trace plots look good and that the posterior distributions of our paramaters look good. +We can see our trace plots look good and that the posterior distributions of our parameters look good. Next let’s look at our latent variables which are `sx` and `sy` or the estimated locations of the fish based on the detection probability. For large/lengthy time periods it will become cumbersome to evaluate the posteriors this way and we recommend inspecting $\hat R$ and ESS. - ``` {r trace plots and post dist latent, fig.cap=""} stan_trace(m$model, pars = c("sx[1,1]", "sx[1,2]", "sy[1,1]", "sy[1,2]")) @@ -262,8 +263,10 @@ stan_dens( linewidth = 0.1 ) ``` + Again, we can see that everything looks good and our model has converged. +### Posteriors Now moving on to the other elements in the outputted object from `COA_*()`. @@ -293,6 +296,7 @@ The fifth element returned, is `data.frame` of the median values for the latent ```{r coas} m$coas ``` + The sixth element returned, is a `data.frame` containing the posterior draws from each non-warm-up iteration from all chains. This contains the posterior distribution for each parameter and latent variable for each individual in each time step. It is unlikely that you will use this object instead we have created further objects that organize this data for specific applications that you are more likely to use. ```{r all draws} @@ -329,7 +333,8 @@ str(yrep) First, we are going to plot the posterior draws for locations. Within the package we have a `sf` object that is the shape of Parry Sound. We can use this to plot the posterior draws of `sx` and `sy` to understand the movement of that of lake trout for 8 hours. -We need to take that `sf` object and transform it into an aeqd projection. +We need to take that `sf` object and transform it into an aeqd projection. + ```{r ps aeqd} ps_aeqd <- ps |> st_transform(aeqd_crs) @@ -402,7 +407,7 @@ p_param <- ggplot( p_param ``` -We can see that they are all quite tightly distributed with `p0` indicating that detection probability at a distance of `0` is between 48 - 54 %, while the `sigma` is around 1 km. +We can see that they are all quite tightly distributed with `p0` indicating that detection probability at a distance of `0` is between 48 - 54 %, while the `sigma` is around 1 km. Considering this is the standard model, `p0` nor other parameters vary. Considering how acoustic telemetry functions, we know that this is an unlikely represenation of these parameters which we can use a time-varying model and a tag-integrated, time-varying model to account for variation in the `p0` and other parameters. Lastly, we can plot the predictive posterior check using `{tidybayes}`. First we need to make detection counts as a vector from our `build_count()` object. @@ -410,7 +415,9 @@ detection counts as a vector from our `build_count()` object. ```{r y_obs} y_obs <- as.vector(ps_count_example[!is.na(ps_count_example)]) ``` + Next we can plot the densities using `ppc_dens_overlay()` from `{bayesplot}`. + ```{r ppc dens, fig.cap=""} ppc <- ppc_dens_overlay(y = y_obs, yrep = yrep$yrep) ppc diff --git a/vignettes/_coa_standard_gaussian_model.Rmd b/vignettes/articles/model_specifications_standard_detection_prob.Rmd similarity index 100% rename from vignettes/_coa_standard_gaussian_model.Rmd rename to vignettes/articles/model_specifications_standard_detection_prob.Rmd diff --git a/vignettes/coa_standard_gaussian_model.Rmd b/vignettes/coa_standard_gaussian_model.Rmd deleted file mode 100644 index fd42eee..0000000 --- a/vignettes/coa_standard_gaussian_model.Rmd +++ /dev/null @@ -1,180 +0,0 @@ ---- -title: "Model specification - Standard Detection Probability" -output: rmarkdown::html_vignette -vignette: > - %\VignetteIndexEntry{Model specification - Standard Detection Probability} - %\VignetteEngine{knitr::rmarkdown} - %\VignetteEncoding{UTF-8} ---- - - - -## Overview - -`COA_Standard_gaussian` is a Bayesian Center of Activity (COA) model that -estimates the latent spatial position of each individual at each discrete -time step from binomial detection counts recorded at an array of fixed -receivers. Detection probability decays with squared distance from the -individual's activity center, which corresponds to a Gaussian (half-normal) -detection function in distance-sampling terms. The "Standard" variant holds -the detection-probability intercept $\alpha_0$ fixed across individuals and -time steps (i.e., a single shared intercept), distinguishing it from -hierarchical or covariate-extended variants of the model. - -## Data - -| Symbol | Stan name | Description | -|---|---|---| -| $N$ | `nind` | Number of individuals | -| $R$ | `nrec` | Number of receivers | -| $T$ | `ntime` | Number of discrete time steps | -| $K$ | `ntrans` | Number of expected transmissions per time step | -| $y_{i,t,r}$ | `y[i, t, r]` | Detections of individual $i$ at receiver $r$ during time step $t$, $y_{i,t,r} \in \{0, 1, \dots, K\}$ | -| $x_r, y_r$ | `recX[r]`, `recY[r]` | Fixed coordinates of receiver $r$ | -| $[x_{\min}, x_{\max}]$ | `xlim` | Spatial bounds of the study area (easting) | -| $[y_{\min}, y_{\max}]$ | `ylim` | Spatial bounds of the study area (northing) | - -for $i = 1, \dots, N$, $t = 1, \dots, T$, $r = 1, \dots, R$. - -## Parameters - -| Symbol | Stan name | Constraint | Description | -|---|---|---|---| -| $\alpha_0$ | `alpha0` | $\alpha_0 \in [-5, 5]$ | Detection-probability intercept (logit scale), shared across all individuals and time steps | -| $\alpha_1$ | `alpha1` | $\alpha_1 > 0$ | Coefficient governing decline in detection probability with squared distance | -| $s_{x,i,t}$ | `sx[i, t]` | $s_{x,i,t} \in [x_{\min}, x_{\max}]$ | Easting coordinate of individual $i$'s activity center at time $t$ | -| $s_{y,i,t}$ | `sy[i, t]` | $s_{y,i,t} \in [y_{\min}, y_{\max}]$ | Northing coordinate of individual $i$'s activity center at time $t$ | - -The activity center for individual $i$ at time $t$ is the latent point -$\mathbf{s}_{i,t} = (s_{x,i,t},\, s_{y,i,t})$. - -## Distance - -For individual $i$ at time $t$, the squared Euclidean distance to receiver -$r$ is - -$$ -d^2_{i,t,r} \;=\; \bigl(x_r - s_{x,i,t}\bigr)^2 \;+\; \bigl(y_r - s_{y,i,t}\bigr)^2 . -$$ - -In the Stan code this is the vector `d2`, computed once per individual–time -combination across all receivers. The `logistic` flag in `transformed data` -is set to `0` for this model, so the working distance term used in the -linear predictor is the squared distance itself, - -$$ -d_{i,t,r} \;=\; d^2_{i,t,r}, -$$ - -rather than its square root. (When `logistic = 1`, an alternate model -variant instead uses $d_{i,t,r} = \sqrt{d^2_{i,t,r}}$; that branch is not -exercised here.) - -## Likelihood - -Detections at each receiver are modeled as binomial counts out of $K$ trials, -with detection probability governed by a logit-linear function of squared -distance from the activity center: - -$$ -y_{i,t,r} \;\sim\; \mathrm{Binomial}\bigl(K,\; p_{i,t,r}\bigr), \qquad -\mathrm{logit}\bigl(p_{i,t,r}\bigr) \;=\; \alpha_0 - \alpha_1 \, d^2_{i,t,r}, -$$ - -equivalently - -$$ -p_{i,t,r} \;=\; \mathrm{logit}^{-1}\!\bigl(\alpha_0 - \alpha_1\, d^2_{i,t,r}\bigr) -\;=\; \frac{1}{1 + \exp\!\bigl[-(\alpha_0 - \alpha_1\, d^2_{i,t,r})\bigr]} . -$$ - -This is implemented via Stan's `binomial_logit`, which parameterizes the -binomial directly on the logit scale for numerical stability: - -$$ -y_{i,t,r} \;\sim\; \mathrm{BinomialLogit}\bigl(K,\; \alpha_0 - \alpha_1\, d^2_{i,t,r}\bigr). -$$ - -The full log-likelihood, summed over all individuals, time steps, and -receivers, is - -$$ -\log L(\alpha_0, \alpha_1, \mathbf{s} \mid \mathbf{y}) -\;=\; \sum_{t=1}^{T} \sum_{i=1}^{N} \sum_{r=1}^{R} -\log \mathrm{BinomialLogit}\bigl(y_{i,t,r} \mid K,\; \alpha_0 - \alpha_1 d^2_{i,t,r}\bigr). -$$ - -## Priors - -$$ -\begin{aligned} -\alpha_0 &\sim \mathrm{Cauchy}(0,\, 2.5), &&\text{truncated to } [-5, 5] \\ -\alpha_1 &\sim \mathrm{Cauchy}(0,\, 2.5), &&\text{truncated to } (0, \infty) \\ -s_{x,i,t} &\sim \mathrm{Uniform}(x_{\min},\, x_{\max}) &&\text{for all } i, t \\ -s_{y,i,t} &\sim \mathrm{Uniform}(y_{\min},\, y_{\max}) &&\text{for all } i, t -\end{aligned} -$$ - -The Cauchy priors on $\alpha_0$ and $\alpha_1$ are weakly informative and -regularize the intercept and decay-rate parameters toward zero while -allowing heavy tails. Because both parameters carry explicit `lower`/`upper` -bounds in the `parameters` block, the priors are implicitly truncated to -those bounds — Stan automatically renormalizes the density over the -constrained support, so no separate normalizing constant needs to be added -by hand. The activity-center coordinates $s_{x,i,t}$ and $s_{y,i,t}$ have no -explicit sampling statement, so they receive Stan's implicit flat (uniform) -prior over their bounded support $[x_{\min}, x_{\max}]$ and -$[y_{\min}, y_{\max}]$, respectively. - -## Joint posterior - -Combining the likelihood and priors, the (unnormalized) joint posterior -density is - -$$ -\begin{aligned} -p(\alpha_0, \alpha_1, \mathbf{s} \mid \mathbf{y}) \;\propto\;\; -& \mathrm{Cauchy}(\alpha_0 \mid 0, 2.5) \cdot \mathrm{Cauchy}(\alpha_1 \mid 0, 2.5) \\[4pt] -\times\;\; & \prod_{t=1}^{T} \prod_{i=1}^{N} \prod_{r=1}^{R} -\mathrm{BinomialLogit}\bigl(y_{i,t,r} \mid K,\; \alpha_0 - \alpha_1 d^2_{i,t,r}(\mathbf{s}_{i,t})\bigr), -\end{aligned} -$$ - -subject to $\alpha_0 \in [-5, 5]$, $\alpha_1 > 0$, -$s_{x,i,t} \in [x_{\min}, x_{\max}]$, and $s_{y,i,t} \in [y_{\min}, y_{\max}]$ -for all $i, t$. - -## Generated quantities - -Two derived quantities are computed from the posterior draws of $\alpha_0$ -and $\alpha_1$: - -**Baseline detection probability.** The detection probability at distance -zero (i.e., when the activity center coincides exactly with a receiver): - -$$ -p_0 \;=\; \mathrm{logit}^{-1}(\alpha_0) \;=\; \frac{1}{1 + e^{-\alpha_0}}. -$$ - -**Detection-decay standard deviation.** The squared-distance decay term -$\alpha_1 d^2$ is the kernel of a Gaussian (half-normal) detection function. -Writing the decay term in the conventional distance-sampling form -$d^2 / (2\sigma^2)$ and matching coefficients, - -$$ -\alpha_1 \;=\; \frac{1}{2\sigma^2} -\quad \Longleftrightarrow \quad -\sigma \;=\; \sqrt{\frac{1}{2\alpha_1}}, -$$ - -so that the detection function can equivalently be written as - -$$ -p_{i,t,r} \;=\; \mathrm{logit}^{-1}\!\left(\alpha_0 - \frac{d^2_{i,t,r}}{2\sigma^2}\right). -$$ - -Here $\sigma$ has units of distance and characterizes the effective spatial -range of detection: it is the standard deviation of the Gaussian decay in -detection log-odds with distance from the activity center, with larger -$\sigma$ corresponding to a more slowly-decaying (longer-range) detection -function. diff --git a/vignettes/ps_estimate_coa_vignette.Rmd b/vignettes/estimate_ps_coa_vignette.Rmd similarity index 82% rename from vignettes/ps_estimate_coa_vignette.Rmd rename to vignettes/estimate_ps_coa_vignette.Rmd index aa0f06f..8516f63 100644 --- a/vignettes/ps_estimate_coa_vignette.Rmd +++ b/vignettes/estimate_ps_coa_vignette.Rmd @@ -12,7 +12,7 @@ vignette: > ## Introduction -This vignette will walk you through the analyses presented in [Winton et al. 2018](https://besjournals.onlinelibrary.wiley.com/doi/full/10.1111/2041-210X.13080), who describes the use of spatial point process models to estimate individual centers of activity (COA) from passive acoustic telemetry data. This vignette walks through how to prepare the data, using the models, and interpretating the results. We will be using the simplest case, which assumes that detection probabilities/receiver detection ranges remain constant over time, to a more complex application of a test-tag integrated model, that incorporates detection data from one or more stationary test transmitters to estimate time-varying detection ranges. The models are fitted in a Bayesian framework using the Stan software ([Carpenter et al. 2017](https://www.jstatsoft.org/article/view/v076i01/0)); code was modified from models found in [Royle et al. 2013](https://www.sciencedirect.com/book/monograph/9780124059399/spatial-capture-recapture) for fitting spatial point process models to data from camera traps. We prefer the Bayesian approach for COA estimation due to the treatment of uncertainty, but realize the longer computational time required may be prohibitive for some applications. We'd also like to note that the models described can support varying degrees of complexity - not all applications will require (or have the data to support) the most complex version of the model. The simpler the model, the shorter the run-time. +This vignette will walk you through the analyses presented in [Winton et al. 2018](https://besjournals.onlinelibrary.wiley.com/doi/full/10.1111/2041-210X.13080), who describes the use of spatial point process models to estimate individual centers of activity (COA) from passive acoustic telemetry data. We will walk through how to prepare the data, using the models, and interpretating the results. We will be using the simplest case, which assumes that detection probabilities/receiver detection ranges remain constant over time, to a more complex application of a test-tag integrated model, that incorporates detection data from one or more stationary test transmitters to estimate time-varying detection ranges. The models are fitted in a Bayesian framework using the Stan software ([Carpenter et al. 2017](https://www.jstatsoft.org/article/view/v076i01/0)); code was modified from models found in [Royle et al. 2013](https://www.sciencedirect.com/book/monograph/9780124059399/spatial-capture-recapture) for fitting spatial point process models to data from camera traps. We prefer the Bayesian approach for COA estimation due to the treatment of uncertainty, but realize the longer computational time required may be prohibitive for some applications. We'd also like to note that the models described can support varying degrees of complexity - not all applications will require (or have the data to support) the most complex version of the model. The simpler the model, the shorter the run-time. We have tried to make the instructions outlined in this vignette user-friendly since we are a group of applied biologists with varying degrees of statistical experience. If some of the statistical notation outlined here or in [Winton et al. 2018](https://besjournals.onlinelibrary.wiley.com/doi/full/10.1111/2041-210X.13080) remains unclear, feel free to contact us with questions for clarification. This is a new package, so if you find bugs, places where code efficiency could be improved, or instances where the documentation could be made more user-friendly, please let us know through the [issues]( https://github.com/trackyverse/TelemetrySpace) on the GitHub repository! @@ -33,6 +33,7 @@ We have several example datasets along with several functions to assist and stre First, we load all the packages needed to carry out the analysis. + ``` r { library(bayesplot) @@ -141,7 +142,7 @@ We can notice that this detection data consists of 5 columns with 577 rows. To b For our detection data, we have a few things we need to do, the first is we need to build a time bin that we will create COAs. For this data we are going to use 1 hour but this time bin can range from 30 mins - 1 day or more and depends on the questions you are asking and the species you are working with. -Let's build our time bins +Let's build our time bins. ``` r @@ -220,7 +221,7 @@ There are a few things to know about running a Bayesian analysis, we suggest rea ### Priors -Bayesian analyses rely on supplying uninformed or informed prior distributions for each parameter (coefficient; predictor) in the model. For the puproses of the deteciton probablity model the following priors are assumed and are not adjustable. To learn more about the structure of the priors and likelihoods see [Standard detection probability](https://telemetryspace.trackyverse.org/articles/coa_standard_gaussian_model.html). +Bayesian analyses rely on supplying uninformed or informed prior distributions for each parameter (coefficient; predictor) in the model. For the puproses of the deteciton probablity model the following priors are assumed and are not adjustable. To learn more about the structure of the priors and likelihoods see [standard detection probability model](https://telemetryspace.trackyverse.org/articles/coa_standard_gaussian_model.html). $$ \begin{aligned} @@ -235,7 +236,7 @@ The Cauchy priors on $\alpha_0$ and $\alpha_1$ are weakly informative and regularize the intercept and decay-rate parameters toward zero while allowing heavy tails. Because both parameters carry explicit `lower`/`upper` bounds in the `parameters` block, the priors are implicitly truncated to -those bounds — Stan automatically renormalizes the density over the +those bounds. Stan automatically renormalizes the density over the constrained support, so no separate normalizing constant needs to be added by hand. The activity-center coordinates $s_{x,i,t}$ and $s_{y,i,t}$ have no explicit sampling statement, so they receive Stan's implicit flat (uniform) @@ -278,8 +279,8 @@ m <- COA_Standard( #> #> SAMPLING FOR MODEL 'COA_Standard_gaussian' NOW (CHAIN 1). #> Chain 1: -#> Chain 1: Gradient evaluation took 0.000967 seconds -#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 9.67 seconds. +#> Chain 1: Gradient evaluation took 0.000997 seconds +#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 9.97 seconds. #> Chain 1: Adjust your expectations accordingly! #> Chain 1: #> Chain 1: @@ -296,15 +297,15 @@ m <- COA_Standard( #> Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling) #> Chain 1: Iteration: 2000 / 2000 [100%] (Sampling) #> Chain 1: -#> Chain 1: Elapsed Time: 12.964 seconds (Warm-up) -#> Chain 1: 12.452 seconds (Sampling) -#> Chain 1: 25.416 seconds (Total) +#> Chain 1: Elapsed Time: 12.6 seconds (Warm-up) +#> Chain 1: 12.083 seconds (Sampling) +#> Chain 1: 24.683 seconds (Total) #> Chain 1: #> #> SAMPLING FOR MODEL 'COA_Standard_gaussian' NOW (CHAIN 2). #> Chain 2: -#> Chain 2: Gradient evaluation took 0.00096 seconds -#> Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 9.6 seconds. +#> Chain 2: Gradient evaluation took 0.000985 seconds +#> Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 9.85 seconds. #> Chain 2: Adjust your expectations accordingly! #> Chain 2: #> Chain 2: @@ -321,15 +322,16 @@ m <- COA_Standard( #> Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling) #> Chain 2: Iteration: 2000 / 2000 [100%] (Sampling) #> Chain 2: -#> Chain 2: Elapsed Time: 11.924 seconds (Warm-up) -#> Chain 2: 10.856 seconds (Sampling) -#> Chain 2: 22.78 seconds (Total) +#> Chain 2: Elapsed Time: 16.519 seconds (Warm-up) +#> Chain 2: 10.421 seconds (Sampling) +#> Chain 2: 26.94 seconds (Total) #> Chain 2: #> warmup sample -#> chain:1 12.964 12.452 -#> chain:2 11.924 10.856 +#> chain:1 12.600 12.083 +#> chain:2 16.519 10.421 ``` +### Convergance and model performance We can inspect the object created with the first object containing the Stan model. To view and work with the model itself call `m$model`. However, as you will see there @@ -350,7 +352,7 @@ summary(m) #> param_draws 10 tbl_df list #> generated_quantities 1 -none- list ``` -Let's look at our trace plots for the model parameters and posterior distributions to ensure the model converged properly. Remember the trace plots should look grassy or caterpillar like. A trace plot is the posterior draw for a given iteration plotted with the iteration number on the x axis and the posterior value for a given parameter or latent variable on the y. We evaluate this for both chains and want to see that both chains are converging on a smiler posterior draw for a given parameter or latent variable. We will first look at parameters of the model. +Let's look at our trace plots for the model parameters and posterior distributions to ensure the model converged properly. Remember the trace plots should look grassy or caterpillar like. A trace plot is the posterior draw for a given iteration plotted with the iteration number on the x-axis and the posterior value for a given parameter or latent variable on the y-axis. We evaluate this for both chains and want to see that both chains are converging on a similar posterior draw for a given parameter or latent variable. We will first look at parameters of the model. ``` r @@ -371,13 +373,12 @@ stan_dens( ![](figure/trace plots and post dist-2.png) -We can see our trace plots look good and that the posterior distributions of our paramaters look good. +We can see our trace plots look good and that the posterior distributions of our parameters look good. Next let’s look at our latent variables which are `sx` and `sy` or the estimated locations of the fish based on the detection probability. For large/lengthy time periods it will become cumbersome to evaluate the posteriors this way and we recommend inspecting $\hat R$ and ESS. - ``` r stan_trace(m$model, pars = c("sx[1,1]", "sx[1,2]", "sy[1,1]", "sy[1,2]")) ``` @@ -395,8 +396,10 @@ stan_dens( ``` ![](figure/trace plots and post dist latent-2.png) + Again, we can see that everything looks good and our model has converged. +### Posteriors Now moving on to the other elements in the outputted object from `COA_*()`. @@ -405,9 +408,9 @@ The second element is a table of parameter estimates, `p0` which is the estimate ``` r m$summary -#> mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat -#> p0 0.5018673 0.000965220 0.01893527 0.4670082 0.4887808 0.5021665 0.5152408 0.5360255 384.8492 1.003026 -#> sigma 0.9926248 0.001377877 0.02835745 0.9404584 0.9723176 0.9916287 1.0118301 1.0472763 423.5588 1.004424 +#> mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat +#> p0 0.5027243 0.0009117563 0.01791754 0.4662184 0.4904595 0.5032605 0.5150034 0.5342045 386.1882 0.9979821 +#> sigma 0.9923695 0.0014063514 0.02888206 0.9370784 0.9739250 0.9918397 1.0117980 1.0490890 421.7633 0.9981151 ``` We can see that the mean detection probability at a distance of 0 m is 0.5 or 50% and that the model converged well for these parameters. @@ -417,7 +420,7 @@ The thrid element returned, is the time required to run the model in minutes. No ``` r m$time -#> [1] 0.8032667 +#> [1] 0.8603833 ``` @@ -427,18 +430,18 @@ The fourth element returned, is a summary of the posterior draws. Returned is th ``` r m$summary_draws #> # A tibble: 21 × 4 -#> variable median q2.5 q97.5 -#> -#> 1 alpha0 0.00867 -0.132 0.144 -#> 2 alpha1 0.508 0.456 0.565 -#> 3 sx[1,1] -2.96 -3.28 -2.67 -#> 4 sx[1,2] -2.96 -3.30 -2.61 -#> 5 sx[1,3] -2.12 -2.30 -1.91 -#> 6 sx[1,4] -2.26 -2.48 -2.04 -#> 7 sx[1,5] -2.77 -3.09 -2.47 -#> 8 sx[1,6] -2.29 -2.52 -2.05 -#> 9 sx[1,7] -2.28 -2.48 -2.07 -#> 10 sx[1,8] -2.98 -3.25 -2.73 +#> variable median q2.5 q97.5 +#> +#> 1 alpha0 0.0130 -0.135 0.137 +#> 2 alpha1 0.508 0.454 0.569 +#> 3 sx[1,1] -2.97 -3.29 -2.68 +#> 4 sx[1,2] -2.94 -3.33 -2.66 +#> 5 sx[1,3] -2.12 -2.32 -1.90 +#> 6 sx[1,4] -2.26 -2.46 -2.05 +#> 7 sx[1,5] -2.77 -3.09 -2.45 +#> 8 sx[1,6] -2.30 -2.51 -2.08 +#> 9 sx[1,7] -2.29 -2.48 -2.07 +#> 10 sx[1,8] -2.98 -3.27 -2.75 #> # ℹ 11 more rows ``` @@ -450,15 +453,16 @@ m$coas #> # A tibble: 8 × 8 #> ind time x y x_lower x_upper y_lower y_upper #> -#> 1 1 1 -2.96 -0.518 -3.28 -2.67 -0.882 -0.179 -#> 2 1 2 -2.96 -0.331 -3.30 -2.61 -0.699 0.0925 -#> 3 1 3 -2.12 0.0652 -2.30 -1.91 -0.182 0.249 -#> 4 1 4 -2.26 0.144 -2.48 -2.04 -0.0993 0.396 -#> 5 1 5 -2.77 -0.273 -3.09 -2.47 -0.578 0.0535 -#> 6 1 6 -2.29 -0.715 -2.52 -2.05 -0.987 -0.461 -#> 7 1 7 -2.28 -0.0781 -2.48 -2.07 -0.279 0.125 -#> 8 1 8 -2.98 -0.473 -3.25 -2.73 -0.766 -0.184 +#> 1 1 1 -2.97 -0.510 -3.29 -2.68 -0.848 -0.181 +#> 2 1 2 -2.94 -0.305 -3.33 -2.66 -0.683 0.0437 +#> 3 1 3 -2.12 0.0433 -2.32 -1.90 -0.170 0.268 +#> 4 1 4 -2.26 0.154 -2.46 -2.05 -0.0796 0.359 +#> 5 1 5 -2.77 -0.258 -3.09 -2.45 -0.593 0.0787 +#> 6 1 6 -2.30 -0.717 -2.51 -2.08 -1.02 -0.456 +#> 7 1 7 -2.29 -0.0833 -2.48 -2.07 -0.295 0.139 +#> 8 1 8 -2.98 -0.472 -3.27 -2.75 -0.728 -0.194 ``` + The sixth element returned, is a `data.frame` containing the posterior draws from each non-warm-up iteration from all chains. This contains the posterior distribution for each parameter and latent variable for each individual in each time step. It is unlikely that you will use this object instead we have created further objects that organize this data for specific applications that you are more likely to use. @@ -466,16 +470,16 @@ The sixth element returned, is a `data.frame` containing the posterior draws fro m$all_estimates #> # A draws_df: 200 iterations, 2 chains, and 21 variables #> alpha0 alpha1 sx[1,1] sx[1,2] sx[1,3] sx[1,4] sx[1,5] sx[1,6] -#> 1 0.0491 0.52 -2.9 -2.9 -2.3 -2.2 -2.7 -2.3 -#> 2 0.1010 0.51 -3.1 -2.8 -2.1 -2.1 -2.8 -2.2 -#> 3 0.0498 0.52 -3.2 -2.6 -2.2 -2.1 -2.8 -2.1 -#> 4 0.0020 0.51 -2.8 -3.0 -2.2 -2.1 -2.9 -2.2 -#> 5 -0.0043 0.49 -2.9 -3.1 -2.1 -2.2 -3.0 -2.2 -#> 6 -0.0358 0.47 -2.9 -3.0 -2.1 -2.3 -2.6 -2.3 -#> 7 0.0061 0.53 -2.8 -3.1 -2.1 -2.2 -2.9 -2.2 -#> 8 0.0050 0.50 -2.7 -3.0 -2.2 -2.3 -2.7 -2.3 -#> 9 -0.0139 0.49 -2.9 -3.0 -2.2 -2.1 -2.9 -2.4 -#> 10 0.0076 0.49 -2.9 -3.0 -2.2 -2.0 -2.9 -2.3 +#> 1 -0.0551 0.49 -3.0 -2.9 -1.9 -2.3 -2.8 -2.4 +#> 2 0.0096 0.48 -2.9 -2.8 -2.1 -2.2 -3.1 -2.1 +#> 3 0.0537 0.49 -2.9 -3.1 -2.1 -2.4 -2.9 -2.3 +#> 4 0.1017 0.59 -2.9 -2.7 -1.9 -2.3 -2.6 -2.2 +#> 5 0.0532 0.47 -3.2 -2.9 -2.4 -2.5 -2.9 -2.5 +#> 6 0.0821 0.48 -2.9 -2.9 -2.2 -2.4 -2.9 -2.3 +#> 7 -0.0076 0.46 -3.1 -3.1 -2.1 -2.4 -2.8 -2.2 +#> 8 0.1247 0.54 -2.8 -3.0 -2.2 -2.2 -2.7 -2.2 +#> 9 -0.0047 0.50 -2.8 -3.0 -2.2 -2.2 -2.7 -2.1 +#> 10 -0.0147 0.48 -3.0 -2.8 -1.9 -2.3 -2.6 -2.2 #> # ... with 390 more draws, and 13 more variables #> # ... hidden reserved variables {'.chain', '.iteration', '.draw'} ``` @@ -490,18 +494,18 @@ Later on we will plot this object. loc_draws <- m$loc_draws loc_draws #> # A tibble: 3,200 × 8 -#> .chain .iteration .draw lp__ fish time x y -#> -#> 1 1 1 1 -1279. 1 1 -2.85 -0.723 -#> 2 1 1 1 -1279. 1 2 -2.87 -0.175 -#> 3 1 1 1 -1279. 1 3 -2.27 0.129 -#> 4 1 1 1 -1279. 1 4 -2.20 0.262 -#> 5 1 1 1 -1279. 1 5 -2.75 -0.190 -#> 6 1 1 1 -1279. 1 6 -2.26 -0.518 -#> 7 1 1 1 -1279. 1 7 -2.32 -0.184 -#> 8 1 1 1 -1279. 1 8 -2.95 -0.400 -#> 9 1 2 2 -1281. 1 1 -3.07 -0.440 -#> 10 1 2 2 -1281. 1 2 -2.80 -0.338 +#> .chain .iteration .draw lp__ fish time x y +#> +#> 1 1 1 1 -1280. 1 1 -3.02 -0.775 +#> 2 1 1 1 -1280. 1 2 -2.95 -0.397 +#> 3 1 1 1 -1280. 1 3 -1.94 0.0106 +#> 4 1 1 1 -1280. 1 4 -2.33 0.181 +#> 5 1 1 1 -1280. 1 5 -2.80 -0.410 +#> 6 1 1 1 -1280. 1 6 -2.38 -0.634 +#> 7 1 1 1 -1280. 1 7 -2.21 -0.0885 +#> 8 1 1 1 -1280. 1 8 -3.07 -0.541 +#> 9 1 2 2 -1284. 1 1 -2.90 -0.821 +#> 10 1 2 2 -1284. 1 2 -2.83 -0.563 #> # ℹ 3,190 more rows ``` @@ -513,18 +517,18 @@ Later on we will plot this object. param_draws <- m$param_draws param_draws #> # A tibble: 6,400 × 10 -#> .chain .iteration .draw lp__ fish time alpha0 alpha1 p0 sigma -#> -#> 1 1 1 1 -1279. 1 1 0.0491 0.519 0.512 0.982 -#> 2 1 1 1 -1279. 1 2 0.0491 0.519 0.512 0.982 -#> 3 1 1 1 -1279. 1 3 0.0491 0.519 0.512 0.982 -#> 4 1 1 1 -1279. 1 4 0.0491 0.519 0.512 0.982 -#> 5 1 1 1 -1279. 1 5 0.0491 0.519 0.512 0.982 -#> 6 1 1 1 -1279. 1 6 0.0491 0.519 0.512 0.982 -#> 7 1 1 1 -1279. 1 7 0.0491 0.519 0.512 0.982 -#> 8 1 1 1 -1279. 1 8 0.0491 0.519 0.512 0.982 -#> 9 1 1 1 -1279. 1 1 0.0491 0.519 0.512 0.982 -#> 10 1 1 1 -1279. 1 2 0.0491 0.519 0.512 0.982 +#> .chain .iteration .draw lp__ fish time alpha0 alpha1 p0 sigma +#> +#> 1 1 1 1 -1280. 1 1 -0.0551 0.485 0.486 1.02 +#> 2 1 1 1 -1280. 1 2 -0.0551 0.485 0.486 1.02 +#> 3 1 1 1 -1280. 1 3 -0.0551 0.485 0.486 1.02 +#> 4 1 1 1 -1280. 1 4 -0.0551 0.485 0.486 1.02 +#> 5 1 1 1 -1280. 1 5 -0.0551 0.485 0.486 1.02 +#> 6 1 1 1 -1280. 1 6 -0.0551 0.485 0.486 1.02 +#> 7 1 1 1 -1280. 1 7 -0.0551 0.485 0.486 1.02 +#> 8 1 1 1 -1280. 1 8 -0.0551 0.485 0.486 1.02 +#> 9 1 1 1 -1280. 1 1 -0.0551 0.485 0.486 1.02 +#> 10 1 1 1 -1280. 1 2 -0.0551 0.485 0.486 1.02 #> # ℹ 6,390 more rows ``` @@ -535,7 +539,7 @@ The ninth element returned is a `list` that contains generated quantities for `y yrep <- m$generated_quantities str(yrep) #> List of 1 -#> $ yrep: int [1:10, 1:640] 1 2 2 0 4 0 0 2 0 1 ... +#> $ yrep: int [1:10, 1:640] 1 0 1 1 1 0 2 0 1 0 ... #> ..- attr(*, "dimnames")=List of 2 #> .. ..$ : chr [1:10] "yrep_1" "yrep_2" "yrep_3" "yrep_4" ... #> .. ..$ : chr [1:640] "tag_1_rec_1_time_1" "tag_1_rec_2_time_1" "tag_1_rec_3_time_1" "tag_1_rec_4_time_1" ... @@ -546,7 +550,8 @@ str(yrep) First, we are going to plot the posterior draws for locations. Within the package we have a `sf` object that is the shape of Parry Sound. We can use this to plot the posterior draws of `sx` and `sy` to understand the movement of that of lake trout for 8 hours. -We need to take that `sf` object and transform it into an aeqd projection. +We need to take that `sf` object and transform it into an aeqd projection. + ``` r ps_aeqd <- ps |> @@ -631,7 +636,7 @@ p_param ![](figure/plot param-1.png) -We can see that they are all quite tightly distributed with `p0` indicating that detection probability at a distance of `0` is between 48 - 54 %, while the `sigma` is around 1 km. +We can see that they are all quite tightly distributed with `p0` indicating that detection probability at a distance of `0` is between 48 - 54 %, while the `sigma` is around 1 km. Considering this is the standard model, `p0` nor other parameters vary. Considering how acoustic telemetry functions, we know that this is an unlikely represenation of these parameters which we can use a time-varying model and a tag-integrated, time-varying model to account for variation in the `p0` and other parameters. Lastly, we can plot the predictive posterior check using `{tidybayes}`. First we need to make detection counts as a vector from our `build_count()` object. @@ -640,8 +645,10 @@ detection counts as a vector from our `build_count()` object. ``` r y_obs <- as.vector(ps_count_example[!is.na(ps_count_example)]) ``` + Next we can plot the densities using `ppc_dens_overlay()` from `{bayesplot}`. + ``` r ppc <- ppc_dens_overlay(y = y_obs, yrep = yrep$yrep) ppc diff --git a/vignettes/figure/plot param-1.png b/vignettes/figure/plot param-1.png index a8664dc..ced4e7f 100644 Binary files a/vignettes/figure/plot param-1.png and b/vignettes/figure/plot param-1.png differ diff --git a/vignettes/figure/post densities combined-1.png b/vignettes/figure/post densities combined-1.png index be94488..ff7c831 100644 Binary files a/vignettes/figure/post densities combined-1.png and b/vignettes/figure/post densities combined-1.png differ diff --git a/vignettes/figure/post densities timestep-1.png b/vignettes/figure/post densities timestep-1.png index df14b41..ff6cdb1 100644 Binary files a/vignettes/figure/post densities timestep-1.png and b/vignettes/figure/post densities timestep-1.png differ diff --git a/vignettes/figure/ppc dens-1.png b/vignettes/figure/ppc dens-1.png index deb17b8..a941272 100644 Binary files a/vignettes/figure/ppc dens-1.png and b/vignettes/figure/ppc dens-1.png differ diff --git a/vignettes/figure/trace plots and post dist latent-1.png b/vignettes/figure/trace plots and post dist latent-1.png index 796c362..8ffebc5 100644 Binary files a/vignettes/figure/trace plots and post dist latent-1.png and b/vignettes/figure/trace plots and post dist latent-1.png differ diff --git a/vignettes/figure/trace plots and post dist latent-2.png b/vignettes/figure/trace plots and post dist latent-2.png index 2e71b54..64cc93d 100644 Binary files a/vignettes/figure/trace plots and post dist latent-2.png and b/vignettes/figure/trace plots and post dist latent-2.png differ diff --git a/vignettes/figure/trace plots and post dist-1.png b/vignettes/figure/trace plots and post dist-1.png index d43b66c..a323fad 100644 Binary files a/vignettes/figure/trace plots and post dist-1.png and b/vignettes/figure/trace plots and post dist-1.png differ diff --git a/vignettes/figure/trace plots and post dist-2.png b/vignettes/figure/trace plots and post dist-2.png index 40e44a6..63a809d 100644 Binary files a/vignettes/figure/trace plots and post dist-2.png and b/vignettes/figure/trace plots and post dist-2.png differ diff --git a/vignettes/precompile.R b/vignettes/precompile.R index 054aca3..7b14a98 100644 --- a/vignettes/precompile.R +++ b/vignettes/precompile.R @@ -31,8 +31,8 @@ knit( ) knit( - "_ps_estimate_coa_vignette.Rmd", - "ps_estimate_coa_vignette.Rmd" + "_estimate_ps_coa_vignette.Rmd", + "estimate_ps_coa_vignette.Rmd" ) knit( "_coa_standard_gaussian_model.Rmd",