fit_bayesPO {bayesPO}

Fit presence-only data using a Bayesian Poisson Process model

Description

The model uses a data augmentation scheme to avoid performing approximations on the likelihood function.

Usage

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,
  ...
)

Arguments

object

Either a bayesPO_model or bayesPO_fit object. If a model, then the model is fit according to specifications. If a fit, then the model used to fit the model is recovered and used to continue the MCMC calculations where the previous one left off.

background

A matrix where the rows are the grid cells for the studied region and the columns are the covariates. NAs must be removed. If the function is being used on a bayesPO_fit object, the background must be exactly the same as the one used in the original fit.

mcmc_setup

A list containing iter to inform the model how many iterations are to be run. The list may optionally contain the objects.

verbose

Set to FALSE to suppress all messages to console.

...

Parameters passed on to specific methods. burnin and thin to inform these instructions as well.

area

A positive number with the studied region's area.

cores

Currently unused.

Details

The background is kept outside of the

Value

An object of class "bayesPO_fit".

See Also

bayesPO_model and bayesPO_fit-class.

Examples

# 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 presence-only 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)