# Estimation of interaction Neyman-Scott point process using auxiliary variable algorithm into Markov chain Monte Carlo. The Bayesian MCMC estimation of parameters for interaction Neyman-Scott point process using auxiliary variable algorithm into Markov chain Monte Carlo. ```r estinternsp( X, control, x_left, x_right, y_bottom, y_top, W_dil, AreaW, AreaMRW, radius ) ``` ## Arguments - `X`: observed point pattern in the `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`: 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 `dilation.owin` to create a dilated window. ## Returns 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. ## 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__mean` and `SD=Prior__SD`. During update, the Lower Bound for each parameter is `_LB` and the Upper Bound is `_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. ## Examples ```r 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) ``` ## 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 `print.output_estinternsp`

estinternsp function