| Title: | Bayesian Models for Data from Unmarked Animals using 'Stan' |
|---|---|
| Description: | Fit Bayesian hierarchical models of animal abundance and occurrence via the 'rstan' package, the R interface to the 'Stan' C++ library. Supported models include single-season occupancy, dynamic occupancy, and N-mixture abundance models. Covariates on model parameters are specified using a formula-based interface similar to package 'unmarked', while also allowing for estimation of random slope and intercept terms. References: Carpenter et al. (2017) <doi:10.18637/jss.v076.i01>; Fiske and Chandler (2011) <doi:10.18637/jss.v043.i10>. |
| Authors: | Ken Kellner [cre, aut] |
| Maintainer: | Ken Kellner <[email protected]> |
| License: | GPL (>=3) |
| Version: | 1.2.7 |
| Built: | 2024-10-30 06:46:55 UTC |
| Source: | https://github.com/biodiverse/ubms |
Extract a Submodel from a ubmsFit Model
## S4 method for signature 'ubmsFit,character,missing,missing' x[i]## S4 method for signature 'ubmsFit,character,missing,missing' x[i]
x |
A |
i |
The name of a submodel to extract |
An object of class ubmsSubmodel.
Extract a ubmsSubmodel From a ubmsSubmodelList Object
## S4 method for signature 'ubmsSubmodelList,character,missing,missing' x[i]## S4 method for signature 'ubmsSubmodelList,character,missing,missing' x[i]
x |
Object of class |
i |
The name of a submodel |
An object of class ubmsSubmodel.
Extract Coefficient Values From a ubmsFit Model
## S4 method for signature 'ubmsFit' coef(object, ...)## S4 method for signature 'ubmsFit' coef(object, ...)
object |
A |
... |
Currently ignored |
A vector of coefficient values for all submodels.
Extracts the pointwise log-likelihood matrix or array from a model.
This is useful as an input for functions in the loo package.
If called on a ubmsFit object, the log-likelihood matrix or array
is calculated using the posterior distributions of the parameters; the
log-likelihood values are not actually saved inside the model object.
If called on a stanfit object, loo::extract_log_lik is used.
In this case, the log-likelihood must be saved as one of the output
parameters in Stan.
extract_log_lik(object, parameter_name = "log_lik", merge_chains = TRUE)extract_log_lik(object, parameter_name = "log_lik", merge_chains = TRUE)
object |
A |
parameter_name |
The name of the log-likelihood parameter in the
Stan output; ignored when |
merge_chains |
If |
A matrix (samples x sites) or array (samples x chains x sites)
Extract samples from a ubmsFit model
## S4 method for signature 'ubmsFit' extract(object, pars, permuted = TRUE, inc_warmup = FALSE, include = TRUE)## S4 method for signature 'ubmsFit' extract(object, pars, permuted = TRUE, inc_warmup = FALSE, include = TRUE)
object |
A |
pars |
An optional character vector providing parameter names of interest. If not specified, all parameters are used |
permuted |
Logical. If |
inc_warmup |
Logical. If |
include |
Logical. If |
If permuted=TRUE, a list; if permuted=FALSE,
an array.
Create a list of ubmsFit models
## S4 method for signature 'ubmsFit' fitList(...)## S4 method for signature 'ubmsFit' fitList(...)
... |
|
An object of class ubmsFitList containing the list of models
Extract fitted values for a given submodel from a ubmsFit object.
Fitted values are calculated separately for each submodel
using the posterior predictive distribution of the latent state z,
following Wright et al. (2019).
## S4 method for signature 'ubmsFit' fitted(object, submodel, draws = NULL, ...)## S4 method for signature 'ubmsFit' fitted(object, submodel, draws = NULL, ...)
object |
A fitted model of class |
submodel |
Submodel to get fitted values for, options are |
draws |
An integer indicating the number of draws to return. The default and maximum number of draws is the size of the posterior sample. |
... |
Currently ignored |
A matrix of fitted values with dimension draws by
observations. Note that calculation of fitted values
for the detection submodel is conditional on , so fitted values
for an observation in a posterior draw where are assigned
value NA (Wright et al. 2019).
Wright, W. J., Irvine, K. M., & Higgs, M. D. (2019). Identifying occupancy model inadequacies: can residuals separately assess detection and presence? Ecology, 100(6), e02703.
Get warmup and sampling time from a ubmsFit object
## S4 method for signature 'ubmsFit' get_elapsed_time(object, ...)## S4 method for signature 'ubmsFit' get_elapsed_time(object, ...)
object |
A |
... |
Arguments passed to |
A matrix with one row per chain and two columns, containing the warmup time and sampling time, respectively, in seconds
Get the Stan code used to run a model as a character string
## S4 method for signature 'ubmsFit' get_stancode(object, ...)## S4 method for signature 'ubmsFit' get_stancode(object, ...)
object |
A |
... |
Arguments passed to |
Pass the result of get_stancode to cat to get the
code in a more readable format. Note that the output in most cases
is Stan code that can be used to fit several types of models, and not
all Stan code will be used in all models.
A character string with the model code
Extract y, the Response Variable, From a ubmsFit Model
## S4 method for signature 'ubmsFit' getY(object)## S4 method for signature 'ubmsFit' getY(object)
object |
A |
A matrix containing the response variable y.
Goodness-of-fit tests for ubmsFit models using posterior predictive
checks
gof(object, draws = NULL, ...) ## S4 method for signature 'ubmsFitOccu' gof(object, draws = NULL, quiet = FALSE, ...) ## S4 method for signature 'ubmsFitAbun' gof(object, draws = NULL, quiet = FALSE, ...)gof(object, draws = NULL, ...) ## S4 method for signature 'ubmsFitOccu' gof(object, draws = NULL, quiet = FALSE, ...) ## S4 method for signature 'ubmsFitAbun' gof(object, draws = NULL, quiet = FALSE, ...)
object |
A fitted model of class |
draws |
Number of draws from the posterior to use in the check |
... |
Currently ignored |
quiet |
If |
An object of class ubmsGOF containing statistics calculated
from the posterior predictive distribution.
gof(ubmsFitOccu): Applies the MacKenzie-Bailey chi-square goodness of fit test for
ocupancy models (MacKenzie and Bailey 2004).
gof(ubmsFitAbun): A goodness-of-fit test for N-mixture type models based on Pearson's chi-square.
MacKenzie, D. I., & Bailey, L. L. (2004). Assessing the fit of site-occupancy models. Journal of Agricultural, Biological, and Environmental Statistics, 9(3), 300-318.
Randomly partition data into K subsets of equal size (by site). Re-fit the model
K times, each time leaving out one of the subsets. Calculate the log-likelihood
for each of the sites that was left out. This function is an alternative
to loo (leave-one-out cross validation).
## S4 method for signature 'ubmsFit' kfold(x, K = 10, folds = NULL, quiet = FALSE, ...)## S4 method for signature 'ubmsFit' kfold(x, K = 10, folds = NULL, quiet = FALSE, ...)
x |
A |
K |
Number of folds into which the data will be partitioned |
folds |
An optional vector with length equal to the number of sites in the data and containing integers from 1 to K, to manually assign sites to folds. You should use this if you plan to compare multiple models, since the folds for each model should be identical. You can use |
quiet |
If |
... |
Currently ignored |
An object of class elpd_generic that is compatible with loo::loo_compare
Leave-one-out Cross Validation
## S4 method for signature 'ubmsFit' loo(x, ..., cores = getOption("mc.cores", 1))## S4 method for signature 'ubmsFit' loo(x, ..., cores = getOption("mc.cores", 1))
x |
A |
... |
Currently ignored |
cores |
Number of cores to use for calculation |
A named list of class loo with estimates of
the expected log predictive density and other parameters used
for model comparison. See ?loo::loo for more information.
Construct a model selection table from a ubmsFitList
## S4 method for signature 'ubmsFitList' modSel(object, ...)## S4 method for signature 'ubmsFitList' modSel(object, ...)
object |
An object of class |
... |
Currently ignored |
A data.frame of model fit information with one row per
model in the input ubmsFitList. Models are ranked in descending
order by expected log pointwise predictive density (elpd).
Get Parameter Names From a ubmsFit Model
## S4 method for signature 'ubmsFit' names(x)## S4 method for signature 'ubmsFit' names(x)
x |
A |
A character vector of parameter names.
Get Names of Models in a ubmsFitList
## S4 method for signature 'ubmsFitList' names(x)## S4 method for signature 'ubmsFitList' names(x)
x |
A |
A character vector of model names.
Get number of Posterior Samples Stored in a ubmsFit Model
## S4 method for signature 'ubmsFit' nsamples(object, ...)## S4 method for signature 'ubmsFit' nsamples(object, ...)
object |
A |
... |
Currently ignored |
An integer representing the number of posterior samples
Generates marginal fixed effects plots of one or more covariates from a
ubmsFit submodel. For each plot, the focal covariate is allowed to
vary across its range (or possible discrete values, for a factor), while
the other covariates are held at their medians or reference levels.
Random effects are ignored.
## S4 method for signature 'ubmsFit' plot_effects( object, submodel, covariate = NULL, level = 0.95, draws = 1000, smooth = NULL, ... ) ## S4 method for signature 'ubmsFit' plot_marginal( object, submodel, covariate = NULL, level = 0.95, draws = 1000, smooth = NULL, ... )## S4 method for signature 'ubmsFit' plot_effects( object, submodel, covariate = NULL, level = 0.95, draws = 1000, smooth = NULL, ... ) ## S4 method for signature 'ubmsFit' plot_marginal( object, submodel, covariate = NULL, level = 0.95, draws = 1000, smooth = NULL, ... )
object |
A fitted model of class |
submodel |
Submodel to get plots for, for example |
covariate |
Plot a specific covariate; provide the name as a string |
level |
Probability mass to include in the uncertainty interval |
draws |
Number of draws from the posterior to use when generating the
plots. If fewer than |
smooth |
Smoother span (f) value used for LOWESS smoothing of the
upper and lower credible interval bounds for a continuous covariate.
No smoothing is done if |
... |
Currently ignored |
A ggplot if a single covariate is plotted, or an object
of class grob if there are multiple covariates/panels
Plot posterior distributions for selected parameters. Posteriors can be represented as density plots or histograms.
## S4 method for signature 'ubmsFit' plot_posteriors(object, pars = NULL, density = FALSE, ...)## S4 method for signature 'ubmsFit' plot_posteriors(object, pars = NULL, density = FALSE, ...)
object |
A fitted model of class |
pars |
A character vector of parameter names to include in the plot
Look at |
density |
If |
... |
Arguments passed to |
A ggplot
Plot residuals for a submodel from a ubmsFit object, for multiple
posterior draws. By default, residuals are plotted against fitted values.
When the submodel has a binomial response (e.g., detection models), regular
residual plots are not typically informative. Instead, the residuals and
fitted values are divided into bins based on fitted value and the averages
are plotted. For a count response (e.g., Poisson), Pearson residuals are
calculated. To plot residuals against values of a particular covariate instead
of the fitted values, supply the name of the covariate (as a string) to the
covariate argument.
## S4 method for signature 'ubmsFit' plot_residuals( object, submodel, covariate = NULL, draws = 9, nbins = NULL, ... )## S4 method for signature 'ubmsFit' plot_residuals( object, submodel, covariate = NULL, draws = 9, nbins = NULL, ... )
object |
A fitted model of class |
submodel |
Submodel to plot residuals for, for example |
covariate |
If specified, plot residuals against values of a covariate.
Covariate name should be provided as a string. If |
draws |
An integer indicating the number of posterior draws to use. Separate plots are generated for each draw, so this number should be relatively small. The default and maximum number of draws is the size of the posterior sample. |
nbins |
For submodels with a binomial response, manually set the number of bins to use |
... |
Currently ignored |
A ggplot of residuals vs. fitted values or covariate values,
with one panel per posterior draw. For binned residual plots, the shaded area
represents plus/minus two standard deviations around the mean residual.
If the model is true, we would expect about 95
fall within this area.
Plot A Map of the State Parameter Based on a Spatial ubms Model
plot_spatial(object, param = c("state", "eta"), sites = TRUE, cell_size = NULL)plot_spatial(object, param = c("state", "eta"), sites = TRUE, cell_size = NULL)
object |
A |
param |
The parameter to plot, either |
sites |
If |
cell_size |
The size of each side of the (square) cells drawn in the map,
in the same units as the coordinates. If |
Plot Residuals For All Submodels in a ubmsFit Model
## S4 method for signature 'ubmsFit,ANY' plot(x, y, ...)## S4 method for signature 'ubmsFit,ANY' plot(x, y, ...)
x |
A |
y |
Currently ignored |
... |
Currently ignored |
A plot object of class gtable with one panel per submodel.
Extract posterior draws of the linear predictor for a ubmsFit
submodel, possibly transformed by the inverse-link function.
## S4 method for signature 'ubmsFit' posterior_linpred( object, transform = FALSE, submodel, newdata = NULL, draws = NULL, re.form = NULL, ... )## S4 method for signature 'ubmsFit' posterior_linpred( object, transform = FALSE, submodel, newdata = NULL, draws = NULL, re.form = NULL, ... )
object |
A fitted model of class |
transform |
Should the linear predictor be transformed using the inverse link function? |
submodel |
The name of the submodel, as a character string, for which to calculate the linear predictor |
newdata |
Optional data frame of newdata to use when calculating the linear predictor. If not provided, the model matrix is used. |
draws |
An integer indicating the number of draws to return. The default and maximum number of draws is the size of the posterior sample. |
re.form |
If |
... |
Currently ignored |
A matrix of simulations from the posterior predictive distribution
of the linear predictor. The dimensions are draws by number of
linear predictor values (e.g., number of sites or number of observations).
Draw from the posterior predictive distribution after fitting a model.
You can draw from the posterior of the observed outcome or
the latent unobserved state .
## S4 method for signature 'ubmsFit' posterior_predict( object, param = c("y", "z"), draws = NULL, re.form = NULL, ... )## S4 method for signature 'ubmsFit' posterior_predict( object, param = c("y", "z"), draws = NULL, re.form = NULL, ... )
object |
A fitted model of class |
param |
Either |
draws |
An integer indicating the number of draws to return. The default and maximum number of draws is the size of the posterior sample. |
re.form |
If |
... |
Currently ignored |
A matrix of simulations from the posterior predictive distribution.
If param = "z", the dimensions are draws by number of sites
(or sites x primary periods in site-major order for dynamic models). If
param = "y", the dimensions are draws by sites x observations
(or sites x primary periods x observations for dynamic models).
This method generates predicted parameter values for the original dataset or a new dataset using the posterior distribution. Standard deviation and a customizable uncertainty interval are also calculated.
## S4 method for signature 'ubmsFit' predict( object, submodel, newdata = NULL, transform = TRUE, re.form = NULL, level = 0.95, ... )## S4 method for signature 'ubmsFit' predict( object, submodel, newdata = NULL, transform = TRUE, re.form = NULL, level = 0.95, ... )
object |
A fitted model of class |
submodel |
Submodel to predict from, for example |
newdata |
Optional data frame, SpatRaster, or RasterStack of covariates to generate predictions from. If not provided (the default), predictions are generated from the original data |
transform |
If |
re.form |
If |
level |
Probability mass to include in the uncertainty interval |
... |
Currently ignored |
If newdata was a data frame: A data frame with one row per
prediction and four columns:
1) Predicted point estimates (posterior means),
2) Standard deviation of the posterior,
3-4) Lower and upper bounds of the specified uncertainty interval
For parameters with more than one dimension, the rows are in site-major order, or site-year-observation for dynamic models.
If newdata was a SpatRaster/RasterStack, returns a SpatRaster/RasterStack
with four layers corresponding to the four columns above with the same projection
as the original SpatRaster/RasterStack.
posterior_linpred, posterior_interval
Specify prior distributions and associated parameters for use in
ubms models.
normal(location = 0, scale = 2.5, autoscale = TRUE) uniform(lower = -5, upper = 5) student_t(df = 1, location = 0, scale = 2.5, autoscale = TRUE) logistic(location = 0, scale = 1) cauchy(location = 0, scale = 2.5, autoscale = TRUE) gamma(shape = 1, rate = 1) laplace(location = 0, scale = 2.5, autoscale = TRUE)normal(location = 0, scale = 2.5, autoscale = TRUE) uniform(lower = -5, upper = 5) student_t(df = 1, location = 0, scale = 2.5, autoscale = TRUE) logistic(location = 0, scale = 1) cauchy(location = 0, scale = 2.5, autoscale = TRUE) gamma(shape = 1, rate = 1) laplace(location = 0, scale = 2.5, autoscale = TRUE)
location |
The mean of the distribution. If setting the priors for regression coefficients, this can be a single value, or multiple values, one per coefficient |
scale |
The standard deviation of the distribution. If setting the priors for regression coefficients, this can be a single value, or multiple values, one per coefficient |
autoscale |
If |
lower |
The lower bound for the uniform distribution |
upper |
The upper bound for the uniform distribution |
df |
The number of degrees of freedom for the Student-t distribution |
shape |
The gamma distribution shape parameter |
rate |
The gamma distribution rate parameter (1/scale) |
A list containing prior settings used internally by ubms.
normal()normal()
Generate posterior draws of occupancy probability for all sites and primary periods, i.e. the projected trajectory (Weir et al. 2009).
projected(object, ...) ## S4 method for signature 'ubmsFitColext' projected(object, draws = NULL, re.form = NULL, ...)projected(object, ...) ## S4 method for signature 'ubmsFitColext' projected(object, draws = NULL, re.form = NULL, ...)
object |
A fitted dynamic occupancy model of class inheriting |
... |
Currently ignored |
draws |
Number of draws from the posterior to use in the check |
re.form |
If |
A matrix of occupancy values from the posterior predictive distribution.
The dimensions are draws by number of sites x primary periods
in site-major order.
Weir LA, Fiske IJ, Royle J. 2009. Trends in Anuran Occupancy from Northeastern States of the North American Amphibian Monitoring Program. Herpetological Conservation and Biology. 4: 389-402.
Extract random effects from a ubmsFit model. Note that this function
works like ranef for merMod objects from lme4, not like
ranef for unmarkedFit objects. To get functionality similar
to that of unmarkedFit, use posterior_predict.
## S4 method for signature 'ubmsFit' ranef(object, submodel, summary = FALSE, add_mean = TRUE, ...)## S4 method for signature 'ubmsFit' ranef(object, submodel, summary = FALSE, add_mean = TRUE, ...)
object |
A fitted model of class |
submodel |
The name of the submodel, as a character string, for which to generate the random effects |
summary |
If |
add_mean |
If |
... |
Currently ignored |
Note: by default this function adds the overall intercept or slope
to the (mean-0) random effect to get the complete random intercept or slope.
In this way the output is more like the output of lme4::coef
and not lme4::ranef. You can turn this off and return just the
mean-0 random effect by setting argument add_mean = FALSE.
If you run ranef on a submodel with a spatial random effect,
the function will return estimates of parameter eta.
If summary=FALSE, a list of random effect values; if
TRUE, a data frame with columns for random effect mean, SD, and
95
Extract residuals for a given submodel from a ubmsFit object.
Residuals are calculated separately for each submodel
using the posterior predictive distribution of the latent state z,
following Wright et al. (2019).
## S4 method for signature 'ubmsFit' residuals(object, submodel, draws = NULL, ...)## S4 method for signature 'ubmsFit' residuals(object, submodel, draws = NULL, ...)
object |
A fitted model of class |
submodel |
Submodel to get residuals for, for example |
draws |
An integer indicating the number of draws to return. The default and maximum number of draws is the size of the posterior sample. |
... |
Currently ignored |
A matrix of residual values with dimension draws by
observations. Note that calculation of residuals
for the detection submodel is conditional on , so residuals
for an observation in a posterior draw where are assigned
value NA (Wright et al. 2019).
Wright, W. J., Irvine, K. M., & Higgs, M. D. (2019). Identifying occupancy model inadequacies: can residuals separately assess detection and presence? Ecology, 100(6), e02703.
A call to RSR in the formula for a state process adds an autocorrelated
spatial random effect to the model in the form of a Restricted Spatial
Regression (RSR). For examples of RSRs applied to ecological data, see
Johnson et al. (2013) and Broms et al. (2014).
The function can also be used outside a formula to visualize the spatial
relationships between sites in your data and choose an appropriate
distance threshold below which two sites will be considered neighbors, and
thus potentially correlated, for the RSR model. For more details, see the example vignette:
vignette("spatial-models", package = "ubms")
RSR(x, y = NULL, threshold, moran_cut = NULL, plot_site = NULL)RSR(x, y = NULL, threshold, moran_cut = NULL, plot_site = NULL)
x |
A vector of coordinates (should be projected) |
y |
An (optional) second vector of coordinates |
threshold |
The distance cutoff below which two sites will be considered neighbors. Should be the same units as the coordinates. |
moran_cut |
The number of eigenvectors to use in the RSR. The possible range of values is between 1 and the number of sites. Smaller numbers will result in faster runtime and smoother map output, and vice-versa. If not provided (the default), the number of eigenvectors will be set as 10% of the number of sites which is usually appropriate. |
plot_site |
If a site number (as an integer) is supplied, the function returns a plot showing that site and its neighbors under the current settings. Useful for deciding what to set your threshold at. |
Either a list of spatial matrices used for the RSR (only useful
internally to ubms), or if plot_site is an integer, a ggplot object.
Broms KM, Johnson DS, Altwegg R, Conquest LL. 2014. Spatial occupancy models applied to atlas data show Southern Ground Hornbills strongly depend on protected areas. Ecological Applications 24: 363-374.
Johnson DS, Conn PB, Hooten MB, Ray JC, Pond BA. 2013. Spatial occupancy models for large data sets. Ecology 94: 801-808.
# Generate some coordinates x <- runif(100,0,10) y <- runif(100,0,10) # Show neighbors of site 10 when threshold is 3 units RSR(x, y, threshold=3, plot_site=10)# Generate some coordinates x <- runif(100,0,10) y <- runif(100,0,10) # Show neighbors of site 10 when threshold is 3 units RSR(x, y, threshold=3, plot_site=10)
This function fits the dynamic occupancy model of MacKenzie et al. (2003).
stan_colext( psiformula = ~1, gammaformula = ~1, epsilonformula = ~1, pformula = ~1, data, prior_intercept_psi = logistic(0, 1), prior_coef_psi = logistic(0, 1), prior_intercept_gamma = logistic(0, 1), prior_coef_gamma = logistic(0, 1), prior_intercept_eps = logistic(0, 1), prior_coef_eps = logistic(0, 1), prior_intercept_det = logistic(0, 1), prior_coef_det = logistic(0, 1), prior_sigma = gamma(1, 1), log_lik = TRUE, ... )stan_colext( psiformula = ~1, gammaformula = ~1, epsilonformula = ~1, pformula = ~1, data, prior_intercept_psi = logistic(0, 1), prior_coef_psi = logistic(0, 1), prior_intercept_gamma = logistic(0, 1), prior_coef_gamma = logistic(0, 1), prior_intercept_eps = logistic(0, 1), prior_coef_eps = logistic(0, 1), prior_intercept_det = logistic(0, 1), prior_coef_det = logistic(0, 1), prior_sigma = gamma(1, 1), log_lik = TRUE, ... )
psiformula |
Right-hand sided formula for the initial probability of occupancy at each site |
gammaformula |
Right-hand sided formula for colonization probability |
epsilonformula |
Right-hand sided formula for extinction probability |
pformula |
Right-hand sided formula for detection probability |
data |
A |
prior_intercept_psi |
Prior distribution for the intercept of the
psi (initial occupancy probability) model; see |
prior_coef_psi |
Prior distribution for the regression coefficients of the psi model |
prior_intercept_gamma |
Prior distribution on intercept for colonization probability |
prior_coef_gamma |
Prior distribution on regression coefficients for colonization probability |
prior_intercept_eps |
Prior distribution on intercept for extinction probability |
prior_coef_eps |
Prior distribution on regression coefficients for extinction probability |
prior_intercept_det |
Prior distribution for the intercept of the detection probability model |
prior_coef_det |
Prior distribution for the regression coefficients of the detection model |
prior_sigma |
Prior distribution on random effect standard deviations |
log_lik |
If |
... |
Arguments passed to the |
ubmsFitColext object describing the model fit.
MacKenzie DI, Nicholas JD, Hines JE, Knutson MG, Franklin AB. 2003. Ecology 84: 2200-2207.
data(frogs) umf <- formatMult(masspcru) umf@y[umf@y > 1] <- 1 #convert counts to presence/absence umf <- umf[1:100,] #Use only 100 sites fit_frog <- stan_colext(~1, ~1, ~1, ~1, umf, chains=3, iter=300)data(frogs) umf <- formatMult(masspcru) umf@y[umf@y > 1] <- 1 #convert counts to presence/absence umf <- umf[1:100,] #Use only 100 sites fit_frog <- stan_colext(~1, ~1, ~1, ~1, umf, chains=3, iter=300)
This function fits the hierarchical distance sampling model of Royle et al. (2004) to line or point transect data recorded in discerete distance intervals.
stan_distsamp( formula, data, keyfun = c("halfnorm", "exp", "hazard"), output = c("density", "abund"), unitsOut = c("ha", "kmsq"), prior_intercept_state = normal(0, 5), prior_coef_state = normal(0, 2.5), prior_intercept_det = normal(0, 5), prior_coef_det = normal(0, 2.5), prior_intercept_scale = normal(0, 2.5), prior_sigma = gamma(1, 1), ... )stan_distsamp( formula, data, keyfun = c("halfnorm", "exp", "hazard"), output = c("density", "abund"), unitsOut = c("ha", "kmsq"), prior_intercept_state = normal(0, 5), prior_coef_state = normal(0, 2.5), prior_intercept_det = normal(0, 5), prior_coef_det = normal(0, 2.5), prior_intercept_scale = normal(0, 2.5), prior_sigma = gamma(1, 1), ... )
formula |
Double right-hand side formula describing covariates of detection and occupancy in that order |
data |
A |
keyfun |
One of the following detection functions:
|
output |
Model either density |
unitsOut |
Units of density. Either |
prior_intercept_state |
Prior distribution for the intercept of the
state (abundance) model; see |
prior_coef_state |
Prior distribution for the regression coefficients of the state model |
prior_intercept_det |
Prior distribution for the intercept of the detection probability model |
prior_coef_det |
Prior distribution for the regression coefficients of the detection model |
prior_intercept_scale |
Prior distribution for the intercept of the scale parameter (i.e., log(scale)) for Hazard-rate models |
prior_sigma |
Prior distribution on random effect standard deviations |
... |
Arguments passed to the |
ubmsFitDistsamp object describing the model fit.
Use of the hazard-rate key function ("hazard")
typically requires a large sample size in order to get good parameter
estimates. If you have a relatively small number of points/transects (<100),
you should be cautious with the resulting models. Check your results against
estimates from unmarked, which doesn't require as much data to get
good estimates of the hazard-rate shape and scale parameters.
Values of 'dist.breaks' in the 'unmarkedFrameDS' should be as small as possible (<10) to facilitate convergence. Consider converting 'unitsIn' from meters to kilometers, for example. See example below.
Royle, J. A., Dawson, D. K., & Bates, S. (2004). Modeling abundance effects in distance sampling. Ecology 85: 1591-1597.
data(issj) #Note use of km instead of m for distance breaks jayUMF <- unmarkedFrameDS(y=as.matrix(issj[,1:3]), siteCovs=issj[,c("elevation","forest")], dist.breaks=c(0,0.1,0.2,0.3), unitsIn="km", survey="point") fm_jay <- stan_distsamp(~1~scale(elevation), jayUMF, chains=3, iter=300)data(issj) #Note use of km instead of m for distance breaks jayUMF <- unmarkedFrameDS(y=as.matrix(issj[,1:3]), siteCovs=issj[,c("elevation","forest")], dist.breaks=c(0,0.1,0.2,0.3), unitsIn="km", survey="point") fm_jay <- stan_distsamp(~1~scale(elevation), jayUMF, chains=3, iter=300)
This function fits the multinomial-Poisson mixture model, useful for data collected via survey methods such as removal or double observer sampling.
stan_multinomPois( formula, data, prior_intercept_state = normal(0, 5), prior_coef_state = normal(0, 2.5), prior_intercept_det = logistic(0, 1), prior_coef_det = logistic(0, 1), prior_sigma = gamma(1, 1), log_lik = TRUE, ... )stan_multinomPois( formula, data, prior_intercept_state = normal(0, 5), prior_coef_state = normal(0, 2.5), prior_intercept_det = logistic(0, 1), prior_coef_det = logistic(0, 1), prior_sigma = gamma(1, 1), log_lik = TRUE, ... )
formula |
Double right-hand side formula describing covariates of detection and abundance in that order |
data |
A |
prior_intercept_state |
Prior distribution for the intercept of the
state (abundance) model; see |
prior_coef_state |
Prior distribution for the regression coefficients of the state model |
prior_intercept_det |
Prior distribution for the intercept of the detection probability model |
prior_coef_det |
Prior distribution for the regression coefficients of the detection model |
prior_sigma |
Prior distribution on random effect standard deviations |
log_lik |
If |
... |
Arguments passed to the |
ubmsFitMultinomPois object describing the model fit.
multinomPois, unmarkedFrameMPois
data(ovendata) ovenFrame <- unmarkedFrameMPois(ovendata.list$data, siteCovs=ovendata.list$covariates, type="removal") oven_fit <- stan_multinomPois(~1~scale(ufc), ovenFrame, chains=3, iter=300)data(ovendata) ovenFrame <- unmarkedFrameMPois(ovendata.list$data, siteCovs=ovendata.list$covariates, type="removal") oven_fit <- stan_multinomPois(~1~scale(ufc), ovenFrame, chains=3, iter=300)
This function fits the single season occupancy model of MacKenzie et al. (2002).
stan_occu( formula, data, prior_intercept_state = logistic(0, 1), prior_coef_state = logistic(0, 1), prior_intercept_det = logistic(0, 1), prior_coef_det = logistic(0, 1), prior_sigma = gamma(1, 1), log_lik = TRUE, ... )stan_occu( formula, data, prior_intercept_state = logistic(0, 1), prior_coef_state = logistic(0, 1), prior_intercept_det = logistic(0, 1), prior_coef_det = logistic(0, 1), prior_sigma = gamma(1, 1), log_lik = TRUE, ... )
formula |
Double right-hand side formula describing covariates of detection and occupancy in that order |
data |
A |
prior_intercept_state |
Prior distribution for the intercept of the
state (occupancy probability) model; see |
prior_coef_state |
Prior distribution for the regression coefficients of the state model |
prior_intercept_det |
Prior distribution for the intercept of the detection probability model |
prior_coef_det |
Prior distribution for the regression coefficients of the detection model |
prior_sigma |
Prior distribution on random effect standard deviations |
log_lik |
If |
... |
Arguments passed to the |
ubmsFitOccu object describing the model fit.
MacKenzie DI, Nichols JD, Lachman GB, Droege S, Royle JA, Langtimm CA. 2002. Estimating site occupancy rates when detection probabilities are less than one. Ecology 83: 2248-2255.
data(frogs) pferUMF <- unmarkedFrameOccu(pfer.bin) #Add some covariates siteCovs(pferUMF) <- data.frame(cov1=rnorm(numSites(pferUMF))) #Fit model (fm <- stan_occu(~1~cov1, pferUMF, chains=3, iter=300))data(frogs) pferUMF <- unmarkedFrameOccu(pfer.bin) #Add some covariates siteCovs(pferUMF) <- data.frame(cov1=rnorm(numSites(pferUMF))) #Fit model (fm <- stan_occu(~1~cov1, pferUMF, chains=3, iter=300))
Fit the occupancy model of Royle and Nichols (2003), which relates probability of detection of the species to the number of individuals available for detection at each site.
stan_occuRN( formula, data, K = 20, prior_intercept_state = normal(0, 5), prior_coef_state = normal(0, 2.5), prior_intercept_det = logistic(0, 1), prior_coef_det = logistic(0, 1), prior_sigma = gamma(1, 1), log_lik = TRUE, ... )stan_occuRN( formula, data, K = 20, prior_intercept_state = normal(0, 5), prior_coef_state = normal(0, 2.5), prior_intercept_det = logistic(0, 1), prior_coef_det = logistic(0, 1), prior_sigma = gamma(1, 1), log_lik = TRUE, ... )
formula |
Double right-hand side formula describing covariates of detection and abundance in that order |
data |
A |
K |
Integer upper index of integration for N-mixture. This should be set high enough so that it does not affect the parameter estimates. Note that computation time will increase with K. |
prior_intercept_state |
Prior distribution for the intercept of the
state (abundance) model; see |
prior_coef_state |
Prior distribution for the regression coefficients of the state model |
prior_intercept_det |
Prior distribution for the intercept of the detection probability model |
prior_coef_det |
Prior distribution for the regression coefficients of the detection model |
prior_sigma |
Prior distribution on random effect standard deviations |
log_lik |
If |
... |
Arguments passed to the |
ubmsFitOccuRN object describing the model fit.
Royle JA, Nichols JD. 2003. Estimating abundance from repeated presence-absence data or point counts. Ecology 84: 777-790.
data(birds) woodthrushUMF <- unmarkedFrameOccu(woodthrush.bin) #Add a site covariate siteCovs(woodthrushUMF) <- data.frame(cov1=rnorm(numSites(woodthrushUMF))) (fm_wood <- stan_occuRN(~1~cov1, woodthrushUMF, chains=3, iter=300))data(birds) woodthrushUMF <- unmarkedFrameOccu(woodthrush.bin) #Add a site covariate siteCovs(woodthrushUMF) <- data.frame(cov1=rnorm(numSites(woodthrushUMF))) (fm_wood <- stan_occuRN(~1~cov1, woodthrushUMF, chains=3, iter=300))
Fit time-to-detection occupancy models of Garrard et al. (2008, 2013). Time-to-detection can be modeled with either an exponential or Weibull distribution.
stan_occuTTD( psiformula = ~1, gammaformula = ~1, epsilonformula = ~1, detformula = ~1, data, ttdDist = c("exp", "weibull"), linkPsi = c("logit"), prior_intercept_state = logistic(0, 1), prior_coef_state = logistic(0, 1), prior_intercept_det = normal(0, 5), prior_coef_det = normal(0, 2.5), prior_intercept_shape = normal(0, 2.5), prior_sigma = gamma(1, 1), log_lik = TRUE, ... )stan_occuTTD( psiformula = ~1, gammaformula = ~1, epsilonformula = ~1, detformula = ~1, data, ttdDist = c("exp", "weibull"), linkPsi = c("logit"), prior_intercept_state = logistic(0, 1), prior_coef_state = logistic(0, 1), prior_intercept_det = normal(0, 5), prior_coef_det = normal(0, 2.5), prior_intercept_shape = normal(0, 2.5), prior_sigma = gamma(1, 1), log_lik = TRUE, ... )
psiformula |
Right-hand sided formula for the initial probability of occupancy at each site. |
gammaformula |
Right-hand sided formula for colonization probability. Currently ignored as dynamic models are not yet supported. |
epsilonformula |
Right-hand sided formula for extinction probability. Currently ignored as dynamic models are not yet supported. |
detformula |
Right-hand sided formula for mean time-to-detection. |
data |
|
ttdDist |
Distribution to use for time-to-detection; either
|
linkPsi |
Link function for the occupancy model. Only option is
|
prior_intercept_state |
Prior distribution for the intercept of the
state (occupancy probability) model; see |
prior_coef_state |
Prior distribution for the regression coefficients of the state model |
prior_intercept_det |
Prior distribution for the intercept of the time-to-detection model |
prior_coef_det |
Prior distribution for the regression coefficients of the time-to-detection model |
prior_intercept_shape |
Prior distribution for the intercept of the shape parameter (i.e., log(shape)) for Weibull TTD models |
prior_sigma |
Prior distribution on random effect standard deviations |
log_lik |
If |
... |
Arguments passed to the |
ubmsFitOccuTTD object describing the model fit.
Garrard, G.E., Bekessy, S.A., McCarthy, M.A. and Wintle, B.A. 2008. When have we looked hard enough? A novel method for setting minimum survey effort protocols for flora surveys. Austral Ecology 33: 986-998.
Garrard, G.E., McCarthy, M.A., Williams, N.S., Bekessy, S.A. and Wintle, B.A. 2013. A general model of detectability using species traits. Methods in Ecology and Evolution 4: 45-52.
Kery, Marc, and J. Andrew Royle. 2016. Applied Hierarchical Modeling in Ecology, Volume 1. Academic Press.
#Simulate data N <- 500; J <- 1 scovs <- data.frame(elev=c(scale(runif(N, 0,100))), forest=runif(N,0,1), wind=runif(N,0,1)) beta_psi <- c(-0.69, 0.71, -0.5) psi <- plogis(cbind(1, scovs$elev, scovs$forest) %*% beta_psi) z <- rbinom(N, 1, psi) Tmax <- 10 #Same survey length for all observations beta_lam <- c(-2, -0.2, 0.7) rate <- exp(cbind(1, scovs$elev, scovs$wind) %*% beta_lam) ttd <- rexp(N, rate) ttd[z==0] <- Tmax #Censor at unoccupied sites ttd[ttd>Tmax] <- Tmax #Censor when ttd was greater than survey length #Build unmarkedFrame umf <- unmarkedFrameOccuTTD(y=ttd, surveyLength=Tmax, siteCovs=scovs) #Fit model (fit <- stan_occuTTD(psiformula=~elev+forest, detformula=~elev+wind, data=umf, chains=3, iter=300))#Simulate data N <- 500; J <- 1 scovs <- data.frame(elev=c(scale(runif(N, 0,100))), forest=runif(N,0,1), wind=runif(N,0,1)) beta_psi <- c(-0.69, 0.71, -0.5) psi <- plogis(cbind(1, scovs$elev, scovs$forest) %*% beta_psi) z <- rbinom(N, 1, psi) Tmax <- 10 #Same survey length for all observations beta_lam <- c(-2, -0.2, 0.7) rate <- exp(cbind(1, scovs$elev, scovs$wind) %*% beta_lam) ttd <- rexp(N, rate) ttd[z==0] <- Tmax #Censor at unoccupied sites ttd[ttd>Tmax] <- Tmax #Censor when ttd was greater than survey length #Build unmarkedFrame umf <- unmarkedFrameOccuTTD(y=ttd, surveyLength=Tmax, siteCovs=scovs) #Fit model (fit <- stan_occuTTD(psiformula=~elev+forest, detformula=~elev+wind, data=umf, chains=3, iter=300))
This function fits the single season N-mixture model of Royle et al. (2004).
stan_pcount( formula, data, K = NULL, mixture = "P", prior_intercept_state = normal(0, 5), prior_coef_state = normal(0, 2.5), prior_intercept_det = logistic(0, 1), prior_coef_det = logistic(0, 1), prior_sigma = gamma(1, 1), log_lik = TRUE, ... )stan_pcount( formula, data, K = NULL, mixture = "P", prior_intercept_state = normal(0, 5), prior_coef_state = normal(0, 2.5), prior_intercept_det = logistic(0, 1), prior_coef_det = logistic(0, 1), prior_sigma = gamma(1, 1), log_lik = TRUE, ... )
formula |
Double right-hand side formula describing covariates of detection and abundance in that order |
data |
A |
K |
Integer upper index of integration for N-mixture. This should be set high enough so that it does not affect the parameter estimates. Note that computation time will increase with K. |
mixture |
Character specifying mixture: "P" is only option currently. |
prior_intercept_state |
Prior distribution for the intercept of the
state (abundance) model; see |
prior_coef_state |
Prior distribution for the regression coefficients of the state model |
prior_intercept_det |
Prior distribution for the intercept of the detection probability model |
prior_coef_det |
Prior distribution for the regression coefficients of the detection model |
prior_sigma |
Prior distribution on random effect standard deviations |
log_lik |
If |
... |
Arguments passed to the |
ubmsFitPcount object describing the model fit.
Royle JA. 2004. N-mixture models for estimating populaiton size from spatially replicated counts. Biometrics 60: 105-108.
data(mallard) mallardUMF <- unmarkedFramePCount(mallard.y, siteCovs=mallard.site) (fm_mallard <- stan_pcount(~1~elev+forest, mallardUMF, K=30, chains=3, iter=300))data(mallard) mallardUMF <- unmarkedFramePCount(mallard.y, siteCovs=mallard.site) (fm_mallard <- stan_pcount(~1~elev+forest, mallardUMF, K=30, chains=3, iter=300))
Extract Summary Statistics from a ubmsFit Model
## S4 method for signature 'ubmsFit' summary(object, submodel, ...)## S4 method for signature 'ubmsFit' summary(object, submodel, ...)
object |
A |
submodel |
Name of submodel to summarize |
... |
Currently ignored |
An object of class data.frame containing summary statistics
for posterior distributions of parameters from the chosen submodel.
Draws traceplots for chains from a ubmsFit object
## S4 method for signature 'ubmsFit' traceplot(object, ...)## S4 method for signature 'ubmsFit' traceplot(object, ...)
object |
A |
... |
Arguments passed to |
A ggplot object.
Generate posterior draws of turnover probability from dynamic occupancy models. Turnover is calculated for each site and each primary period after the first.
turnover(object, ...) ## S4 method for signature 'ubmsFitColext' turnover(object, draws, re.form = NULL, ...)turnover(object, ...) ## S4 method for signature 'ubmsFitColext' turnover(object, draws, re.form = NULL, ...)
object |
A fitted dynamic occupancy model of class inheriting |
... |
Currently ignored |
draws |
Number of draws from the posterior to use in the check |
re.form |
If |
A matrix of turnover values from the posterior predictive distribution.
The dimensions are draws by number of sites x (primary periods - 1)
in site-major order.
Extractors for ubmsFitList objects Extract parts of ubmsFitList objects.
## S4 method for signature 'ubmsFitList' x$name ## S4 method for signature 'ubmsFitList,numeric,missing' x[[i]] ## S4 method for signature 'ubmsFitList,numeric,missing,missing' x[i]## S4 method for signature 'ubmsFitList' x$name ## S4 method for signature 'ubmsFitList,numeric,missing' x[[i]] ## S4 method for signature 'ubmsFitList,numeric,missing,missing' x[i]
x |
A list of |
name, i
|
The names or indices of |
A ubmsFit object or list of such objects.
Widely Applicable Information Criterion (WAIC)
## S4 method for signature 'ubmsFit' waic(x, ...)## S4 method for signature 'ubmsFit' waic(x, ...)
x |
A |
... |
Currently ignored |
An object of class waic containing an estimate of WAIC and
other parameters useful for model comparison. See ?loo::waic for
more information.