estintp {binspp}

Estimation of Thomas-type cluster point process with complex inhomogeneities

Description

The Bayesian MCMC estimation of parameters for Thomas-type cluster point process with inhomogeneity is performed in any of the following parts: (i) distribution of parent points, (ii) mean number of points in a cluster, (iii) cluster spread. The process is observed in the observation window W which is a union of aligned rectangles, aligned with the coordinate axes. The inhomogeneities are described through a parametric model depending on covariates. The estimation algorithm is described in Dvořák, Remeš, Beránek & Mrkvička (2022) (doi:10.48550/arXiv.2205.07946).

Usage

estintp(
  X,
  control,
  x_left,
  x_right,
  y_bottom,
  y_top,
  W_dil,
  z_beta = NULL,
  z_alpha = NULL,
  z_omega = NULL,
  verbose = TRUE
)

Arguments

X

observed point pattern in the spatstat.geom::ppp() format of the spatstat package.

control

list specifying various tuning constants for the MCMC estimation. See also Details.

x_left

vector describing the observation window, contains the lower x-coordinate of the corners of each rectangle.

x_right

vector describing the observation window, contains the higher x-coordinate of the corners of each rectangle.

y_bottom

vector describing the observation window, contains the smaller y-coordinate of the corners of each rectangle.

y_top

vector describing the observation window, contains the higher y-coordinate of the corners of each rectangle.

W_dil

the observation window dilated by the assumed maximal cluster radius.

z_beta

list of covariates describing the intensity function of the parent process, each covariate being a pixel image as used in the spatstat package.

z_alpha

list of covariates describing the location-dependent mean number of points in a cluster, each covariate being a pixel image as used in the spatstat package.

z_omega

list of covariates describing the location-dependent scale of a cluster, each covariate being a pixel image as used in the spatstat package.

verbose

logical (TRUE or FALSE). For suppressing information messages to the console set value to FALSE. Defaults to TRUE.

Value

The output of the function is given by the list containing the parameter estimates along with the 2.5% and 97.5% quantiles of the posterior distributions. Also, several auxiliary objects are included in the list which are needed for the print_outputs() and plot_outputs() functions.

Details

Parametric model

The model for the intensity function of the parent process is the following: f(u) = kappa * exp(beta_1 * z\_beta_1(u) + … + beta_k * z\_beta_k(u)), where (kappa, beta_1, …, beta_k) is the vector of parameters and z\_beta = (z\_beta_1, …, z\_beta_k) is the list of covariates. Note that choosing k = 0 is acceptable, resulting in a homogeneous distribution of parents. In such a case z_beta must be an empty list or NULL. Furthermore, the list z_beta must contain named covariates in order to properly function with the function spatstat.model::ppm() from the spatstat package which is used in the first step to estimate the parameters (kappa, beta_1, …, beta_k). Note that due to identifiability issues the covariate lists z_beta and z_alpha must be disjoint.

The model for the mean number of points in a cluster corresponding to the parent at location u is the following: g(u) = exp(alpha + alpha_1 * z\_alpha_1(u) + … + alpha_l * z\_alpha_l(u)), where (alpha, alpha_1, …, alpha_l) is the vector of parameters and z\_alpha = (z\_alpha_1, …, z\_alpha_l) is the list of covariates. Note that choosing l = 0 is acceptable, resulting in a constant model. In such a case z_alpha must be an empty list or NULL. Note that due to identifiability issues the covariate lists z_beta and z_alpha must be disjoint.

The model for the scale of a cluster corresponding to the parent at location u is the following: h(u) = exp(omega + omega_1 * z\_omega_1(u) + … + omega_m * z\_omega_m(u)), where
(omega, omega_1, …, omega_m) is the vector of parameters and
z\_omega = (z\_omega_1, …, z\_omega_m) is the list of covariates. Note that choosing m = 0 is acceptable, resulting in a constant model. In such a case z_omega must be an empty list or NULL.

Observation window and its dilation

