From 42e7533e5cf2b816f99e29f634a6537646fdd1cc Mon Sep 17 00:00:00 2001 From: Ken Kellner Date: Tue, 16 Feb 2021 21:42:34 -0500 Subject: Add new ubms spatial models post --- .gitignore | 1 + src/plos-one.csl | 16 +++ src/references.bib | 188 +++++++++++++++++++++++++++ src/ubms-spatial.Rmd | 347 ++++++++++++++++++++++++++++++++++++++++++++++++++ src/ubms-vignette.Rmd | 10 +- 5 files changed, 557 insertions(+), 5 deletions(-) create mode 100644 src/plos-one.csl create mode 100644 src/references.bib create mode 100644 src/ubms-spatial.Rmd diff --git a/.gitignore b/.gitignore index 43c785d..70e754d 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,4 @@ *.html *.xml data/* +*.Rds diff --git a/src/plos-one.csl b/src/plos-one.csl new file mode 100644 index 0000000..190a801 --- /dev/null +++ b/src/plos-one.csl @@ -0,0 +1,16 @@ + + diff --git a/src/references.bib b/src/references.bib new file mode 100644 index 0000000..d04362e --- /dev/null +++ b/src/references.bib @@ -0,0 +1,188 @@ +@article{Bates_2015, + author = {Bates, Douglas and Mächler, Martin and Bolker, Ben and Walker, Steve}, + title = {Fitting Linear Mixed-Effects Models Using lme4}, + journal = {Journal of Statistical Software}, + year = {2015}, + number = {1}, + volume = {67}, + doi = {10.18637/jss.v067.i01} +} + +@article{Carpenter_2017, + author = {Carpenter, Bob and Gelman, Andrew and Hoffman, Matthew D. and Lee, Daniel and Goodrich, Ben and Betancourt, Michael and Brubaker, Marcus and Guo, Jiqiang and Li, Peter and Riddell, Allen}, + title = {Stan: A Probabilistic Programming Language}, + journal = {Journal of Statistical Software}, + publisher = {Foundation for Open Access Statistic}, + year = {2017}, + number = {1}, + volume = {76}, + url = {http://dx.doi.org/10.18637/jss.v076.i01}, + doi = {10.18637/jss.v076.i01} +} + +@article{Chandler2011, + author = {Chandler, Richard B. and Royle, J. Andrew and King, David I.}, + title = {Inference about density and temporary emigration in unmarked populations}, + journal = {Ecology}, + year = {2011}, + number = {7}, + volume = {92}, + pages = {1429--1435}, + url = {http://dx.doi.org/10.1890/10-2433.1}, + doi = {10.1890/10-2433.1}, + file = {:Chandler2011b.pdf:PDF}, + keyword = {Chestnut-sided Warbler, Dendroica pensylvanica, N-mixture model, USA, White Mountain National Forest, detection probability, hierarchical, models, population density, spot-mapping, temporary emigration, unmarked populations} +} + +@article{Fiske_2011, + author = {Fiske, Ian and Chandler, Richard and others, }, + title = {Unmarked: an {R} package for fitting hierarchical models of wildlife occurrence and abundance}, + journal = {Journal of Statistical Software}, + year = {2011}, + number = {10}, + volume = {43}, + pages = {1--23}, + doi = {10.18637/jss.v043.i10} +} + +@book{Gelman2007, + author = {Gelman, Andrew and Hill, Jennifer}, + title = {Data analysis using regression and multilevel/hierarchical models}, + publisher = {Cambridge University Press}, + year = {2007}, + pages = {651}, + url = {http://dx.doi.org/10.2277/0521867061, http://dx.doi.org/10.2277/0521867061}, + doi = {10.2277/0521867061}, + address = {New York, NY}, + file = {:Gelman2007.pdf:PDF}, + isbn = {0521867061}, + pmid = {14341096} +} + +@article{Hoffman_2014, + author = {Hoffman, Matthew D and Gelman, Andrew}, + title = {The {N}o-{U}-{T}urn sampler: adaptively setting path lengths in {H}amiltonian {M}onte {C}arlo.}, + journal = {Journal of Machine Learning Research}, + year = {2014}, + number = {1}, + volume = {15}, + pages = {1593--1623} +} + +@article{MacKenzie_2002, + author = {MacKenzie, Darryl I. and Nichols, James D. and Lachman, Gideon B. and Droege, Sam and Royle, J. Andrew and Langtimm, Catherine A.}, + title = {Estimating site occupancy rates when detection probabilities are less than one}, + journal = {Ecology}, + publisher = {Wiley}, + year = {2002}, + month = {aug}, + number = {8}, + volume = {83}, + pages = {2248--2255}, + url = {http://dx.doi.org/10.1890/0012-9658(2002)083[2248:esorwd]2.0.co;2}, + doi = {10.1890/0012-9658(2002)083[2248:esorwd]2.0.co;2} +} + +@article{MacKenzie_2004a, + author = {MacKenzie, Darryl I. and Bailey, Larissa L.}, + title = {Assessing the fit of site-occupancy models}, + journal = {Journal of Agricultural, Biological, and Environmental Statistics}, + year = {2004}, + number = {3}, + volume = {9}, + pages = {300--318}, + doi = {10.1198/108571104x3361} +} + +@article{Schmid_2004, + author = {Schmid, Hans and Zbinden, Niklaus and Keller, Verena}, + title = {Überwachung der bestandsentwicklung häufiger brutvögel in der schweiz}, + journal = {Swiss Ornithological Institute Sempach Switzerland}, + year = {2004} +} + +@article{Vehtari_2017, + author = {Vehtari, Aki and Gelman, Andrew and Gabry, Jonah}, + title = {Practical {B}ayesian model evaluation using leave-one-out cross-validation and {WAIC}}, + journal = {Statistics and Computing}, + year = {2017}, + number = {5}, + volume = {27}, + pages = {1413--1432}, + doi = {10.1007/s11222-016-9696-4} +} + +@article{Wright_2019, + author = {Wright, Wilson J. and Irvine, Kathryn M. and Higgs, Megan D.}, + title = {Identifying occupancy model inadequacies: can residuals separately assess detection and presence?}, + journal = {Ecology}, + year = {2019}, + pages = {e02703}, + doi = {10.1002/ecy.2703} +} + +@article{Yackulic_2020, + author = {Yackulic, Charles B. and Dodrill, Michael and Dzul, Maria and Sanderlin, Jamie S. and Reid, Janice A.}, + title = {A need for speed in {B}ayesian population models: a practical guide to marginalizing and recovering discrete latent states}, + journal = {Ecological Applications}, + year = {2020}, + number = {5}, + volume = {30}, + doi = {10.1002/eap.2112} +} + +@article{Broms_2014, +author = {Broms, Kristin M. and Johnson, Devin S. and Altwegg, Res and Conquest, Loveday L.}, +title = {Spatial occupancy models applied to atlas data show Southern Ground Hornbills strongly depend on protected areas}, +journal = {Ecological Applications}, +volume = {24}, +number = {2}, +pages = {363-374}, +doi = {https://doi.org/10.1890/12-2151.1}, +url = {https://esajournals.onlinelibrary.wiley.com/doi/abs/10.1890/12-2151.1}, +year = {2014} +} + + +@article{Clark_2019, +author = {Clark, Allan E. and Altwegg, Res}, +title = {Efficient Bayesian analysis of occupancy models with logit link functions}, +journal = {Ecology and Evolution}, +volume = {9}, +number = {2}, +pages = {756-768}, +doi = {https://doi.org/10.1002/ece3.4850}, +url = {https://onlinelibrary.wiley.com/doi/abs/10.1002/ece3.4850}, +year = {2019} +} + + +@article{Hodges_2010, + title={Adding spatially-correlated errors can mess up the fixed effect you love}, + author={Hodges, James S and Reich, Brian J}, + journal={The American Statistician}, + volume={64}, + number={4}, + pages={325--334}, + year={2010}, + publisher={Taylor \& Francis} +} + +@article{Johnson_2013, +author = {Johnson, Devin S. and Conn, Paul B. and Hooten, Mevin B. and Ray, Justina C. and Pond, Bruce A.}, +title = {Spatial occupancy models for large data sets}, +journal = {Ecology}, +volume = {94}, +number = {4}, +pages = {801-808}, +doi = {https://doi.org/10.1890/12-0564.1}, +url = {https://esajournals.onlinelibrary.wiley.com/doi/abs/10.1890/12-0564.1}, +year = {2013} +} + +@book{Banerjee_2014, + title={Hierarchical modeling and analysis for spatial data}, + author={Banerjee, Sudipto and Carlin, Bradley P and Gelfand, Alan E}, + year={2014}, + publisher={CRC press} +} diff --git a/src/ubms-spatial.Rmd b/src/ubms-spatial.Rmd new file mode 100644 index 0000000..0d2d611 --- /dev/null +++ b/src/ubms-spatial.Rmd @@ -0,0 +1,347 @@ +--- +title: Spatial Models in ubms +date: 2021-02-15 +bibliography: references.bib +link-citations: yes +csl: plos-one.csl +output: + html_document: + css: style.css + highlight: pygments +--- + +# Introduction + +In many studies of animal and plant occurrence and abundance, sites may not be entirely independent. +For example, sites may be close together, or site size may be smaller than a typical animal's home range. +In these cases we expect spatial autocorrelation among nearby sites. +Ideally, this autocorrelation would be accounted for in models. + +A solution to this problem is so-called spatial occupancy models, which include spatial random effects. +A common approach to handling spatial autocorrelation is the use of an intrinsic conditional autoregressive (ICAR) random effect [@Banerjee_2014]. +The growing accessibility of Bayesian approaches has facilitated fitting these models. +However, fitting ICAR models is computationally intensive and the ICAR random effect may be correlated with fixed covariates included in the model, making inference difficult [@Hodges_2010]. +To solve these issues, Johnson et al. [-@Johnson_2013] extended restricted spatial regression (RSR) to occupancy models. +In an RSR model, the random effect is constructed in such a way that it is uncorrelated with the fixed covariates, and it can be significantly less computationally intensive than ICAR to estimate [@Broms_2014]. +Broms et al. [@Broms_2014] provide a nice overview of the concepts. +Several recent studies have shown that RSR provides better and faster inference than ICAR [e.g. @Broms_2014; @Clark_2019]. +At least two specialized R packages ([`stocc`](https://cran.r-project.org/web/packages/stocc/index.html) and [`Rcppocc`](https://github.com/AllanClark/Rcppocc)) are available for fitting ICAR and RSR occupancy models. + +It is now possible to fit RSRs in Stan with `ubms` as well, allowing those familiar with `unmarked`/`ubms` syntax and workflows to easily fit spatial models. +Currently only single-season occupancy models are supported, but support will soon be extended to other single-season model types. +In addition to fitting RSRs, `ubms` also provides tools for applying and visualizing fitted spatial models. +For now these features are available in a dev branch of `ubms` available on [Github](https://github.com/kenkellner/ubms/tree/RSR). + +# Example + +We will adapt the demo code and dataset available in R package `stocc` to fit a spatial occupancy model, which will also facilitate a comparison between `stocc` and `ubms` results. +The `stocc` package was an extremely valuable guide when implementing RSRs in `ubms`, and contains many other useful features related to spatial occupancy modeling. + +```{r,eval=FALSE} +library(ubms) +library(stocc) +``` + +```{r, echo=FALSE, message=FALSE, warning=FALSE} +suppressMessages(library(ubms)) +suppressMessages(library(stocc)) +``` + +## Format the input data + +The sample dataset that comes with `stocc` is in long format. +We need to convert it to the wide format that `unmarked` and `ubms` expect. + + +```{r} +# Load datasets +data(habData) # site covariate data +data(visitData) # observation cov data + +# Visits per site +nobs <- table(visitData$site, visitData$obs) + +# Dimensions of y matrix +J <- max(nobs) # max obs per site +M <- nrow(habData) # no. of sites + +# Create blank matrices to hold wide-format versions +y <- detCov1 <- detCov2 <- matrix(NA, nrow=M, ncol=J) + +# Iterate over sites +for (i in 1:M){ + # Data from site i + sub <- visitData[visitData$site==i,] + # Put this data into our wide matrices + subJ <- nrow(sub) + if(subJ == 0) next + y[i,1:subJ] <- sub$obs + detCov1[i,1:subJ] <- sub$detCov1 + detCov2[i,1:subJ] <- as.character(sub$detCov2) +} + +# Make site and obs cov data frames +site_cov <- habData +obs_cov <- data.frame(detCov1=as.vector(t(detCov1)), detCov2=as.vector(t(detCov2))) + +# Convert detCov2 back into a factor (as it was originally) +obs_cov$detCov2 <- factor(obs_cov$detCov2, levels=c("0","1","2","3")) +``` + +In order to fit a spatial model, we need coordinates included in our site covariates. +These should be in some projected coordinate system. +In this dataset, coordinates are contained in columns `x` and `y`: + +```{r} +head(site_cov) +``` + +Site and observation covariates are otherwise handled the same way as normal for `unmarked` and `ubms` models. +Here's a quick visualization of how `habCov2` varies across the study area. + +```{r,fig.width=6} +library(ggplot2) +#Which sites were sampled at least once? +sampled <- apply(y, 1, function(x) any(!is.na(x))) + +ggplot(data=site_cov, aes(x=x, y=y)) + + geom_tile(aes(fill=habCov2)) + + geom_point(data=site_cov[sampled,], size=0.5) + + scale_fill_gradientn(colors=heat.colors(10)) + + theme_bw(base_size=12) + + theme(panel.grid=element_blank()) +``` + +The black dots represent locations that were sampled at least once. +Note that not all grid cells with available site covariates were actually sampled; in fact most of them were not. +For typical models, `ubms` would just ignore a site that was present in the dataset but never sampled (i.e., `y` all `NA`) since these sites cannot contribute any information. +However for spatial models, generating predictive maps across areas that may not have been sampled is a typical goal. +To facilitate this, spatial models in `ubms` do not discard sites which were never sampled but that have available site covariates. +As you'd expect, you simply include these by including a row of all missing values in `y` and your observation covariates, as appropriate. +For example, in this dataset, sites 1 and 2 were sampled but site 3 was not: + +```{r} +y[1:3,] +round(detCov1[1:3,],3) +``` + +The final step in data formatting is to construct the `unmarkedFrame`. + +```{r, warning=FALSE} +umf <- unmarkedFrameOccu(y=y, siteCovs=site_cov, obsCovs=obs_cov) +``` + +## Choose RSR options + +Next, we need to decide on settings for the RSR. +Specifically, we need to decide the distance threshold below which two sites will be considered neighbors and thus potentially correlated with each other. +For example, if site A and site B are 0.5 distance units apart and we set the threshold at 1, A and B will be considered neighbors. +If A and B are 1.5 distance units apart, they will not be neighbors. +The distance units are defined by the units of the coordinates you provided. + +To visualize this, we can use the `RSR()` function, which we will use again later when fitting the model. +The RSR() function requires as input at least one coordinate vector (typically you will have 2, e.g. `x` and `y`) and a value for the `threshold` argument. +If we set the `plot_site` argument to an integer, `RSR` will return a map highlighting that site and its neighbors under the current threshold setting: + +```{r, fig.width=6} +with(site_cov, RSR(x, y, threshold=1, plot_site=300)) +``` + +With `threshold=1` and our coordinate grid, site 300 has just 4 neighbors. +Now consider `threshold=10`: + +```{r, fig.width=6} +with(site_cov, RSR(x, y, threshold=10, plot_site=300)) +``` + +The threshold choice depends on your modeling goals and study system, but it is probably best to start with sites having a relatively small number of neighbors. +We will set `threshold=1` for this example. + +A final argument available in `RSR()` is `moran_cut`. +This option controls the number of eigenvectors that are used when calculating the spatial random effect (possible values are 1 to the number of sites). +Generally, with fewer eigenvectors, the model will run faster and the result will be smoother. +Previous studies have shown that inference is not particularly sensitive to this value, and recommended it be set to 10% of the number of sites [@Broms_2014; @Clark_2019]. +This is the default in `ubms` so we will leave it alone. + +## Fit the model + +We are finally ready to fit our spatial model. +The formula syntax is identical to normal `ubms` and (and `unmarked`) models, except that we will include a call to `RSR()` with our desired options in our formula for `psi`. +We also include two fixed detection and two fixed occupancy covariates. + +```{r} +#Double formula: first part is for detection, second for occupancy +form <- ~detCov1 + detCov2 ~habCov1 + habCov2 + RSR(x, y, threshold=1) +``` + +We will then fit the model in Stan with `stan_occu()`, first setting our desired number of parallel cores to use to 3. +This will take Stan a few minutes with default MCMC settings. + +```{r, eval=FALSE} +options(mc.cores=3) +fit_ubms <- stan_occu(form, umf, chains=3) +``` + +```{r, message=FALSE, echo=FALSE, warning=FALSE} +if(file.exists("fit_ubms.Rds")){ + fit_ubms <- readRDS("fit_ubms.Rds") +} else { + options(mc.cores=3) + fit_ubms <- stan_occu(form, umf, chains=3, refresh=0) + saveRDS(fit_ubms, "fit_ubms.Rds") +} +``` + +You will likely get some warnings about effective sample size. +The solution is to run the model for more iterations, but this is good enough for an example. + +## Examine results + +Look at the parameter estimates: + +```{r} +fit_ubms +``` + +The precision term (tau, $\tau$) associated with the spatial random effect is the final row of the occupancy model estimates. +Smaller values of $\tau$ imply greater variability in the spatial random effect. + +There also appears to be a strong effect of `habCov2` on occupancy. +We can quickly visualize marginal effects of fixed covariates with `plot_marginal`. + +```{r, fig.width=5} +plot_marginal(fit_ubms, "state") +``` + +You can extract predicted occupancy values for each site (including sites never actually sampled) using the `predict` function: + +```{r, message=FALSE} +# "state" to get state/occupancy process, as opposed to "det" +psi <- predict(fit_ubms, "state") +head(psi) +``` + +## Plotting + +It is much more interesting to see how predicted occupancy varies across the study area. +`ubms` has a built-in function for generating such a map, `plot_spatial`. + +```{r,message=FALSE, fig.width=6} +plot_spatial(fit_ubms) +``` + +The plot also shows the location of sampled sites and whether the species was ever detected at each one. + +You can also look at the spatial distribution of random effects $\eta$: + +```{r,message=FALSE, fig.width=6} +plot_spatial(fit_ubms, "eta") +``` + +## Model selection + +It is important to confirm that the spatial random effect is actually improving the predictive power of the model; if it isn't, we ought to remove it. +To check this, first fit a similar model without the RSR: +```{r, eval=FALSE} +fit_nonspatial <- stan_occu(~detCov1 + detCov2 ~habCov1 + habCov2, umf, chains=3) +``` + +```{r, echo=FALSE, message=FALSE, warning=FALSE} +if(file.exists("fit_nonspatial.Rds")){ + fit_nonspatial <- readRDS("fit_nonspatial.Rds") +} else { + options(mc.cores=3) + fit_nonspatial <- stan_occu(~detCov1 + detCov2 ~habCov1 + habCov2, umf, + chains=3, refresh=0) + saveRDS(fit_nonspatial, "fit_nonspatial.Rds") +} +``` + +Combine the models to compare into a `fitList` (note: they must have the same number of chains and MCMC iterations): + +```{r} +fl <- fitList(fit_ubms, fit_nonspatial) +``` + +Then, compare the spatial and non-spatial models with leave-one-out cross-validation (LOO)[@Vehtari_2017] using the `modSel()` function: + +```{r, warning=FALSE} +round(modSel(fl),2) +``` + +A measure of predictive accuracy (`elpd`) is shown for each model, and the model with the highest predictive accuracy is listed first. +The difference in `elpd` between the top and model and each subsequent model (`elpd_diff`), along with a measure of the uncertainty of that difference (`se_diff`), are also shown. +If the magnitude of `elpd_diff` is large relative to the associated uncertainty, we can conclude that there is a difference in predictive power between the two models. +In this case `elpd_diff` between the spatial (`fit_ubms`) and non-spatial (`fit_nonspatial`) models is several times the size of `se_diff`, so including the spatial random effect appears to improve model predictive accuracy. + +# Compare with stocc + +The implementation of spatial models in `ubms` owes much to the `stocc` package. +The two packages give essentially identical results. +To show this, we will fit the same model in `stocc` and compare the occupancy predictions with `ubms`. +The code for fitting the model in `stocc` is taken directly from the package demo. + +```{r, eval=FALSE} +names <- list( + visit=list(site="site",obs="obs"), + site=list(site="site", coords=c("x","y")) + ) + +fit_stocc <- spatial.occupancy( + detection.model = ~ detCov1 + detCov2, + occupancy.model = ~ habCov1 + habCov2, + spatial.model = list(model="rsr", threshold=1, moran.cut=0.1*nrow(habData)), + so.data = make.so.data(visitData, habData, names), + prior = list(a.tau=0.5, b.tau=0.00005, Q.b=0.1, Q.g=0.1), + control = list(burnin=1000/5, iter=4000/5, thin=5) + ) +``` + +```{r, echo=FALSE,message=FALSE, warning=FALSE} +names <- list( + visit=list(site="site",obs="obs"), + site=list(site="site", coords=c("x","y")) + ) + +if(file.exists("fit_stocc.Rds")){ + fit_stocc <- readRDS("fit_stocc.Rds") +} else { + # Takes a few minutes + fit_stocc <- spatial.occupancy( + detection.model = ~ detCov1 + detCov2, + occupancy.model = ~ habCov1 + habCov2, + spatial.model = list(model="rsr", threshold=1, moran.cut=0.1*nrow(habData)), + so.data = make.so.data(visitData, habData, names), + prior = list(a.tau=0.5, b.tau=0.00005, Q.b=0.1, Q.g=0.1), + control = list(burnin=1000/5, iter=4000/5, thin=5) + ) + saveRDS(fit_stocc, "fit_stocc.Rds") +} +``` + +Here's a similar plot of predicted occupancy across the study area: + +```{r, fig.width=6} +# Predicted occupancy +psi_stocc <- fit_stocc$occupancy.df$psi.est + +# Sites that were sampled at least once +sampled <- apply(getY(umf), 1, function(x) any(!is.na(x))) +plot_data <- cbind(habData, psi=psi_stocc) + +ggplot(data=plot_data, aes(x=x,y=y)) + + geom_point(aes(col=psi), size=5, pch=15) + + scale_color_gradientn(colors=terrain.colors(10)) + + geom_point(data=plot_data[sampled,],aes(x=x,y=y), size=0.5) + + labs(col="psi") + + theme_bw(base_size=12) + + theme(panel.grid=element_blank(), + strip.background=element_rect("white")) +``` + +As you can see, the predicted occupancy maps for `ubms` and `stocc` are almost identical. +Note, however, that coefficient estimates may not be identical because `stocc` uses a probit link function, while `ubms` uses a logit link. + +# References + +
diff --git a/src/ubms-vignette.Rmd b/src/ubms-vignette.Rmd index 8d328e3..870a6c7 100644 --- a/src/ubms-vignette.Rmd +++ b/src/ubms-vignette.Rmd @@ -22,7 +22,7 @@ knitr::opts_chunk$set(message=FALSE, warning=FALSE) ## What is ubms? -`ubms` is an R package for fitting models of wildlife occurrence and abundance to with datasets where animals are not individually marked. +`ubms` is an R package for fitting models of wildlife occurrence and abundance with datasets where animals are not individually marked. It provides a nearly identical interface to the popular `unmarked` package [@Fiske_2011]. Instead of using maximum likelihood to fit models (as with `unmarked`), models are fit in a Bayesian framework using [Stan](https://mc-stan.org/) [@Carpenter_2017]. @@ -150,7 +150,7 @@ options(mc.cores=3) ```{r, echo=FALSE} library(ubms) -(fit_stan <- stan_occu(~1~1, data=umf, chains=3, iter=500, refresh=0)) +(fit_stan <- stan_occu(~1~1, data=umf, chains=3, iter=1000, refresh=0)) ``` ### Compare results @@ -232,13 +232,13 @@ fit_global <- stan_occu(~scale(date)~scale(forest)+scale(ele), data=umf, ```{r, echo=FALSE} fit_null <- fit_stan fit_sc <- stan_occu(~1~scale(forest)+scale(ele), data=umf, - chains=3, iter=500, refresh=0) -fit_oc <- stan_occu(~scale(date)~1, data=umf, chains=3, iter=500, refresh=0) + chains=3, iter=1000, refresh=0) +fit_oc <- stan_occu(~scale(date)~1, data=umf, chains=3, iter=1000, refresh=0) ``` ```{r, warning=TRUE, echo=FALSE} fit_global <- stan_occu(~scale(date)~scale(forest)+scale(ele), data=umf, - chains=3, iter=500, refresh=0) + chains=3, iter=1000, refresh=0) ``` The `fit_global` model gave us some warnings about the effective sample size (`n_eff`) along with a suggested solution. -- cgit v1.2.3