estinternsp {binspp} | R Documentation |
Estimation of interaction Neyman-Scott point process using auxiliary variable algorithm into Markov chain Monte Carlo.
Description
The Bayesian MCMC estimation of parameters for interaction Neyman-Scott point process using auxiliary variable algorithm into Markov chain Monte Carlo.
Usage
estinternsp(
X,
control,
x_left,
x_right,
y_bottom,
y_top,
W_dil,
AreaW,
AreaMRW,
radius
)
Arguments
X |
observed point pattern in the |
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 |
observation window dilated by the assumed maximal cluster radius. |
AreaW |
area of a window. |
AreaMRW |
area of a dilated window. |
radius |
The radius of dilation, passed as an argument to |
Details
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 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.
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), SamplingFreq (sampling frequency for estimating the posterior distributions). Additionally, the hyperparameters for the prior distributions should be given, see below.
Prior distributions and hyperparameters
The prior distribution for each parameter (alpha, omega, kappa, theta1, theta2) is normal with respective mean=Prior_<parameter>_mean
and SD=Prior_<parameter>_SD
.
During update, the Lower Bound for each parameter is <parameter>_LB
and the Upper Bound is <parameter>_UB
.
alpha controls the expected number of occurrences or events around specified focal points or regions, omega controls the width of cluster events activity and kappa controls the overall intensity of the parent process.
theta1 controls the shape of the interaction function and theta2 controls the shape of the interaction function and provides the location of the peak value.
Value
The output of the function is given by the list containing the parameter estimates along with the 2.5% and 97.5% confidence interval of the posterior distributions, cluster centers, control, W_dil
and parameters.
References
Park, J., Chang, W., & Choi, B. (2022). An interaction Neyman-Scott point process model for coronavirus disease-19. Spatial Statistics, 47, 100561.
See Also
Examples
library(spatstat)
# library(spatstat.geom)
library(fields)
library(mvtnorm)
library(binspp)
# Generate example window
W <- owin(xrange = c(0, 100), yrange = c(0, 50)) # example window (rectangle)
radius = 2
W_dil = dilation.owin(W, radius)
AreaW <- area(W)
AreaMRW <- area(W_dil)
x_left = W$xrange[1]
x_right = W$xrange[2]
y_bottom = W$yrange[1]
y_top = W$yrange[2]
# True parameters
alpha <- 15
omega <- 1.5
kappa <- 0.001
theta1 <- 3
theta2 <- 10
# Generate parents process
CCC <- rpoispp(kappa, win = W_dil)
CC <- t(rbind(CCC$x,CCC$y))
CCdist <- rdist(CC)
r <- binspp::coeff(c(theta1,theta2))
r1 <- r[1]; r2 <- r[2]; t1 <- theta1; t2 <- theta2; t3 <- 0.5; R <- 0;
rho0sum <- rep(0,dim(CC)[1]) # row sum of rho matrix
res <- binspp::pCClik2(c(theta1,theta2), CC)
rho0sum <- res$rhosum
likelihoodprev <- res$likelihood
# this is just for example, for a real estimate
# use value of at least niter = 10000
CC.true <- AuxVarGen(kappa, c(theta1,theta2), likelihoodprev, rho0sum, CC,
AreaMRW, W_dil, 100)
# Generate offsprings given parents
gaus <- function(n, omega) {
matrix(rnorm(2 * n, mean=0, sd=omega), ncol=2)
}
parents <- ppp(CC.true[[1]][,1],CC.true[[1]][,2],window=W)
np <- npoints(parents)
csize <- qpois(runif(np, min = dpois(0, alpha)), alpha)
noff <- sum(csize)
xparent <- parents$x
yparent <- parents$y
x0 <- rep.int(xparent, csize)
y0 <- rep.int(yparent, csize)
dd <- gaus(noff, omega)
xy <- xy.coords(dd)
dx <- xy$x; dy <- xy$y
xoff <- x0 + dx; yoff <- y0 + dy
result <- ppp(xoff, yoff, window = W_dil, check = FALSE, marks = NULL)
X <- cbind(result$x,result$y) # Generate example data
# this is just for example, for a real estimate
# use values of at least NStep = 20000 and BurnIn = 5000
control = list(NStep = 50, BurnIn = 25, SamplingFreq = 10,
Prior_alpha_mean = 15, Prior_alpha_SD = 2, alpha_LB = 10, alpha_UB = 20,
Prior_omega_mean = 1.5, Prior_omega_SD = 0.2, omega_LB = 1, omega_UB = 2,
Prior_kappa_mean = 0.001, Prior_kappa_SD = 0.0001, kappa_LB = 0.0005,
kappa_UB = 0.002,
Prior_theta1_mean = 3, Prior_theta1_SD = 0.2, theta1_LB = 2.5,
theta1_UB = 3.5,
Prior_theta2_mean = 10, Prior_theta2_SD = 1, theta2_LB = 8,
theta2_UB = 12)
Output = estinternsp(X, control,
x_left, x_right, y_bottom, y_top,
W_dil,
AreaW, AreaMRW, radius)