The observation window must be provided as the union of aligned rectangles, aligned with the coordinate axes. This, however, allows the analysis of point patterns observed in rather irregular regions by approximating the region by a union of aligned rectangles. The structure of the vectors x_left, x_right, y_bottom and y_top is such that the first rectangle is constructed using the function spatstat.geom::owin() from the spatstat package as owin(c(x_left[1], x_right[1]), c(y_bottom[1], y_top[1])), and similarly for the other rectangles. Naturally, a rectangular window can be used and in such a case the vectors x_left to y_top each contain a single element.

Covariates

The covariates must be provided as pixel images of the class spatstat.geom::im() used in the spatstat package. It is recommended that all the covariates have the same pixel resolution. However, it is necessary that all the covariates in the list z_beta have the same resolution, all the covariates in the list z_alpha have the same resolution and all the covariates in the list z_omega have the same distribution. The covariates must be provided in the dilated observation window W_dil, with NA values at pixels lying outside W_dil.

Control

The control list must contain the following elements: NStep (the required number of MCMC iterations to be performed), BurnIn (burn-in, how many iterations at the beginning of the chain will be disregarded when computing the resulting estimates – note that this choice can be updated after the computation without re-running the chain, see the function re_estimate()), SamplingFreq (sampling frequency for estimating the posterior distributions). Additionally, the hyperparameters for the prior distributions should be given, see below. Note that some default values for the hyperparameters are provided but it is strongly encouraged that the hyperparameter values are given by the user, based on the actual knowledge of the problem at hand.

Prior distributions and hyperparameters

The prior distribution for alpha is normal with mean = Prior_alpha_mean and
SD = Prior_alpha_SD.

The prior distribution for the vector (alpha_1, …, alpha_l) is centered normal with diagonal variance matrix and the vector of SDs = Prior_alphavec_SD.

The prior distribution for omega is normal with mean = Prior_omega_mean and
SD = Prior_omega_SD.

The prior distribution for the vector (omega_1, …, omega_m) is centered normal with diagonal variance matrix and the vector of SDs = Prior_omegavec_SD.

The hyperparameters should be provided in the control list. However, the following default choices are applied if the hyperparameter values are not provided by user or are given as NULL:
Prior_alpha_mean = 3, Prior_alpha_SD = 2, Prior_omega_mean = log(sqrt(area(W) / 20)), Prior_omega_SD = log(3 + sqrt(area(W) / 40)), Prior_alphavec_SD[i] = 2 / max(z_alpha_i), Prior_omegavec_SD[i] = 2 / max(z_omega_i) * log(3 + sqrt(area(W) / 20)).

Output

The output of the function is given by the list containing the parameter estimates along with the 2.5% and 97.5% quantiles of the posterior distributions. Also, several auxiliary objects are included in the list which are needed for the print_outputs() and plot_outputs() functions.

Examples

library(spatstat)
# Prepare the dataset:
X = trees_N4
x_left = x_left_N4
x_right = x_right_N4
y_bottom = y_bottom_N4
y_top = y_top_N4

z_beta = list(refor = cov_refor, slope = cov_slope)
z_alpha = list(tmi = cov_tmi, tdensity = cov_tdensity)
z_omega = list(slope = cov_slope, reserv = cov_reserv)

# Determine the union of rectangles:
W = NULL
W = owin(c(x_left[1], x_right[1]), c(y_bottom[1], y_top[1]))
if (length(x_left) >= 2) {
  for (i in 2:length(x_left)) {
    W2 = owin(c(x_left[i], x_right[i]), c(y_bottom[i], y_top[i]))
    W = union.owin(W, W2)
  }
}

# Dilated observation window:
W_dil = dilation.owin(W, 100)


# User-specified hyperparameters for prior distributions:
control = list(NStep = 100, BurnIn = 50, SamplingFreq = 5,
    Prior_alpha_mean = 3, Prior_alpha_SD = 2, Prior_omega_mean = 5.5,
    Prior_omega_SD = 5, Prior_alphavec_SD = c(4.25, 0.012),
    Prior_omegavec_SD = c(0.18,0.009))

# MCMC estimation:
Output = estintp(X, control, x_left, x_right, y_bottom, y_top,
    W_dil, z_beta, z_alpha, z_omega, verbose = FALSE)

# Text output + series of figures:
print_outputs(Output)
plot_outputs(Output)