fit_bayesPO  R Documentation 
The model uses a data augmentation scheme to avoid performing approximations on the likelihood function.
fit_bayesPO(
object,
background,
mcmc_setup = list(iter = 5000),
verbose = TRUE,
...
)
## S4 method for signature 'bayesPO_model,matrix'
fit_bayesPO(
object,
background,
mcmc_setup,
verbose = TRUE,
area = 1,
cores = 1,
...
)
## S4 method for signature 'bayesPO_fit,matrix'
fit_bayesPO(
object,
background,
mcmc_setup = list(iter = object$mcmc_setup$iter),
verbose = TRUE,
cores = 1,
...
)
object 
Either a 
background 
A matrix where the rows are the grid cells for the studied
region and the columns are the covariates. 
mcmc_setup 
A list containing 
verbose 
Set to 
... 
Parameters passed on to specific methods.

area 
A positive number with the studied region's area. 
cores 
Currently unused. 
The background is kept outside of the
An object of class "bayesPO_fit"
.
bayesPO_model
and bayesPO_fitclass
.
# This code is replicated from the vignette.
## Not run:
beta < c(1, 2) # Intercept = 1. Only one covariate
delta < c(3, 4) # Intercept = 3. Only one covariate
lambdaStar < 1000
total_points < rpois(1, lambdaStar)
random_points < cbind(runif(total_points), runif(total_points))
grid_size < 50
# Find covariate values to explain the species occurrence.
# We give them a Gaussian spatial structure.
reg_grid < as.matrix(expand.grid(seq(0, 1, len = grid_size), seq(0, 1, len = grid_size)))
Z < MASS::mvrnorm(1, rep(0, total_points + grid_size * grid_size),
3 * exp(as.matrix(dist(rbind(random_points, reg_grid))) / 0.2))
Z1 < Z[1:total_points]; Z2 < Z[(1:total_points)]
# Thin the points by comparing the retaining probabilities with uniforms
# in the log scale to find the occurrences
occurrences < log(runif(total_points)) <= log1p(exp(beta[1]  beta[2] * Z1))
n_occurrences < sum(occurrences)
occurrences_points < random_points[occurrences,]
occurrences_Z < Z1[occurrences]
# Find covariate values to explain the observation bias.
# Additionally create a regular grid to plot the covariate later.
W < MASS::mvrnorm(1, rep(0, n_occurrences + grid_size * grid_size),
2 * exp(as.matrix(dist(rbind(occurrences_points, reg_grid))) / 0.3))
W1 < W[1:n_occurrences]; W2 < W[(1:n_occurrences)]
# Find the presenceonly observations.
po_sightings < log(runif(n_occurrences)) <= log1p(exp(delta[1]  delta[2] * W1))
n_po < sum(po_sightings)
po_points < occurrences_points[po_sightings, ]
po_Z < occurrences_Z[po_sightings]
po_W < W1[po_sightings]
jointPrior < prior(
NormalPrior(rep(0, 2), 10 * diag(2)), # Beta
NormalPrior(rep(0, 2), 10 * diag(2)), # Delta
GammaPrior(0.00001, 0.00001) # LambdaStar
)
model < bayesPO_model(po = cbind(po_Z, po_W),
intensitySelection = 1, observabilitySelection = 2,
intensityLink = "logit", observabilityLink = "logit",
initial_values = 2, joint_prior = jointPrior)
bkg < cbind(Z2, W2) # Create background
fit < fit_bayesPO(model, bkg, area = 1, mcmc_setup = list(burnin = 1000, iter = 2000))
summary(fit)
# Rhat upper CI values are above 1.1. More iterations are needed, so...
fit2 < fit_bayesPO(fit, bkg, mcmc_setup = list(iter = 10000))
summary(fit2)
mcmc_trace(fit2)
mcmc_dens(fit2)
## End(Not run)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.