diff --git a/DESCRIPTION b/DESCRIPTION index 883ffc5..b2077f4 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -27,6 +27,7 @@ Imports: stats, tidyr Suggests: + bayesplot, ggplot2, ggpubr, hexbin, diff --git a/vignettes/_ps_estimate_coa_vignette.Rmd b/vignettes/_ps_estimate_coa_vignette.Rmd index c9be830..e98877b 100644 --- a/vignettes/_ps_estimate_coa_vignette.Rmd +++ b/vignettes/_ps_estimate_coa_vignette.Rmd @@ -12,40 +12,50 @@ vignette: > ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, - comment = "#>" + comment = "#>", + fig.width = 10, + fig.height = 7 ) ``` ## 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 progresses walks through how to prepare the data, using the model, 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 its 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. 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. -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 the paper 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! +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! ## Models + For simplicity purposes of this vignette, model specifications can be found in the following documentation: + +- [Standard detection probability](https://telemetryspace.trackyverse.org/articles/coa_standard_gaussian_model.html) (i.e., `COA_standard()`) +- [Time-varying detection probability]() (i.e., `COA_TimeVarying ()`) +- [Tag-integrated time-varying detection probability]() (i.e., `COA_TagInt()`) + ## Data preparation -To run spatial point process/detection probablity models in Stan we need to do a couple of things to our detection and receiver location data. This includes first providing the model with all the necessary information as well as in the proper format. Stan likes to handle data in vectors and arrays, which we don't often use in R, that are then provided to Stan in a `list`. To learn more about Stan, please click this [link](https://mc-stan.org/). Stan is written in C++, making it quite fast, with Stan programs following a very systematic order, see [link to understand a Stan program](). Stan uses [Markov chain Monte Carlo (MCMC)](https://en.wikipedia.org/wiki/Markov_chain_Monte_Carlo) and [No U-Turn Sampler (NUTS)](https://mc-stan.org/docs/2_18/stan-users-guide/sampling-difficulties-with-problematic-priors.html) to improve efficiency. We will now walk through how to setup the data to be able to run the models. - -We have serveral 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 efficent to calcuate distances among receivers, however, the distance between receivers needs to be on the same x and y plain and ranges need to be between 0.2 - 15 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 (*Salvlinus namaycush*) that was implanted with an acoustic transmitter in Parry Sound which is a large embayment of Geogian Bay, Lake Huron. +To run spatial point process/detection probability models in Stan we need to do a couple of things to our detection and receiver location data. This includes first providing the model with all the necessary information as well as in the proper format. Stan likes to handle data using vectors and arrays, which we often use in R but they are configured differently (e.g., `data.frame()` is a bunch vectors of different classes, while a `matrix` is a bunch of vectors of the same class). These vectors and arrays need to be provided to Stan in a `list`. To learn more about Stan, please click this [link](https://mc-stan.org/). Stan is written in C++, making it quite fast, with Stan programs following a very systematic order, see [link to understand a Stan program]( https://mc-stan.org/docs/reference-manual/blocks.html). Stan uses a [Hamilton Monte Carlo](https://mc-stan.org/docs/2_18/reference-manual/hamiltonian-monte-carlo.html) which is a [Markov chain Monte Carlo (MCMC)](https://en.wikipedia.org/wiki/Markov_chain_Monte_Carlo) that obtains a sequence of random samples whose distribution converges to a target probability distribution that is difficult to sample directly. Stan also employs a [No U-Turn Sampler (NUTS)](https://mc-stan.org/docs/2_18/stan-users-guide/sampling-difficulties-with-problematic-priors.html) to improve efficiency. We will now walk-through how-to setup the data to be able to run the models. -First we load all the packages needed to carry out the analysis. +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) library(ggplot2) library(rstan) library(sf) library(TelemetrySpace) + library(tidyr) } ``` ### Receiver Locations -Next we will look at the location of the receivers in Parry Sound. To build a aqed projection we need the locations to -be in meters which we can transform these locations to a meter focused coordinate reference systme (crs), in this case UTMs. You will notice that `ps_rec_loc` is a `data.frame`. We need to make it into an `sf` object to then transforme it into UTMs and then build the aeqd projection. We +Next, we will look at the location of the receivers in Parry Sound. To build an aqed projection we need the locations to +be in meters which we can transform these locations to a meter focused coordinate reference system (crs), in this case a projected coordinate system (e.g., UTMs). You will notice that `ps_rec_loc` is a `data.frame`. We need to make it into an `sf` object which is in a geographical coordinate reference system (e.g., WGS 84; degrees longitude and latitude). We need to transform it into a projected coordinate system, in this case we used UTMs. We then can build the aeqd projection. We can do this using `st_as_sf()` from [{sf}](https://r-spatial.github.io/sf/). + ```{r recs inspect} # first look at the example head(ps_rec_loc) @@ -55,31 +65,33 @@ str(ps_rec_loc) ps_rec_loc_sf <- ps_rec_loc |> st_as_sf(coords = c("deploy_long", "deploy_lat"), crs = 4326) -# next transform it into the correct utm +# next transform it into projected coordinate system (e.g., utm ESPG:32617) ps_rec_loc_utm <- ps_rec_loc_sf |> st_transform(32617) ``` -Next we will build an aeqd projection for this receiver array that we can then transform the locations of the receivers into. +Next, we will build an aeqd projection for this receiver array that we can then transform the locations of the receivers into. ```{r aeqd} aeqd_crs <- build_aeqd(ps_rec_loc_utm) ``` -We then can transform the receiver locations into our aeqd project +We then can transform the receiver locations into our aeqd project and make sure they are in the correct order. ```{r transform aqed} ps_rec_loc_aeqd <- ps_rec_loc_sf |> st_transform(aeqd_crs) |> (\(.) .[order(.$station_no), ])() ``` -Now that we the receiver locations transformed we need to first index them appropiately, this is because Stan will not be able to handel +Now that we have the receiver locations transformed, we need to first index them appropriately, this is because Stan will not be able to handle the `station_no` but instead can handle a numerical index value to identify the receivers. + ```{r index rec} ps_rec_loc_aeqd$rec <- 1:nrow(ps_rec_loc_aeqd) ``` -Next we will transform this into a `data.frame` the contains two `vectors` that we can supply to the Stan model. We also need to build +Next, we will transform this into a `data.frame` the contains two columns (i.e.,`vectors`) that we can supply to the Stan model. We also need to build the boundary box of the receiver array. The argument `buffer` in `build_bbox()` assumes a 1 km buffer but this can be changed depending on how you would like the boundary box to be created. + ```{r build rec} rec_loc_vec <- build_rec_coords(ps_rec_loc_aeqd) @@ -88,7 +100,7 @@ rec_limits <- build_bbox(rec_loc_vec) ### Detection Data -Now that we have the receiver locations in a format that Stan can handle, we are going to prepare the detection data. First lets look at the detection data. +Now that we have the receiver locations in a format that Stan can handle, we are going to prepare the detection data. First let’s look at the detection data. ```{r inspect det data} head(ps_det_example) str(ps_det_example) @@ -107,7 +119,7 @@ head(ps_det_example) You will notice both a `POSIXct` column that is called `time_bin` and a numerical column called `time`. This `time` column is a numerical index of the time bins. -We have a few more things we need to do to prepare the data, first we need to add in the numerical index of the receiver values that we created in the first section. We can do this by using `merge()` from base or we want we could use `left_join()` from `{dplyr}` or `merge.data.table()` from `{data.table}`. +We have a few more things we need to do to prepare the data, first we need to add in the numerical index of the receiver values that we created in the first section. We can do this by using `merge()` from base or we could use `left_join()` from `{dplyr}` or `merge.data.table()` from `{data.table}`. ```{r merge rec} ps_det_example <- merge( ps_det_example, @@ -116,9 +128,9 @@ ps_det_example <- merge( ) head(ps_det_example) ``` -We are starting to get somewhere, you can see that we now have the time and recevier index values lined up. The last two things we need to do is first determine the numer of individuals, in this case it will be 1 and the expetected number of detections for invidiaul for each receiver for each time bin, and the number of transmissions this is where `min_delay` and `max_delay` come into play. We do suggest running the model for each invidual with modeling large time periods for example 1 year's worth of data being broken down into 7 day chunks. +We are starting to get somewhere; you can see that we now have the time and receiver index values lined up. The last two things we need to do is first determine the number of individuals, in this case it will be 1 and the expected number of detections for an individual for each receiver for each time bin, and the number of transmissions this is where `min_delay` and `max_delay` come into play. -Lets build the count data and the number of time steps in the data. +Let’s build the count data and the number of time steps in the data. ```{r build count} ps_count_example <- build_counts( @@ -131,20 +143,21 @@ ps_count_example <- build_counts( time_steps <- build_tstep(ps_count_example) ``` -Now we can create the number of inviduals and the number of transmisisons expected within the time bin. +Now we can create the number of individuals and the number of transmissions expected within the time bin. ```{r build nind and trans} nind <- length(unique(ps_det_example$tag_serial_no)) ntrans <- build_ntrans(ps_det_example) ``` -We can finally move on to running the model, however, we need to understand a bit more about Bayesian frameworks before we proceed. +We can finally move on to running the model; however, we need to understand a bit more about Bayesian frameworks before we proceed. ## Bayesian Analysis -We can now estimate the center of activity for lake trout in Parry Sound, Lake Huron. +We can now estimate the center of activity for lake trout in Parry Sound +which is a large embayment of Georgian Bay, Lake Huron. -There are a few things to know about running a Bayesian analysis, I suggest reading these resources: +There are a few things to know about running a Bayesian analysis, we suggest reading these resources: 1. [Basics of Bayesian Statistics - Book](https://statswithr.github.io/book/) 2. [A Student’s Guide to Bayesian Statistics by Ben Lambert](https://sites.math.rutgers.edu/~zeilberg/EM20/Lambert.pdf) @@ -155,14 +168,34 @@ There are a few things to know about running a Bayesian analysis, I suggest read ### 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. +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). + +$$ +\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. ### Model convergence -It is important to always run the model with at least 2 chains. If the model does not converge you can try to increase the following: +It is important to always run the model with at least 2 chains. If the model does not converge you can try to increase or decreasing the following: -1. The amount of samples that are discarded; this can be controlled by the argument `warmup`. +1. The number of samples that are discarded; this can be controlled by the argument `warmup`. 2. The number of iterative samples retained; this can be controlled by the argument `iter`. @@ -171,10 +204,10 @@ It is important to always run the model with at least 2 chains. If the model doe 4. The `adapt_delta` value using `control = list(adapt_delta = 0.95)`. -When assessing the model we want $\hat R$ to be 1 or within 0.05 of 1, which indicates that the variance among and within chains are equal (see [{rstan} documentation on $\hat R$](https://mc-stan.org/rstan/reference/Rhat.html)), a high value for effective sample size (ESS), trace plots to look "grassy" or "caterpillar like," and posterior distributions to look relatively normal. +When assessing the model, we want $\hat R$ to be 1 or within 0.05 of 1, which indicates that the variance among and within chains are equal (see [{rstan} documentation on $\hat R$](https://mc-stan.org/rstan/reference/Rhat.html)), a high value for effective sample size (ESS), trace plots to look "grassy" or "caterpillar like," and posterior distributions to look relatively normal. ## Model -We can now run the standard point process/detection probablity model. +We can now run the standard point process/detection probability model. ```{r run coa_standard} m <- COA_Standard( nind = nind, @@ -192,12 +225,15 @@ m <- COA_Standard( ``` -We can inspect the model object with the first object containing the Stan model object (accessible via `m$model`, which will allow you to use `rstan` plotting tools and diagnostic plots - see rstan documentation for details). +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 +Are a bunch of other objects in the objected created. These objects have pulled important information from the Stan model. + + ```{r summary of model object} summary(m) ``` - -Let's look at our trace plots for the model parameters and posterior distrubitons to ensure the model converged proeprly. +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. ``` {r trace plots and post dist} stan_trace(m$model, pars = c("alpha0", "alpha1", "p0", "sigma")) @@ -212,10 +248,10 @@ stan_dens( We can see our trace plots look good and that the posterior distributions of our paramaters look good. -Next lets look at our latent variables which are `sx` and `sy` or the estiamted locations of the fish based on -the detection probablity. For large/lengthy time periods it will become cumbersome to evaluate the posteriors +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} stan_trace(m$model, pars = c("sx[1,1]", "sx[1,2]", "sy[1,1]", "sy[1,2]")) @@ -226,59 +262,74 @@ stan_dens( linewidth = 0.1 ) ``` -Again we can see that everything looks good and our model has converted. +Again, we can see that everything looks good and our model has converged. + +Now moving on to the other elements in the outputted object from `COA_*()`. -Now moving on to the other elements in the outputed object from `COA_*()`. -The second element is a table of parameter estimates, `p0` which is the estimated detection probality at distance of 0 m and the standard deviation (`sigma`). Their mean, standard error of the mean, and the associated quantiles from the posterior distriubtions. The table also includes the effective sample size and the Rhat statistic (which should be between 0.95 and 1.05). +The second element is a table of parameter estimates, `p0` which is the estimated detection probability at distance of 0 m and $\sigma$ which is units of distance and characterizes the effective spatial range of detection: it is the standard deviation of the Gaussian or logistic decay in detection log-odds with distance from the activity center. A larger $\sigma$ corresponds to a more slow-decaying (longer-range) detection function (`sigma`). The table contains their mean, standard error of the mean, and the associated quantiles from the posterior distributions. The table also includes the effective sample size and the $\hat R$ statistic (which should be between 0.95 and 1.05). ```{r summary from COA} m$summary ``` -We can see that the mean detection probablity at a distance of 0 m is 0.5 or 50% and that the model converged well for this paramaters. +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. -The thrid element returned, is the time required to run the model in minutes. Note that Stan will automatically detect and use multiple cores. If the computer used to run this has multiple cores, the time returned will be longer than the actual run time. This is because it will sum the time for each core. To return the realized run time, divide `fit$time` by the number of cores. +The thrid element returned, is the time required to run the model in minutes. Note that Stan will automatically detect and use multiple cores. If the computer used to run this has multiple cores, the time returned will be longer than the actual run time. This is because it will sum the time for each core. To return the realized run time, divide `m$time` by the number of cores. ```{r runtime } m$time ``` -The fouth element returned, is a summary of the posterior draws. Returned is the `variable`, `median`, and the 2.5 and 97.5% quartiles (e.g., `q2.5` and `q97.5`). This summary `data.frame` provides you the ability to inspect all elements, however, most of the time a you will be using this information that has been filtered and displayed in a more user frinedly manner in futher objects returned by `COA_*()` functions. +The fourth element returned, is a summary of the posterior draws. Returned is the `variable`, `median`, and the 2.5 and 97.5% quartiles (e.g., `q2.5` and `q97.5`). This summary `data.frame` provides you the ability to inspect all elements, however, most of the time you will be using this data in a filtered and transformed manner created in other objects in the model output. These objects are more user friendly for plotting and understanding the results. ```{r summary_draws} m$summary_draws ``` -The fourth returns `data.frame` of the median values for the latent variables `sx` and `sy` which are the estimated location of the fish. The median values can be viewed as the center of activity while the posterior distibution is the uncertainity around that estimate. The `data.frame` has a few important compoents `ind` is the index number of the invidiaul, `time` is the index number for a given time bin, `x` and `y` are the median values for a given individual and time, while `x_*` and `y_*` are the 2.5 and 97.5% quantiles for those posteriors, this effectivively is viewed as the 95% credible interval (Bayesian version of a confidence interval) for each coordinate. +The fifth element returned, is `data.frame` of the median values for the latent variables `sx` and `sy` which are the estimated location of the fish. The median values can be viewed as the center of activity while the posterior distribution is the overall probability of occurrence around that estimate. The `data.frame` has a few important components `ind` is the index number of the individual, `time` is the index number for a given time bin, `x` and `y` are the median values for a given individual and time, while `x_*` and `y_*` are the 2.5 and 97.5% quantiles for those posteriors, this effectively is viewed as the 95% credible interval, (i.e., Bayesian version of a confidence interval) for each coordinate. ```{r coas} m$coas ``` - -The fith element (`m$all_estimates`) contains the posterior draws from each non-warm-up iteration from all chains. If you're not familiar with Bayesian terminology all you need to know is that this contains disturbtions for each paramter and latent variable for each individual in each time step. It is unlikely that you will use this object instead we have created futher objects below that you are more likely to use. +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} m$all_estimates ``` -The sixth element (`m$loc_draws`), contains an extracted and transformed posterior draws -for the latent variables `sx` and `sy`. This object we will use to plot our estimates of -centers of activity. There are several other columns besides, the draw (i.e, `x`, `y`), the -`fish`/ind index number and the `time` index number. We have `.chain`, `.itteration`, `.draw` which identifies the chain, iteration, and draw that posterior draw came from. We also -have the column `lp__` which is the unnormalized log posterior density of your model, used for diagnosing sampling efficiency, estimating approximations, and comparing models. -We are going to inspect this element and then plot it. +The seventh element returned is a `data.frame` that contains the extracted and transformed posterior draws for the latent variables `sx` and `sy`. This object we will use to plot our estimates of centers of activity. There are several other columns besides, the draws (i.e, `x`, `y`), the `fish`/ind index number and the `time` index number. We have `.chain`, `.iteration`, `.draw` which identifies the chain, iteration, and draw that posterior draw came from. We also have the column `lp__` which is the unnormalized log posterior density of your model, used for diagnosing sampling efficiency, estimating approximations, and comparing models. + +Later on we will plot this object. + ```{r loc draws} loc_draws <- m$loc_draws loc_draws ``` -Next 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 that we plot this onto so we can understand where the fish was during each time period. +The eighth element returned, is a `data.frame` that contains the contains the extracted and transformed posterior draws for the model parameters `alpha0`, `alpah1`, `p0`, and `sigma`. There are several other columns besides, the draws, the `fish`/ind index number and the `time` index number. We have `.chain`, `.iteration`, `.draw` which identifies the chain, iteration, and draw that posterior draw came from. We also have the column `lp__` which is the unnormalized log posterior density of your model, used for diagnosing sampling efficiency, estimating approximations, and comparing models. + +Later on we will plot this object. +```{r param draws} +param_draws <- m$param_draws +param_draws +``` + +The ninth element returned is a `list` that contains generated quantities for `y` from the model. This `y` is different than `sy` and is the estimated counts of detections for each individual at each time step at each receiver. The object `yrep` which is a 2-dimensional array that contains 10 draws (`ndraws` argument controls the number of draws here) and 640 posterior counts of detections the tag number, receiver number and time bin. We will assess to see how well this matches the existing data's distribution as using [posterior probability](https://en.wikipedia.org/wiki/Posterior_probability) which is a type of conditional probability that results from updating the prior probability with information summarized by the likelihood. + +``` {r yrep} +yrep <- m$generated_quantities +yrep +``` + +## Plotting + +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. -First though 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) @@ -303,7 +354,68 @@ p <- ggplot() + p1 <- p + facet_wrap(~time) +``` -p +Let's first look at the posterior densities for each time step + +```{r post densities timestep} p1 ``` + +Next let's look at the densities combined to gather a better understanding of the movement over the 8 timesteps. + +```{r post densities combined} +p +``` + +Next we can plot the posterior distributions for `alpha0`, `alpha1`, `p0`, and `sigma`. We need to first make `fish` and `time` a `character` and then to make plotting easier we need to make the `data.frame` in long format. + +```{r plot params} +param_draws$fish <- as.character(param_draws$fish) +param_draws$time <- as.character(param_draws$time) + +param_draws_long <- param_draws |> + pivot_longer( + cols = -c(".chain", ".iteration", ".draw", "lp__", "fish", "time"), + names_to = "param", + values_to = "est" + ) +``` + +We can now plot those posterior distributions of the model parameters. + +```{r plot param} +p_param <- ggplot( + data = param_draws_long, + aes(x = time, y = est, fill = fish), +) + + geom_violin() + + facet_wrap(~param, scale = "free_y") + + scale_fill_manual(name = "Fish ID", values = "#2768F5") + + theme_bw() + + theme( + panel.grid = element_blank(), + strip.background = element_blank() + ) + + labs(x = "Time Bin", y = "Estimate") + +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. + +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. + +```{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} +ppc <- ppc_dens_overlay(y = y_obs, yrep = yrep$yrep) +ppc +``` + +We can see our predictive posterior check distribution for 10 draws closely lines up with +our observed detection counts indicating that the model fits the data well. + +Congratulations we have successfully run a Bayesian standard point processing/detection probability model for one Lake Trout for 8, 1 hour time bins (8 hrs) in Parry Sound, a large embayment of Georgian Bay, Lake Huron. We can take the results and write up how this individual behaved for a given time period based on the model results! \ No newline at end of file diff --git a/vignettes/figure/plot param-1.png b/vignettes/figure/plot param-1.png new file mode 100644 index 0000000..dc1ccb7 Binary files /dev/null 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 new file mode 100644 index 0000000..641f57e Binary files /dev/null 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 new file mode 100644 index 0000000..0230527 Binary files /dev/null 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 new file mode 100644 index 0000000..51636f6 Binary files /dev/null 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 015f8c4..6528036 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 0a3d947..45eeb61 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 68f579a..93d3f9a 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 45ef393..d51b73b 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/ps_estimate_coa_vignette.Rmd b/vignettes/ps_estimate_coa_vignette.Rmd index d47679e..a6d3137 100644 --- a/vignettes/ps_estimate_coa_vignette.Rmd +++ b/vignettes/ps_estimate_coa_vignette.Rmd @@ -12,37 +12,45 @@ 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 progresses walks through how to prepare the data, using the model, 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 its 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. 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. -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 the paper 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! +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! ## Models + For simplicity purposes of this vignette, model specifications can be found in the following documentation: + +- [Standard detection probability](https://telemetryspace.trackyverse.org/articles/coa_standard_gaussian_model.html) (i.e., `COA_standard()`) +- [Time-varying detection probability]() (i.e., `COA_TimeVarying ()`) +- [Tag-integrated time-varying detection probability]() (i.e., `COA_TagInt()`) + ## Data preparation -To run spatial point process/detection probablity models in Stan we need to do a couple of things to our detection and receiver location data. This includes first providing the model with all the necessary information as well as in the proper format. Stan likes to handle data in vectors and arrays, which we don't often use in R, that are then provided to Stan in a `list`. To learn more about Stan, please click this [link](https://mc-stan.org/). Stan is written in C++, making it quite fast, with Stan programs following a very systematic order, see [link to understand a Stan program](). Stan uses [Markov chain Monte Carlo (MCMC)](https://en.wikipedia.org/wiki/Markov_chain_Monte_Carlo) and [No U-Turn Sampler (NUTS)](https://mc-stan.org/docs/2_18/stan-users-guide/sampling-difficulties-with-problematic-priors.html) to improve efficiency. We will now walk through how to setup the data to be able to run the models. +To run spatial point process/detection probability models in Stan we need to do a couple of things to our detection and receiver location data. This includes first providing the model with all the necessary information as well as in the proper format. Stan likes to handle data using vectors and arrays, which we often use in R but they are configured differently (e.g., `data.frame()` is a bunch vectors of different classes, while a `matrix` is a bunch of vectors of the same class). These vectors and arrays need to be provided to Stan in a `list`. To learn more about Stan, please click this [link](https://mc-stan.org/). Stan is written in C++, making it quite fast, with Stan programs following a very systematic order, see [link to understand a Stan program]( https://mc-stan.org/docs/reference-manual/blocks.html). Stan uses a [Hamilton Monte Carlo](https://mc-stan.org/docs/2_18/reference-manual/hamiltonian-monte-carlo.html) which is a [Markov chain Monte Carlo (MCMC)](https://en.wikipedia.org/wiki/Markov_chain_Monte_Carlo) that obtains a sequence of random samples whose distribution converges to a target probability distribution that is difficult to sample directly. Stan also employs a [No U-Turn Sampler (NUTS)](https://mc-stan.org/docs/2_18/stan-users-guide/sampling-difficulties-with-problematic-priors.html) to improve efficiency. We will now walk-through how-to setup the data to be able to run the models. -We have serveral 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 efficent to calcuate distances among receivers, however, the distance between receivers needs to be on the same x and y plain and ranges need to be between 0.2 - 15 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 (*Salvlinus namaycush*) that was implanted with an acoustic transmitter in Parry Sound which is a large embayment of Geogian Bay, Lake Huron. - -First we load all the packages needed to carry out the analysis. +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 { + library(bayesplot) library(ggplot2) library(rstan) library(sf) library(TelemetrySpace) + library(tidyr) } ``` ### Receiver Locations -Next we will look at the location of the receivers in Parry Sound. To build a aqed projection we need the locations to -be in meters which we can transform these locations to a meter focused coordinate reference systme (crs), in this case UTMs. You will notice that `ps_rec_loc` is a `data.frame`. We need to make it into an `sf` object to then transforme it into UTMs and then build the aeqd projection. We +Next, we will look at the location of the receivers in Parry Sound. To build an aqed projection we need the locations to +be in meters which we can transform these locations to a meter focused coordinate reference system (crs), in this case a projected coordinate system (e.g., UTMs). You will notice that `ps_rec_loc` is a `data.frame`. We need to make it into an `sf` object which is in a geographical coordinate reference system (e.g., WGS 84; degrees longitude and latitude). We need to transform it into a projected coordinate system, in this case we used UTMs. We then can build the aeqd projection. We can do this using `st_as_sf()` from [{sf}](https://r-spatial.github.io/sf/). + ``` r # first look at the example head(ps_rec_loc) @@ -65,35 +73,37 @@ str(ps_rec_loc) ps_rec_loc_sf <- ps_rec_loc |> st_as_sf(coords = c("deploy_long", "deploy_lat"), crs = 4326) -# next transform it into the correct utm +# next transform it into projected coordinate system (e.g., utm ESPG:32617) ps_rec_loc_utm <- ps_rec_loc_sf |> st_transform(32617) ``` -Next we will build an aeqd projection for this receiver array that we can then transform the locations of the receivers into. +Next, we will build an aeqd projection for this receiver array that we can then transform the locations of the receivers into. ``` r aeqd_crs <- build_aeqd(ps_rec_loc_utm) #> ✔ Successfully built "+proj=aeqd +lon_0=-80.124804 +lat_0=45.333008 +x_0=0 +y_0=0 +datum=WGS84 +units=km" ``` -We then can transform the receiver locations into our aeqd project +We then can transform the receiver locations into our aeqd project and make sure they are in the correct order. ``` r ps_rec_loc_aeqd <- ps_rec_loc_sf |> st_transform(aeqd_crs) |> (\(.) .[order(.$station_no), ])() ``` -Now that we the receiver locations transformed we need to first index them appropiately, this is because Stan will not be able to handel +Now that we have the receiver locations transformed, we need to first index them appropriately, this is because Stan will not be able to handle the `station_no` but instead can handle a numerical index value to identify the receivers. + ``` r ps_rec_loc_aeqd$rec <- 1:nrow(ps_rec_loc_aeqd) ``` -Next we will transform this into a `data.frame` the contains two `vectors` that we can supply to the Stan model. We also need to build +Next, we will transform this into a `data.frame` the contains two columns (i.e.,`vectors`) that we can supply to the Stan model. We also need to build the boundary box of the receiver array. The argument `buffer` in `build_bbox()` assumes a 1 km buffer but this can be changed depending on how you would like the boundary box to be created. + ``` r rec_loc_vec <- build_rec_coords(ps_rec_loc_aeqd) @@ -103,26 +113,28 @@ rec_limits <- build_bbox(rec_loc_vec) ### Detection Data -Now that we have the receiver locations in a format that Stan can handle, we are going to prepare the detection data. First lets look at the detection data. +Now that we have the receiver locations in a format that Stan can handle, we are going to prepare the detection data. First let’s look at the detection data. ``` r head(ps_det_example) -#> # A tibble: 6 × 5 -#> detection_timestamp_utc station_no tag_serial_no min_delay max_delay -#> -#> 1 2024-05-03 21:01:53 PSM-007 1594061 190 290 -#> 2 2024-05-03 21:01:54 PSM-003 1594061 190 290 -#> 3 2024-05-03 21:05:25 PSM-004 1594061 190 290 -#> 4 2024-05-03 21:05:25 PSM-007 1594061 190 290 -#> 5 2024-05-03 21:05:25 PSM-008 1594061 190 290 -#> 6 2024-05-03 21:05:26 PSM-003 1594061 190 290 +#> station_no detection_timestamp_utc tag_serial_no min_delay max_delay time_bin time rec.x rec.y +#> 1 PSM-001 2024-05-04 02:09:48 1594061 190 290 2024-05-04 02:00:00 6 1 1 +#> 2 PSM-001 2024-05-04 04:05:03 1594061 190 290 2024-05-04 04:00:00 8 1 1 +#> 3 PSM-001 2024-05-04 03:37:50 1594061 190 290 2024-05-04 03:00:00 7 1 1 +#> 4 PSM-002 2024-05-04 03:11:04 1594061 190 290 2024-05-04 03:00:00 7 2 2 +#> 5 PSM-002 2024-05-04 02:40:09 1594061 190 290 2024-05-04 02:00:00 6 2 2 +#> 6 PSM-002 2024-05-03 22:57:37 1594061 190 290 2024-05-03 22:00:00 2 2 2 str(ps_det_example) -#> tibble [577 × 5] (S3: tbl_df/tbl/data.frame) -#> $ detection_timestamp_utc: POSIXct[1:577], format: "2024-05-03 21:01:53" "2024-05-03 21:01:54" "2024-05-03 21:05:25" "2024-05-03 21:05:25" ... -#> $ station_no : chr [1:577] "PSM-007" "PSM-003" "PSM-004" "PSM-007" ... -#> $ tag_serial_no : chr [1:577] "1594061" "1594061" "1594061" "1594061" ... -#> $ min_delay : num [1:577] 190 190 190 190 190 190 190 190 190 190 ... -#> $ max_delay : num [1:577] 290 290 290 290 290 290 290 290 290 290 ... +#> 'data.frame': 577 obs. of 9 variables: +#> $ station_no : chr "PSM-001" "PSM-001" "PSM-001" "PSM-002" ... +#> $ detection_timestamp_utc: POSIXct, format: "2024-05-04 02:09:48" "2024-05-04 04:05:03" "2024-05-04 03:37:50" "2024-05-04 03:11:04" ... +#> $ tag_serial_no : chr "1594061" "1594061" "1594061" "1594061" ... +#> $ min_delay : num 190 190 190 190 190 190 190 190 190 190 ... +#> $ max_delay : num 290 290 290 290 290 290 290 290 290 290 ... +#> $ time_bin : POSIXct, format: "2024-05-04 02:00:00" "2024-05-04 04:00:00" "2024-05-04 03:00:00" "2024-05-04 03:00:00" ... +#> $ time : int 6 8 7 7 6 2 6 3 3 8 ... +#> $ rec.x : int 1 1 1 2 2 2 2 2 2 2 ... +#> $ rec.y : int 1 1 1 2 2 2 2 2 2 2 ... ``` We can notice that this detection data consists of 5 columns with 577 rows. To better understand what each field is you can run `?ps_det_example` to review the full documentation. @@ -135,20 +147,18 @@ Let's build our time bins ``` r ps_det_example <- build_time_bin(ps_det_example, unit = "1 hour") head(ps_det_example) -#> # A tibble: 6 × 7 -#> detection_timestamp_utc station_no tag_serial_no min_delay max_delay time_bin time -#> -#> 1 2024-05-03 21:01:53 PSM-007 1594061 190 290 2024-05-03 21:00:00 1 -#> 2 2024-05-03 21:01:54 PSM-003 1594061 190 290 2024-05-03 21:00:00 1 -#> 3 2024-05-03 21:05:25 PSM-004 1594061 190 290 2024-05-03 21:00:00 1 -#> 4 2024-05-03 21:05:25 PSM-007 1594061 190 290 2024-05-03 21:00:00 1 -#> 5 2024-05-03 21:05:25 PSM-008 1594061 190 290 2024-05-03 21:00:00 1 -#> 6 2024-05-03 21:05:26 PSM-003 1594061 190 290 2024-05-03 21:00:00 1 +#> station_no detection_timestamp_utc tag_serial_no min_delay max_delay time_bin time rec.x rec.y +#> 1 PSM-007 2024-05-03 21:01:53 1594061 190 290 2024-05-03 21:00:00 1 7 7 +#> 2 PSM-003 2024-05-03 21:01:54 1594061 190 290 2024-05-03 21:00:00 1 3 3 +#> 3 PSM-004 2024-05-03 21:05:25 1594061 190 290 2024-05-03 21:00:00 1 4 4 +#> 4 PSM-007 2024-05-03 21:05:25 1594061 190 290 2024-05-03 21:00:00 1 7 7 +#> 5 PSM-008 2024-05-03 21:05:25 1594061 190 290 2024-05-03 21:00:00 1 8 8 +#> 6 PSM-003 2024-05-03 21:05:26 1594061 190 290 2024-05-03 21:00:00 1 3 3 ``` You will notice both a `POSIXct` column that is called `time_bin` and a numerical column called `time`. This `time` column is a numerical index of the time bins. -We have a few more things we need to do to prepare the data, first we need to add in the numerical index of the receiver values that we created in the first section. We can do this by using `merge()` from base or we want we could use `left_join()` from `{dplyr}` or `merge.data.table()` from `{data.table}`. +We have a few more things we need to do to prepare the data, first we need to add in the numerical index of the receiver values that we created in the first section. We can do this by using `merge()` from base or we could use `left_join()` from `{dplyr}` or `merge.data.table()` from `{data.table}`. ``` r ps_det_example <- merge( @@ -157,17 +167,17 @@ ps_det_example <- merge( by = "station_no" ) head(ps_det_example) -#> station_no detection_timestamp_utc tag_serial_no min_delay max_delay time_bin time rec -#> 1 PSM-001 2024-05-04 02:09:48 1594061 190 290 2024-05-04 02:00:00 6 1 -#> 2 PSM-001 2024-05-04 04:05:03 1594061 190 290 2024-05-04 04:00:00 8 1 -#> 3 PSM-001 2024-05-04 03:37:50 1594061 190 290 2024-05-04 03:00:00 7 1 -#> 4 PSM-002 2024-05-04 03:11:04 1594061 190 290 2024-05-04 03:00:00 7 2 -#> 5 PSM-002 2024-05-04 02:40:09 1594061 190 290 2024-05-04 02:00:00 6 2 -#> 6 PSM-002 2024-05-03 22:57:37 1594061 190 290 2024-05-03 22:00:00 2 2 +#> station_no detection_timestamp_utc tag_serial_no min_delay max_delay time_bin time rec.x rec.y rec +#> 1 PSM-001 2024-05-04 02:09:48 1594061 190 290 2024-05-04 02:00:00 6 1 1 1 +#> 2 PSM-001 2024-05-04 04:05:03 1594061 190 290 2024-05-04 04:00:00 8 1 1 1 +#> 3 PSM-001 2024-05-04 03:37:50 1594061 190 290 2024-05-04 03:00:00 7 1 1 1 +#> 4 PSM-002 2024-05-04 03:11:04 1594061 190 290 2024-05-04 03:00:00 7 2 2 2 +#> 5 PSM-002 2024-05-04 02:40:09 1594061 190 290 2024-05-04 02:00:00 6 2 2 2 +#> 6 PSM-002 2024-05-03 22:57:37 1594061 190 290 2024-05-03 22:00:00 2 2 2 2 ``` -We are starting to get somewhere, you can see that we now have the time and recevier index values lined up. The last two things we need to do is first determine the numer of individuals, in this case it will be 1 and the expetected number of detections for invidiaul for each receiver for each time bin, and the number of transmissions this is where `min_delay` and `max_delay` come into play. We do suggest running the model for each invidual with modeling large time periods for example 1 year's worth of data being broken down into 7 day chunks. +We are starting to get somewhere; you can see that we now have the time and receiver index values lined up. The last two things we need to do is first determine the number of individuals, in this case it will be 1 and the expected number of detections for an individual for each receiver for each time bin, and the number of transmissions this is where `min_delay` and `max_delay` come into play. -Lets build the count data and the number of time steps in the data. +Let’s build the count data and the number of time steps in the data. ``` r @@ -182,7 +192,7 @@ time_steps <- build_tstep(ps_count_example) #> ✔ Successfully built the number of time steps 8 ``` -Now we can create the number of inviduals and the number of transmisisons expected within the time bin. +Now we can create the number of individuals and the number of transmissions expected within the time bin. ``` r nind <- length(unique(ps_det_example$tag_serial_no)) @@ -192,13 +202,14 @@ ntrans <- build_ntrans(ps_det_example) #> "mean delay". ``` -We can finally move on to running the model, however, we need to understand a bit more about Bayesian frameworks before we proceed. +We can finally move on to running the model; however, we need to understand a bit more about Bayesian frameworks before we proceed. ## Bayesian Analysis -We can now estimate the center of activity for lake trout in Parry Sound, Lake Huron. +We can now estimate the center of activity for lake trout in Parry Sound +which is a large embayment of Georgian Bay, Lake Huron. -There are a few things to know about running a Bayesian analysis, I suggest reading these resources: +There are a few things to know about running a Bayesian analysis, we suggest reading these resources: 1. [Basics of Bayesian Statistics - Book](https://statswithr.github.io/book/) 2. [A Student’s Guide to Bayesian Statistics by Ben Lambert](https://sites.math.rutgers.edu/~zeilberg/EM20/Lambert.pdf) @@ -209,14 +220,34 @@ There are a few things to know about running a Bayesian analysis, I suggest read ### 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. +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). + +$$ +\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. ### Model convergence -It is important to always run the model with at least 2 chains. If the model does not converge you can try to increase the following: +It is important to always run the model with at least 2 chains. If the model does not converge you can try to increase or decreasing the following: -1. The amount of samples that are discarded; this can be controlled by the argument `warmup`. +1. The number of samples that are discarded; this can be controlled by the argument `warmup`. 2. The number of iterative samples retained; this can be controlled by the argument `iter`. @@ -225,10 +256,10 @@ It is important to always run the model with at least 2 chains. If the model doe 4. The `adapt_delta` value using `control = list(adapt_delta = 0.95)`. -When assessing the model we want $\hat R$ to be 1 or within 0.05 of 1, which indicates that the variance among and within chains are equal (see [{rstan} documentation on $\hat R$](https://mc-stan.org/rstan/reference/Rhat.html)), a high value for effective sample size (ESS), trace plots to look "grassy" or "caterpillar like," and posterior distributions to look relatively normal. +When assessing the model, we want $\hat R$ to be 1 or within 0.05 of 1, which indicates that the variance among and within chains are equal (see [{rstan} documentation on $\hat R$](https://mc-stan.org/rstan/reference/Rhat.html)), a high value for effective sample size (ESS), trace plots to look "grassy" or "caterpillar like," and posterior distributions to look relatively normal. ## Model -We can now run the standard point process/detection probablity model. +We can now run the standard point process/detection probability model. ``` r m <- COA_Standard( @@ -247,8 +278,8 @@ m <- COA_Standard( #> #> SAMPLING FOR MODEL 'COA_Standard_gaussian' NOW (CHAIN 1). #> Chain 1: -#> Chain 1: Gradient evaluation took 0.001044 seconds -#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 10.44 seconds. +#> Chain 1: Gradient evaluation took 0.001088 seconds +#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 10.88 seconds. #> Chain 1: Adjust your expectations accordingly! #> Chain 1: #> Chain 1: @@ -265,15 +296,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.552 seconds (Warm-up) -#> Chain 1: 14.249 seconds (Sampling) -#> Chain 1: 26.801 seconds (Total) +#> Chain 1: Elapsed Time: 13.128 seconds (Warm-up) +#> Chain 1: 10.399 seconds (Sampling) +#> Chain 1: 23.527 seconds (Total) #> Chain 1: #> #> SAMPLING FOR MODEL 'COA_Standard_gaussian' NOW (CHAIN 2). #> Chain 2: -#> Chain 2: Gradient evaluation took 0.001059 seconds -#> Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 10.59 seconds. +#> Chain 2: Gradient evaluation took 0.001036 seconds +#> Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 10.36 seconds. #> Chain 2: Adjust your expectations accordingly! #> Chain 2: #> Chain 2: @@ -290,20 +321,21 @@ m <- COA_Standard( #> Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling) #> Chain 2: Iteration: 2000 / 2000 [100%] (Sampling) #> Chain 2: -#> Chain 2: Elapsed Time: 12.656 seconds (Warm-up) -#> Chain 2: 8.811 seconds (Sampling) -#> Chain 2: 21.467 seconds (Total) -#> Chain 2: -#> Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable. -#> Running the chains for more iterations may help. See -#> https://mc-stan.org/misc/warnings.html#bulk-ess +#> Chain 2: Elapsed Time: 12.355 seconds (Warm-up) +#> Chain 2: 12.772 seconds (Sampling) +#> Chain 2: 25.127 seconds (Total) +#> Chain 2: #> warmup sample -#> chain:1 12.552 14.249 -#> chain:2 12.656 8.811 +#> chain:1 13.128 10.399 +#> chain:2 12.355 12.772 ``` -We can inspect the model object with the first object containing the Stan model object (accessible via `m$model`, which will allow you to use `rstan` plotting tools and diagnostic plots - see rstan documentation for details). +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 +Are a bunch of other objects in the objected created. These objects have pulled important information from the Stan model. + + ``` r summary(m) @@ -318,8 +350,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 distrubitons to ensure the model converged proeprly. +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. ``` r @@ -342,11 +373,11 @@ stan_dens( We can see our trace plots look good and that the posterior distributions of our paramaters look good. -Next lets look at our latent variables which are `sx` and `sy` or the estiamted locations of the fish based on -the detection probablity. For large/lengthy time periods it will become cumbersome to evaluate the posteriors +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]")) ``` @@ -364,32 +395,33 @@ stan_dens( ``` ![plot of chunk trace plots and post dist latent](figure/trace plots and post dist latent-2.png) -Again we can see that everything looks good and our model has converted. +Again, we can see that everything looks good and our model has converged. + +Now moving on to the other elements in the outputted object from `COA_*()`. -Now moving on to the other elements in the outputed object from `COA_*()`. -The second element is a table of parameter estimates, `p0` which is the estimated detection probality at distance of 0 m and the standard deviation (`sigma`). Their mean, standard error of the mean, and the associated quantiles from the posterior distriubtions. The table also includes the effective sample size and the Rhat statistic (which should be between 0.95 and 1.05). +The second element is a table of parameter estimates, `p0` which is the estimated detection probability at distance of 0 m and $\sigma$ which is units of distance and characterizes the effective spatial range of detection: it is the standard deviation of the Gaussian or logistic decay in detection log-odds with distance from the activity center. A larger $\sigma$ corresponds to a more slow-decaying (longer-range) detection function (`sigma`). The table contains their mean, standard error of the mean, and the associated quantiles from the posterior distributions. The table also includes the effective sample size and the $\hat R$ statistic (which should be between 0.95 and 1.05). ``` r m$summary -#> mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat -#> p0 0.5019460 0.0008856739 0.01633549 0.4723767 0.4899777 0.5012662 0.5133962 0.5319465 340.1860 0.9987742 -#> sigma 0.9949453 0.0014203483 0.02650908 0.9404368 0.9776532 0.9957883 1.0124727 1.0467824 348.3371 1.0023169 +#> mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat +#> p0 0.5021174 0.0009529285 0.01918709 0.4684838 0.4877082 0.5017623 0.5153259 0.5396254 405.4131 1.000994 +#> sigma 0.9952465 0.0014471561 0.02871756 0.9363534 0.9754131 0.9963352 1.0155073 1.0547302 393.7897 1.001954 ``` -We can see that the mean detection probablity at a distance of 0 m is 0.5 or 50% and that the model converged well for this paramaters. +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. -The thrid element returned, is the time required to run the model in minutes. Note that Stan will automatically detect and use multiple cores. If the computer used to run this has multiple cores, the time returned will be longer than the actual run time. This is because it will sum the time for each core. To return the realized run time, divide `fit$time` by the number of cores. +The thrid element returned, is the time required to run the model in minutes. Note that Stan will automatically detect and use multiple cores. If the computer used to run this has multiple cores, the time returned will be longer than the actual run time. This is because it will sum the time for each core. To return the realized run time, divide `m$time` by the number of cores. ``` r m$time -#> [1] 0.8044667 +#> [1] 0.8109 ``` -The fouth element returned, is a summary of the posterior draws. Returned is the `variable`, `median`, and the 2.5 and 97.5% quartiles (e.g., `q2.5` and `q97.5`). This summary `data.frame` provides you the ability to inspect all elements, however, most of the time a you will be using this information that has been filtered and displayed in a more user frinedly manner in futher objects returned by `COA_*()` functions. +The fourth element returned, is a summary of the posterior draws. Returned is the `variable`, `median`, and the 2.5 and 97.5% quartiles (e.g., `q2.5` and `q97.5`). This summary `data.frame` provides you the ability to inspect all elements, however, most of the time you will be using this data in a filtered and transformed manner created in other objects in the model output. These objects are more user friendly for plotting and understanding the results. ``` r @@ -397,20 +429,20 @@ m$summary_draws #> # A tibble: 21 × 4 #> variable median q2.5 q97.5 #> -#> 1 alpha0 0.00506 -0.111 0.128 -#> 2 alpha1 0.504 0.456 0.565 -#> 3 sx[1,1] -2.97 -3.30 -2.70 -#> 4 sx[1,2] -2.95 -3.33 -2.66 -#> 5 sx[1,3] -2.13 -2.34 -1.93 -#> 6 sx[1,4] -2.26 -2.51 -2.03 -#> 7 sx[1,5] -2.77 -3.07 -2.50 -#> 8 sx[1,6] -2.30 -2.53 -2.08 -#> 9 sx[1,7] -2.29 -2.49 -2.08 -#> 10 sx[1,8] -3.01 -3.27 -2.72 +#> 1 alpha0 0.00705 -0.126 0.159 +#> 2 alpha1 0.504 0.449 0.570 +#> 3 sx[1,1] -2.97 -3.23 -2.68 +#> 4 sx[1,2] -2.94 -3.32 -2.61 +#> 5 sx[1,3] -2.13 -2.35 -1.92 +#> 6 sx[1,4] -2.27 -2.50 -2.05 +#> 7 sx[1,5] -2.79 -3.15 -2.49 +#> 8 sx[1,6] -2.30 -2.53 -2.06 +#> 9 sx[1,7] -2.29 -2.50 -2.08 +#> 10 sx[1,8] -2.98 -3.26 -2.73 #> # ℹ 11 more rows ``` -The fourth returns `data.frame` of the median values for the latent variables `sx` and `sy` which are the estimated location of the fish. The median values can be viewed as the center of activity while the posterior distibution is the uncertainity around that estimate. The `data.frame` has a few important compoents `ind` is the index number of the invidiaul, `time` is the index number for a given time bin, `x` and `y` are the median values for a given individual and time, while `x_*` and `y_*` are the 2.5 and 97.5% quantiles for those posteriors, this effectivively is viewed as the 95% credible interval (Bayesian version of a confidence interval) for each coordinate. +The fifth element returned, is `data.frame` of the median values for the latent variables `sx` and `sy` which are the estimated location of the fish. The median values can be viewed as the center of activity while the posterior distribution is the overall probability of occurrence around that estimate. The `data.frame` has a few important components `ind` is the index number of the individual, `time` is the index number for a given time bin, `x` and `y` are the median values for a given individual and time, while `x_*` and `y_*` are the 2.5 and 97.5% quantiles for those posteriors, this effectively is viewed as the 95% credible interval, (i.e., Bayesian version of a confidence interval) for each coordinate. ``` r @@ -418,68 +450,350 @@ m$coas #> # A tibble: 8 × 8 #> ind time x y x_lower x_upper y_lower y_upper #> -#> 1 1 1 -2.97 -0.510 -3.30 -2.70 -0.832 -0.181 -#> 2 1 2 -2.95 -0.311 -3.33 -2.66 -0.713 0.107 -#> 3 1 3 -2.13 0.0464 -2.34 -1.93 -0.174 0.275 -#> 4 1 4 -2.26 0.152 -2.51 -2.03 -0.101 0.369 -#> 5 1 5 -2.77 -0.247 -3.07 -2.50 -0.645 0.0998 -#> 6 1 6 -2.30 -0.730 -2.53 -2.08 -1.01 -0.437 -#> 7 1 7 -2.29 -0.0731 -2.49 -2.08 -0.308 0.144 -#> 8 1 8 -3.01 -0.474 -3.27 -2.72 -0.760 -0.174 +#> 1 1 1 -2.97 -0.515 -3.23 -2.68 -0.903 -0.163 +#> 2 1 2 -2.94 -0.326 -3.32 -2.61 -0.744 0.0513 +#> 3 1 3 -2.13 0.0554 -2.35 -1.92 -0.192 0.249 +#> 4 1 4 -2.27 0.145 -2.50 -2.05 -0.0818 0.383 +#> 5 1 5 -2.79 -0.285 -3.15 -2.49 -0.602 0.0360 +#> 6 1 6 -2.30 -0.728 -2.53 -2.06 -1.06 -0.445 +#> 7 1 7 -2.29 -0.0921 -2.50 -2.08 -0.295 0.129 +#> 8 1 8 -2.98 -0.476 -3.26 -2.73 -0.748 -0.160 ``` - -The fith element (`m$all_estimates`) contains the posterior draws from each non-warm-up iteration from all chains. If you're not familiar with Bayesian terminology all you need to know is that this contains disturbtions for each paramter and latent variable for each individual in each time step. It is unlikely that you will use this object instead we have created futher objects below that you are more likely to use. +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 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.086 0.57 -2.6 -2.7 -2.0 -2.2 -2.6 -2.2 -#> 2 -0.138 0.49 -3.0 -2.9 -2.1 -2.3 -2.6 -2.3 -#> 3 -0.096 0.43 -2.9 -3.1 -2.2 -2.2 -3.0 -2.8 -#> 4 -0.024 0.51 -3.0 -3.0 -2.3 -2.3 -2.8 -2.5 -#> 5 0.032 0.48 -2.9 -2.9 -2.0 -2.1 -2.7 -2.4 -#> 6 -0.056 0.51 -2.9 -2.8 -2.1 -2.1 -2.8 -2.3 -#> 7 -0.077 0.46 -3.0 -3.0 -2.3 -2.2 -2.7 -2.3 -#> 8 0.043 0.55 -2.9 -2.8 -2.2 -2.1 -2.5 -2.2 -#> 9 -0.047 0.51 -3.0 -2.7 -2.2 -2.3 -3.0 -2.1 -#> 10 -0.036 0.47 -3.0 -3.1 -2.0 -2.4 -2.9 -2.3 +#> alpha0 alpha1 sx[1,1] sx[1,2] sx[1,3] sx[1,4] sx[1,5] sx[1,6] +#> 1 -0.11145 0.45 -3.2 -3.1 -2.2 -2.2 -2.7 -2.4 +#> 2 0.10130 0.49 -2.7 -3.0 -2.1 -2.2 -2.8 -2.3 +#> 3 0.08593 0.53 -3.0 -2.8 -2.2 -2.5 -3.0 -2.4 +#> 4 0.00502 0.48 -2.8 -3.0 -2.0 -2.2 -2.9 -2.0 +#> 5 -0.05169 0.48 -3.1 -3.0 -2.2 -2.5 -2.8 -2.2 +#> 6 -0.00041 0.52 -3.1 -2.8 -2.2 -2.4 -2.7 -2.1 +#> 7 -0.00792 0.48 -3.1 -3.1 -2.2 -2.4 -2.7 -2.2 +#> 8 -0.03658 0.47 -3.4 -3.4 -2.0 -2.4 -3.0 -2.5 +#> 9 0.00985 0.53 -2.9 -3.0 -2.2 -2.4 -2.7 -2.3 +#> 10 0.01182 0.54 -2.7 -2.8 -2.0 -2.2 -3.1 -2.0 #> # ... with 390 more draws, and 13 more variables #> # ... hidden reserved variables {'.chain', '.iteration', '.draw'} ``` -The sixth element (`m$loc_draws`), contains an extracted and transformed posterior draws -for the latent variables `sx` and `sy`. This object we will use to plot our estimates of -centers of activity. There are several other columns besides, the draw (i.e, `x`, `y`), the -`fish`/ind index number and the `time` index number. We have `.chain`, `.itteration`, `.draw` which identifies the chain, iteration, and draw that posterior draw came from. We also -have the column `lp__` which is the unnormalized log posterior density of your model, used for diagnosing sampling efficiency, estimating approximations, and comparing models. -We are going to inspect this element and then plot it. +The seventh element returned is a `data.frame` that contains the extracted and transformed posterior draws for the latent variables `sx` and `sy`. This object we will use to plot our estimates of centers of activity. There are several other columns besides, the draws (i.e, `x`, `y`), the `fish`/ind index number and the `time` index number. We have `.chain`, `.iteration`, `.draw` which identifies the chain, iteration, and draw that posterior draw came from. We also have the column `lp__` which is the unnormalized log posterior density of your model, used for diagnosing sampling efficiency, estimating approximations, and comparing models. + +Later on we will plot this object. + ``` r loc_draws <- m$loc_draws loc_draws #> # A tibble: 3,200 × 8 -#> .chain .iteration .draw lp__ fish time x y -#> -#> 1 1 1 1 -1282. 1 1 -2.65 -0.518 -#> 2 1 1 1 -1282. 1 2 -2.71 -0.0455 -#> 3 1 1 1 -1282. 1 3 -1.97 -0.0654 -#> 4 1 1 1 -1282. 1 4 -2.21 0.0812 -#> 5 1 1 1 -1282. 1 5 -2.63 -0.319 -#> 6 1 1 1 -1282. 1 6 -2.25 -0.881 -#> 7 1 1 1 -1282. 1 7 -2.09 -0.0597 -#> 8 1 1 1 -1282. 1 8 -3.00 -0.385 -#> 9 1 2 2 -1281. 1 1 -2.96 -0.355 -#> 10 1 2 2 -1281. 1 2 -2.89 -0.489 +#> .chain .iteration .draw lp__ fish time x y +#> +#> 1 1 1 1 -1281. 1 1 -3.25 -0.609 +#> 2 1 1 1 -1281. 1 2 -3.08 -0.391 +#> 3 1 1 1 -1281. 1 3 -2.20 0.00656 +#> 4 1 1 1 -1281. 1 4 -2.23 0.122 +#> 5 1 1 1 -1281. 1 5 -2.75 -0.342 +#> 6 1 1 1 -1281. 1 6 -2.43 -0.610 +#> 7 1 1 1 -1281. 1 7 -2.31 -0.202 +#> 8 1 1 1 -1281. 1 8 -3.08 -0.503 +#> 9 1 2 2 -1289. 1 1 -2.75 -0.698 +#> 10 1 2 2 -1289. 1 2 -2.96 -0.546 #> # ℹ 3,190 more rows ``` -Next 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 that we plot this onto so we can understand where the fish was during each time period. +The eighth element returned, is a `data.frame` that contains the contains the extracted and transformed posterior draws for the model parameters `alpha0`, `alpah1`, `p0`, and `sigma`. There are several other columns besides, the draws, the `fish`/ind index number and the `time` index number. We have `.chain`, `.iteration`, `.draw` which identifies the chain, iteration, and draw that posterior draw came from. We also have the column `lp__` which is the unnormalized log posterior density of your model, used for diagnosing sampling efficiency, estimating approximations, and comparing models. + +Later on we will plot this object. + +``` r +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 -1281. 1 1 -0.111 0.446 0.472 1.06 +#> 2 1 1 1 -1281. 1 2 -0.111 0.446 0.472 1.06 +#> 3 1 1 1 -1281. 1 3 -0.111 0.446 0.472 1.06 +#> 4 1 1 1 -1281. 1 4 -0.111 0.446 0.472 1.06 +#> 5 1 1 1 -1281. 1 5 -0.111 0.446 0.472 1.06 +#> 6 1 1 1 -1281. 1 6 -0.111 0.446 0.472 1.06 +#> 7 1 1 1 -1281. 1 7 -0.111 0.446 0.472 1.06 +#> 8 1 1 1 -1281. 1 8 -0.111 0.446 0.472 1.06 +#> 9 1 1 1 -1281. 1 1 -0.111 0.446 0.472 1.06 +#> 10 1 1 1 -1281. 1 2 -0.111 0.446 0.472 1.06 +#> # ℹ 6,390 more rows +``` + +The ninth element returned is a `list` that contains generated quantities for `y` from the model. This `y` is different than `sy` and is the estimated counts of detections for each individual at each time step at each receiver. The object `yrep` which is a 2-dimensional array that contains 10 draws (`ndraws` argument controls the number of draws here) and 640 posterior counts of detections the tag number, receiver number and time bin. We will assess to see how well this matches the existing data's distribution as using [posterior probability](https://en.wikipedia.org/wiki/Posterior_probability) which is a type of conditional probability that results from updating the prior probability with information summarized by the likelihood. + + +``` r +yrep <- m$generated_quantities +yrep +#> $yrep +#> tag_1_rec_1_time_1 tag_1_rec_2_time_1 tag_1_rec_3_time_1 tag_1_rec_4_time_1 tag_1_rec_5_time_1 tag_1_rec_6_time_1 +#> yrep_1 0 3 3 9 3 8 +#> tag_1_rec_7_time_1 tag_1_rec_8_time_1 tag_1_rec_9_time_1 tag_1_rec_10_time_1 tag_1_rec_11_time_1 tag_1_rec_12_time_1 +#> yrep_1 10 4 1 1 0 0 +#> tag_1_rec_13_time_1 tag_1_rec_14_time_1 tag_1_rec_15_time_1 tag_1_rec_16_time_1 tag_1_rec_17_time_1 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_18_time_1 tag_1_rec_19_time_1 tag_1_rec_20_time_1 tag_1_rec_21_time_1 tag_1_rec_22_time_1 +#> yrep_1 1 3 3 3 0 +#> tag_1_rec_23_time_1 tag_1_rec_24_time_1 tag_1_rec_25_time_1 tag_1_rec_26_time_1 tag_1_rec_27_time_1 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_28_time_1 tag_1_rec_29_time_1 tag_1_rec_30_time_1 tag_1_rec_31_time_1 tag_1_rec_32_time_1 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_33_time_1 tag_1_rec_34_time_1 tag_1_rec_35_time_1 tag_1_rec_36_time_1 tag_1_rec_37_time_1 +#> yrep_1 0 1 0 0 0 +#> tag_1_rec_38_time_1 tag_1_rec_39_time_1 tag_1_rec_40_time_1 tag_1_rec_41_time_1 tag_1_rec_42_time_1 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_43_time_1 tag_1_rec_44_time_1 tag_1_rec_45_time_1 tag_1_rec_46_time_1 tag_1_rec_47_time_1 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_48_time_1 tag_1_rec_49_time_1 tag_1_rec_50_time_1 tag_1_rec_51_time_1 tag_1_rec_52_time_1 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_53_time_1 tag_1_rec_54_time_1 tag_1_rec_55_time_1 tag_1_rec_56_time_1 tag_1_rec_57_time_1 +#> yrep_1 0 0 1 0 0 +#> tag_1_rec_58_time_1 tag_1_rec_59_time_1 tag_1_rec_60_time_1 tag_1_rec_61_time_1 tag_1_rec_62_time_1 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_63_time_1 tag_1_rec_64_time_1 tag_1_rec_65_time_1 tag_1_rec_66_time_1 tag_1_rec_67_time_1 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_68_time_1 tag_1_rec_69_time_1 tag_1_rec_70_time_1 tag_1_rec_71_time_1 tag_1_rec_72_time_1 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_73_time_1 tag_1_rec_74_time_1 tag_1_rec_75_time_1 tag_1_rec_76_time_1 tag_1_rec_77_time_1 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_78_time_1 tag_1_rec_79_time_1 tag_1_rec_80_time_1 tag_1_rec_1_time_2 tag_1_rec_2_time_2 tag_1_rec_3_time_2 +#> yrep_1 0 0 0 0 0 5 +#> tag_1_rec_4_time_2 tag_1_rec_5_time_2 tag_1_rec_6_time_2 tag_1_rec_7_time_2 tag_1_rec_8_time_2 tag_1_rec_9_time_2 +#> yrep_1 5 6 6 10 2 0 +#> tag_1_rec_10_time_2 tag_1_rec_11_time_2 tag_1_rec_12_time_2 tag_1_rec_13_time_2 tag_1_rec_14_time_2 +#> yrep_1 1 0 0 0 0 +#> tag_1_rec_15_time_2 tag_1_rec_16_time_2 tag_1_rec_17_time_2 tag_1_rec_18_time_2 tag_1_rec_19_time_2 +#> yrep_1 0 0 0 1 2 +#> tag_1_rec_20_time_2 tag_1_rec_21_time_2 tag_1_rec_22_time_2 tag_1_rec_23_time_2 tag_1_rec_24_time_2 +#> yrep_1 3 1 4 0 1 +#> tag_1_rec_25_time_2 tag_1_rec_26_time_2 tag_1_rec_27_time_2 tag_1_rec_28_time_2 tag_1_rec_29_time_2 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_30_time_2 tag_1_rec_31_time_2 tag_1_rec_32_time_2 tag_1_rec_33_time_2 tag_1_rec_34_time_2 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_35_time_2 tag_1_rec_36_time_2 tag_1_rec_37_time_2 tag_1_rec_38_time_2 tag_1_rec_39_time_2 +#> yrep_1 0 0 1 0 0 +#> tag_1_rec_40_time_2 tag_1_rec_41_time_2 tag_1_rec_42_time_2 tag_1_rec_43_time_2 tag_1_rec_44_time_2 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_45_time_2 tag_1_rec_46_time_2 tag_1_rec_47_time_2 tag_1_rec_48_time_2 tag_1_rec_49_time_2 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_50_time_2 tag_1_rec_51_time_2 tag_1_rec_52_time_2 tag_1_rec_53_time_2 tag_1_rec_54_time_2 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_55_time_2 tag_1_rec_56_time_2 tag_1_rec_57_time_2 tag_1_rec_58_time_2 tag_1_rec_59_time_2 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_60_time_2 tag_1_rec_61_time_2 tag_1_rec_62_time_2 tag_1_rec_63_time_2 tag_1_rec_64_time_2 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_65_time_2 tag_1_rec_66_time_2 tag_1_rec_67_time_2 tag_1_rec_68_time_2 tag_1_rec_69_time_2 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_70_time_2 tag_1_rec_71_time_2 tag_1_rec_72_time_2 tag_1_rec_73_time_2 tag_1_rec_74_time_2 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_75_time_2 tag_1_rec_76_time_2 tag_1_rec_77_time_2 tag_1_rec_78_time_2 tag_1_rec_79_time_2 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_80_time_2 tag_1_rec_1_time_3 tag_1_rec_2_time_3 tag_1_rec_3_time_3 tag_1_rec_4_time_3 tag_1_rec_5_time_3 +#> yrep_1 0 0 0 7 9 2 +#> tag_1_rec_6_time_3 tag_1_rec_7_time_3 tag_1_rec_8_time_3 tag_1_rec_9_time_3 tag_1_rec_10_time_3 tag_1_rec_11_time_3 +#> yrep_1 2 5 7 7 1 0 +#> tag_1_rec_12_time_3 tag_1_rec_13_time_3 tag_1_rec_14_time_3 tag_1_rec_15_time_3 tag_1_rec_16_time_3 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_17_time_3 tag_1_rec_18_time_3 tag_1_rec_19_time_3 tag_1_rec_20_time_3 tag_1_rec_21_time_3 +#> yrep_1 1 4 4 8 2 +#> tag_1_rec_22_time_3 tag_1_rec_23_time_3 tag_1_rec_24_time_3 tag_1_rec_25_time_3 tag_1_rec_26_time_3 +#> yrep_1 5 0 4 1 0 +#> tag_1_rec_27_time_3 tag_1_rec_28_time_3 tag_1_rec_29_time_3 tag_1_rec_30_time_3 tag_1_rec_31_time_3 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_32_time_3 tag_1_rec_33_time_3 tag_1_rec_34_time_3 tag_1_rec_35_time_3 tag_1_rec_36_time_3 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_37_time_3 tag_1_rec_38_time_3 tag_1_rec_39_time_3 tag_1_rec_40_time_3 tag_1_rec_41_time_3 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_42_time_3 tag_1_rec_43_time_3 tag_1_rec_44_time_3 tag_1_rec_45_time_3 tag_1_rec_46_time_3 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_47_time_3 tag_1_rec_48_time_3 tag_1_rec_49_time_3 tag_1_rec_50_time_3 tag_1_rec_51_time_3 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_52_time_3 tag_1_rec_53_time_3 tag_1_rec_54_time_3 tag_1_rec_55_time_3 tag_1_rec_56_time_3 +#> yrep_1 0 0 3 2 1 +#> tag_1_rec_57_time_3 tag_1_rec_58_time_3 tag_1_rec_59_time_3 tag_1_rec_60_time_3 tag_1_rec_61_time_3 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_62_time_3 tag_1_rec_63_time_3 tag_1_rec_64_time_3 tag_1_rec_65_time_3 tag_1_rec_66_time_3 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_67_time_3 tag_1_rec_68_time_3 tag_1_rec_69_time_3 tag_1_rec_70_time_3 tag_1_rec_71_time_3 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_72_time_3 tag_1_rec_73_time_3 tag_1_rec_74_time_3 tag_1_rec_75_time_3 tag_1_rec_76_time_3 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_77_time_3 tag_1_rec_78_time_3 tag_1_rec_79_time_3 tag_1_rec_80_time_3 tag_1_rec_1_time_4 tag_1_rec_2_time_4 +#> yrep_1 0 0 0 0 2 1 +#> tag_1_rec_3_time_4 tag_1_rec_4_time_4 tag_1_rec_5_time_4 tag_1_rec_6_time_4 tag_1_rec_7_time_4 tag_1_rec_8_time_4 +#> yrep_1 6 4 1 3 3 9 +#> tag_1_rec_9_time_4 tag_1_rec_10_time_4 tag_1_rec_11_time_4 tag_1_rec_12_time_4 tag_1_rec_13_time_4 +#> yrep_1 3 1 0 0 0 +#> tag_1_rec_14_time_4 tag_1_rec_15_time_4 tag_1_rec_16_time_4 tag_1_rec_17_time_4 tag_1_rec_18_time_4 +#> yrep_1 0 0 0 3 4 +#> tag_1_rec_19_time_4 tag_1_rec_20_time_4 tag_1_rec_21_time_4 tag_1_rec_22_time_4 tag_1_rec_23_time_4 +#> yrep_1 6 10 2 1 3 +#> tag_1_rec_24_time_4 tag_1_rec_25_time_4 tag_1_rec_26_time_4 tag_1_rec_27_time_4 tag_1_rec_28_time_4 +#> yrep_1 1 1 0 0 0 +#> tag_1_rec_29_time_4 tag_1_rec_30_time_4 tag_1_rec_31_time_4 tag_1_rec_32_time_4 tag_1_rec_33_time_4 +#> yrep_1 0 0 0 1 0 +#> tag_1_rec_34_time_4 tag_1_rec_35_time_4 tag_1_rec_36_time_4 tag_1_rec_37_time_4 tag_1_rec_38_time_4 +#> yrep_1 1 1 0 0 0 +#> tag_1_rec_39_time_4 tag_1_rec_40_time_4 tag_1_rec_41_time_4 tag_1_rec_42_time_4 tag_1_rec_43_time_4 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_44_time_4 tag_1_rec_45_time_4 tag_1_rec_46_time_4 tag_1_rec_47_time_4 tag_1_rec_48_time_4 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_49_time_4 tag_1_rec_50_time_4 tag_1_rec_51_time_4 tag_1_rec_52_time_4 tag_1_rec_53_time_4 +#> yrep_1 0 0 0 1 0 +#> tag_1_rec_54_time_4 tag_1_rec_55_time_4 tag_1_rec_56_time_4 tag_1_rec_57_time_4 tag_1_rec_58_time_4 +#> yrep_1 1 0 0 0 0 +#> tag_1_rec_59_time_4 tag_1_rec_60_time_4 tag_1_rec_61_time_4 tag_1_rec_62_time_4 tag_1_rec_63_time_4 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_64_time_4 tag_1_rec_65_time_4 tag_1_rec_66_time_4 tag_1_rec_67_time_4 tag_1_rec_68_time_4 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_69_time_4 tag_1_rec_70_time_4 tag_1_rec_71_time_4 tag_1_rec_72_time_4 tag_1_rec_73_time_4 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_74_time_4 tag_1_rec_75_time_4 tag_1_rec_76_time_4 tag_1_rec_77_time_4 tag_1_rec_78_time_4 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_79_time_4 tag_1_rec_80_time_4 tag_1_rec_1_time_5 tag_1_rec_2_time_5 tag_1_rec_3_time_5 tag_1_rec_4_time_5 +#> yrep_1 0 0 1 7 4 8 +#> tag_1_rec_5_time_5 tag_1_rec_6_time_5 tag_1_rec_7_time_5 tag_1_rec_8_time_5 tag_1_rec_9_time_5 tag_1_rec_10_time_5 +#> yrep_1 3 3 8 5 1 0 +#> tag_1_rec_11_time_5 tag_1_rec_12_time_5 tag_1_rec_13_time_5 tag_1_rec_14_time_5 tag_1_rec_15_time_5 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_16_time_5 tag_1_rec_17_time_5 tag_1_rec_18_time_5 tag_1_rec_19_time_5 tag_1_rec_20_time_5 +#> yrep_1 0 0 2 1 3 +#> tag_1_rec_21_time_5 tag_1_rec_22_time_5 tag_1_rec_23_time_5 tag_1_rec_24_time_5 tag_1_rec_25_time_5 +#> yrep_1 0 1 1 0 0 +#> tag_1_rec_26_time_5 tag_1_rec_27_time_5 tag_1_rec_28_time_5 tag_1_rec_29_time_5 tag_1_rec_30_time_5 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_31_time_5 tag_1_rec_32_time_5 tag_1_rec_33_time_5 tag_1_rec_34_time_5 tag_1_rec_35_time_5 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_36_time_5 tag_1_rec_37_time_5 tag_1_rec_38_time_5 tag_1_rec_39_time_5 tag_1_rec_40_time_5 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_41_time_5 tag_1_rec_42_time_5 tag_1_rec_43_time_5 tag_1_rec_44_time_5 tag_1_rec_45_time_5 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_46_time_5 tag_1_rec_47_time_5 tag_1_rec_48_time_5 tag_1_rec_49_time_5 tag_1_rec_50_time_5 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_51_time_5 tag_1_rec_52_time_5 tag_1_rec_53_time_5 tag_1_rec_54_time_5 tag_1_rec_55_time_5 +#> yrep_1 0 1 0 2 0 +#> tag_1_rec_56_time_5 tag_1_rec_57_time_5 tag_1_rec_58_time_5 tag_1_rec_59_time_5 tag_1_rec_60_time_5 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_61_time_5 tag_1_rec_62_time_5 tag_1_rec_63_time_5 tag_1_rec_64_time_5 tag_1_rec_65_time_5 +#> yrep_1 0 0 0 0 1 +#> tag_1_rec_66_time_5 tag_1_rec_67_time_5 tag_1_rec_68_time_5 tag_1_rec_69_time_5 tag_1_rec_70_time_5 +#> yrep_1 1 0 0 0 0 +#> tag_1_rec_71_time_5 tag_1_rec_72_time_5 tag_1_rec_73_time_5 tag_1_rec_74_time_5 tag_1_rec_75_time_5 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_76_time_5 tag_1_rec_77_time_5 tag_1_rec_78_time_5 tag_1_rec_79_time_5 tag_1_rec_80_time_5 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_1_time_6 tag_1_rec_2_time_6 tag_1_rec_3_time_6 tag_1_rec_4_time_6 tag_1_rec_5_time_6 tag_1_rec_6_time_6 +#> yrep_1 2 5 6 7 9 3 +#> tag_1_rec_7_time_6 tag_1_rec_8_time_6 tag_1_rec_9_time_6 tag_1_rec_10_time_6 tag_1_rec_11_time_6 tag_1_rec_12_time_6 +#> yrep_1 5 6 2 0 1 0 +#> tag_1_rec_13_time_6 tag_1_rec_14_time_6 tag_1_rec_15_time_6 tag_1_rec_16_time_6 tag_1_rec_17_time_6 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_18_time_6 tag_1_rec_19_time_6 tag_1_rec_20_time_6 tag_1_rec_21_time_6 tag_1_rec_22_time_6 +#> yrep_1 2 2 3 0 0 +#> tag_1_rec_23_time_6 tag_1_rec_24_time_6 tag_1_rec_25_time_6 tag_1_rec_26_time_6 tag_1_rec_27_time_6 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_28_time_6 tag_1_rec_29_time_6 tag_1_rec_30_time_6 tag_1_rec_31_time_6 tag_1_rec_32_time_6 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_33_time_6 tag_1_rec_34_time_6 tag_1_rec_35_time_6 tag_1_rec_36_time_6 tag_1_rec_37_time_6 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_38_time_6 tag_1_rec_39_time_6 tag_1_rec_40_time_6 tag_1_rec_41_time_6 tag_1_rec_42_time_6 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_43_time_6 tag_1_rec_44_time_6 tag_1_rec_45_time_6 tag_1_rec_46_time_6 tag_1_rec_47_time_6 +#> yrep_1 1 1 0 0 0 +#> tag_1_rec_48_time_6 tag_1_rec_49_time_6 tag_1_rec_50_time_6 tag_1_rec_51_time_6 tag_1_rec_52_time_6 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_53_time_6 tag_1_rec_54_time_6 tag_1_rec_55_time_6 tag_1_rec_56_time_6 tag_1_rec_57_time_6 +#> yrep_1 0 2 0 0 0 +#> tag_1_rec_58_time_6 tag_1_rec_59_time_6 tag_1_rec_60_time_6 tag_1_rec_61_time_6 tag_1_rec_62_time_6 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_63_time_6 tag_1_rec_64_time_6 tag_1_rec_65_time_6 tag_1_rec_66_time_6 tag_1_rec_67_time_6 +#> yrep_1 0 0 0 0 1 +#> tag_1_rec_68_time_6 tag_1_rec_69_time_6 tag_1_rec_70_time_6 tag_1_rec_71_time_6 tag_1_rec_72_time_6 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_73_time_6 tag_1_rec_74_time_6 tag_1_rec_75_time_6 tag_1_rec_76_time_6 tag_1_rec_77_time_6 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_78_time_6 tag_1_rec_79_time_6 tag_1_rec_80_time_6 tag_1_rec_1_time_7 tag_1_rec_2_time_7 tag_1_rec_3_time_7 +#> yrep_1 0 0 0 2 6 6 +#> tag_1_rec_4_time_7 tag_1_rec_5_time_7 tag_1_rec_6_time_7 tag_1_rec_7_time_7 tag_1_rec_8_time_7 tag_1_rec_9_time_7 +#> yrep_1 2 3 6 12 8 6 +#> tag_1_rec_10_time_7 tag_1_rec_11_time_7 tag_1_rec_12_time_7 tag_1_rec_13_time_7 tag_1_rec_14_time_7 +#> yrep_1 1 0 0 0 0 +#> tag_1_rec_15_time_7 tag_1_rec_16_time_7 tag_1_rec_17_time_7 tag_1_rec_18_time_7 tag_1_rec_19_time_7 +#> yrep_1 0 0 0 4 2 +#> tag_1_rec_20_time_7 tag_1_rec_21_time_7 tag_1_rec_22_time_7 tag_1_rec_23_time_7 tag_1_rec_24_time_7 +#> yrep_1 6 1 3 1 1 +#> tag_1_rec_25_time_7 tag_1_rec_26_time_7 tag_1_rec_27_time_7 tag_1_rec_28_time_7 tag_1_rec_29_time_7 +#> yrep_1 0 1 0 0 0 +#> tag_1_rec_30_time_7 tag_1_rec_31_time_7 tag_1_rec_32_time_7 tag_1_rec_33_time_7 tag_1_rec_34_time_7 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_35_time_7 tag_1_rec_36_time_7 tag_1_rec_37_time_7 tag_1_rec_38_time_7 tag_1_rec_39_time_7 +#> yrep_1 1 0 0 0 0 +#> tag_1_rec_40_time_7 tag_1_rec_41_time_7 tag_1_rec_42_time_7 tag_1_rec_43_time_7 tag_1_rec_44_time_7 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_45_time_7 tag_1_rec_46_time_7 tag_1_rec_47_time_7 tag_1_rec_48_time_7 tag_1_rec_49_time_7 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_50_time_7 tag_1_rec_51_time_7 tag_1_rec_52_time_7 tag_1_rec_53_time_7 tag_1_rec_54_time_7 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_55_time_7 tag_1_rec_56_time_7 tag_1_rec_57_time_7 tag_1_rec_58_time_7 tag_1_rec_59_time_7 +#> yrep_1 2 0 1 0 0 +#> tag_1_rec_60_time_7 tag_1_rec_61_time_7 tag_1_rec_62_time_7 tag_1_rec_63_time_7 tag_1_rec_64_time_7 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_65_time_7 tag_1_rec_66_time_7 tag_1_rec_67_time_7 tag_1_rec_68_time_7 tag_1_rec_69_time_7 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_70_time_7 tag_1_rec_71_time_7 tag_1_rec_72_time_7 tag_1_rec_73_time_7 tag_1_rec_74_time_7 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_75_time_7 tag_1_rec_76_time_7 tag_1_rec_77_time_7 tag_1_rec_78_time_7 tag_1_rec_79_time_7 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_80_time_7 tag_1_rec_1_time_8 tag_1_rec_2_time_8 tag_1_rec_3_time_8 tag_1_rec_4_time_8 tag_1_rec_5_time_8 +#> yrep_1 0 2 0 5 7 10 +#> tag_1_rec_6_time_8 tag_1_rec_7_time_8 tag_1_rec_8_time_8 tag_1_rec_9_time_8 tag_1_rec_10_time_8 tag_1_rec_11_time_8 +#> yrep_1 4 8 3 0 0 0 +#> tag_1_rec_12_time_8 tag_1_rec_13_time_8 tag_1_rec_14_time_8 tag_1_rec_15_time_8 tag_1_rec_16_time_8 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_17_time_8 tag_1_rec_18_time_8 tag_1_rec_19_time_8 tag_1_rec_20_time_8 tag_1_rec_21_time_8 +#> yrep_1 0 0 1 2 0 +#> tag_1_rec_22_time_8 tag_1_rec_23_time_8 tag_1_rec_24_time_8 tag_1_rec_25_time_8 tag_1_rec_26_time_8 +#> yrep_1 0 1 0 0 0 +#> tag_1_rec_27_time_8 tag_1_rec_28_time_8 tag_1_rec_29_time_8 tag_1_rec_30_time_8 tag_1_rec_31_time_8 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_32_time_8 tag_1_rec_33_time_8 tag_1_rec_34_time_8 tag_1_rec_35_time_8 tag_1_rec_36_time_8 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_37_time_8 tag_1_rec_38_time_8 tag_1_rec_39_time_8 tag_1_rec_40_time_8 tag_1_rec_41_time_8 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_42_time_8 tag_1_rec_43_time_8 tag_1_rec_44_time_8 tag_1_rec_45_time_8 tag_1_rec_46_time_8 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_47_time_8 tag_1_rec_48_time_8 tag_1_rec_49_time_8 tag_1_rec_50_time_8 tag_1_rec_51_time_8 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_52_time_8 tag_1_rec_53_time_8 tag_1_rec_54_time_8 tag_1_rec_55_time_8 tag_1_rec_56_time_8 +#> yrep_1 0 0 1 0 0 +#> tag_1_rec_57_time_8 tag_1_rec_58_time_8 tag_1_rec_59_time_8 tag_1_rec_60_time_8 tag_1_rec_61_time_8 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_62_time_8 tag_1_rec_63_time_8 tag_1_rec_64_time_8 tag_1_rec_65_time_8 tag_1_rec_66_time_8 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_67_time_8 tag_1_rec_68_time_8 tag_1_rec_69_time_8 tag_1_rec_70_time_8 tag_1_rec_71_time_8 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_72_time_8 tag_1_rec_73_time_8 tag_1_rec_74_time_8 tag_1_rec_75_time_8 tag_1_rec_76_time_8 +#> yrep_1 0 0 0 0 0 +#> tag_1_rec_77_time_8 tag_1_rec_78_time_8 tag_1_rec_79_time_8 tag_1_rec_80_time_8 +#> yrep_1 0 0 0 0 +#> [ reached 'max' / getOption("max.print") -- omitted 9 rows ] +``` + +## Plotting + +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. -First though 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 |> @@ -506,14 +820,82 @@ p <- ggplot() + p1 <- p + facet_wrap(~time) +``` + +Let's first look at the posterior densities for each time step + +``` r +p1 +``` + +![plot of chunk post densities timestep](figure/post densities timestep-1.png) + +Next let's look at the densities combined to gather a better understanding of the movement over the 8 timesteps. + + +``` r p ``` -![plot of chunk plot posteriors](figure/plot posteriors-1.png) +![plot of chunk post densities combined](figure/post densities combined-1.png) + +Next we can plot the posterior distributions for `alpha0`, `alpha1`, `p0`, and `sigma`. We need to first make `fish` and `time` a `character` and then to make plotting easier we need to make the `data.frame` in long format. + ``` r -p1 +param_draws$fish <- as.character(param_draws$fish) +param_draws$time <- as.character(param_draws$time) + +param_draws_long <- param_draws |> + pivot_longer( + cols = -c(".chain", ".iteration", ".draw", "lp__", "fish", "time"), + names_to = "param", + values_to = "est" + ) ``` -![plot of chunk plot posteriors](figure/plot posteriors-2.png) +We can now plot those posterior distributions of the model parameters. + + +``` r +p_param <- ggplot( + data = param_draws_long, + aes(x = time, y = est, fill = fish), +) + + geom_violin() + + facet_wrap(~param, scale = "free_y") + + scale_fill_manual(name = "Fish ID", values = "#2768F5") + + theme_bw() + + theme( + panel.grid = element_blank(), + strip.background = element_blank() + ) + + labs(x = "Time Bin", y = "Estimate") + +p_param +``` + +![plot of chunk plot 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. + +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. + + +``` 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 +``` + +![plot of chunk ppc dens](figure/ppc dens-1.png) + +We can see our predictive posterior check distribution for 10 draws closely lines up with +our observed detection counts indicating that the model fits the data well. + +Congratulations we have successfully run a Bayesian standard point processing/detection probability model for one Lake Trout for 8, 1 hour time bins (8 hrs) in Parry Sound, a large embayment of Georgian Bay, Lake Huron. We can take the results and write up how this individual behaved for a given time period based on the model results!