fit_pompp {pompp}

Fit presence-only data using the Presence-Only for Marked Point Process model

Description

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

Usage

fit_pompp(
  object,
  background,
  mcmc_setup = list(iter = 5000),
  verbose = TRUE,
  ...
)

## S4 method for signature 'pompp_model,matrix'
fit_pompp(
  object,
  background,
  neighborhoodSize = 20,
  mcmc_setup,
  verbose = TRUE,
  area = 1,
  cores = parallel::detectCores(),
  ...
)

checkFormatBackground(object, background)

Arguments

object

Either a pompp_model or pompp_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 pompp_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.

neighborhoodSize

Number of neighbors to use in the NNGP method.

area

A positive number with the studied region's area.

cores

Number of cores to pass to OpenMP.

Details

The background is kept outside of the

Value

An object of class "pompp_fit".

See Also

pompp_model and pompp_fit-class.

Examples

beta <- c(-1, 2) # Intercept = -1. Only one covariate
delta <- c(3, 4) # Intercept = 3. Only one covariate
lambdaStar <- 1000
gamma <- 2
mu <- 5

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)]
S <- MASS::mvrnorm(1, rep(0, n_occurrences),
  2 * exp(-as.matrix(dist(occurrences_points)) / 0.3))

# Find the presence-only observations.
marks <- exp(mu + S + rnorm(n_occurrences, 0, 0.3))
po_sightings <- log(runif(n_occurrences)) <= -log1p(exp(-delta[1] - delta[2] * W1 - gamma * S))
n_po <- sum(po_sightings)
po_points <- occurrences_points[po_sightings, ]
po_Z <- occurrences_Z[po_sightings]
po_W <- W1[po_sightings]
po_marks <- marks[po_sightings]

jointPrior <- prior(
  NormalPrior(rep(0, 2), 10 * diag(2)), # Beta
NormalPrior(rep(0, 3), 10 * diag(3)), # Delta
GammaPrior(0.00001, 0.00001), # LambdaStar
NormalPrior(0, 100), GammaPrior(0.001, 0.001) # Marks
)

model <- pompp_model(po = cbind(po_Z, po_W, po_points, po_marks),
intensitySelection = 1, observabilitySelection = 2, marksSelection = 5,
                    coordinates = 3:4,
                    intensityLink = "logit", observabilityLink = "logit",
                    initial_values = 2, joint_prior = jointPrior)

bkg <- cbind(Z2, W2, reg_grid) # Create background

# Be prepared to wait a long time for this
fit <- fit_pompp(model, bkg, neighborhoodSize = 20, area = 1,
  mcmc_setup = list(burnin = 1000, iter = 2000), cores = 1)

summary(fit)

# Rhat upper CI values are above 1.1. More iterations are needed...