| Title: | Bayesian Inference for Neyman-Scott Point Processes |
|---|---|
| Description: | The Bayesian MCMC estimation of parameters for Thomas-type cluster point process with various inhomogeneities. It allows for inhomogeneity in (i) distribution of parent points, (ii) mean number of points in a cluster, (iii) cluster spread. The package also allows for the Bayesian MCMC algorithm for the homogeneous generalized Thomas process. The cluster size is allowed to have a variance that is greater or less than the expected value (cluster sizes are over or under dispersed). Details are described in Dvořák, Remeš, Beránek & Mrkvička (2022) <arXiv: 10.48550/arXiv.2205.07946>. |
| Authors: | Mrkvicka Tomas [aut], Dvorak Jiri [aut], Beranek Ladislav [aut], Remes Radim [aut, cre] |
| Maintainer: | Remes Radim <[email protected]> |
| License: | GPL-3 |
| Version: | 0.1.26 |
| Built: | 2025-03-12 04:46:38 UTC |
| Source: | https://github.com/tomasmrkvicka/binspp |
The Bayesian MCMC estimation of parameters for Thomas-type cluster point process with various inhomogeneities. It allows for inhomogeneity in (i) distribution of parent points, (ii) mean number of points in a cluster, (iii) cluster spread. The package also allows for the Bayesian MCMC algorithm for the homogeneous generalized Thomas process. The cluster size is allowed to have a variance that is greater or less than the expected value (cluster sizes are over or under dispersed). Details are described in Dvořák, Remeš, Beránek & Mrkvička (2022) (doi:10.48550/arXiv.2205.07946).
License: GPL-3
Tomas Mrkvicka <[email protected]> (author), Jiri Dvorak <[email protected]> (author), Ladislav Beranek <[email protected]> (author), Radim Remes <[email protected]> (author, creator)
Anderson, C. Mrkvička T. (2020). Inference for cluster point processes with over- or under-dispersed cluster sizes, Statistics and computing 30, 1573–1590, doi:10.1007/s11222-020-09960-8.
Kopecký J., Mrkvička T. (2016). On the Bayesian estimation for the stationary Neyman-Scott point processes, Applications of Mathematics 61/4, 503-514. Available from: https://am.math.cas.cz/am61-4/9.html.
Dvořák, J., Remeš, R., Beránek, L., Mrkvička, T. (2022). binspp: An R Package for Bayesian Inference for Neyman-Scott Point Processes with Complex Inhomogeneity Structure. arXiv. doi:10.48550/ARXIV.2205.07946.
Covariate for data trees_N4.
cov_reforcov_refor
An object of class im with 208 rows and 374 columns.
Covariate for data trees_N4.
cov_reservcov_reserv
An object of class im with 208 rows and 374 columns.
Covariate for data trees_N4.
cov_slopecov_slope
An object of class im with 208 rows and 374 columns.
Covariate for data trees_N4.
cov_tdensitycov_tdensity
An object of class im with 208 rows and 374 columns.
Covariate for data trees_N4.
cov_tmicov_tmi
An object of class im with 208 rows and 374 columns.
Bayesian MCMC estimation of parameters of generalized Thomas process. The cluster size is allowed to have a variance that is greater or less than the expected value (cluster sizes are over or under dispersed).
estgtp( X, kappa0 = exp(a_kappa + ((b_kappa^2)/2)), omega0 = exp(a_omega + ((b_omega^2)/2)), lambda0 = (l_lambda + u_lambda)/2, theta0 = exp(a_theta + ((b_theta^2)/2)), skappa, somega, dlambda, stheta, smove, a_kappa, b_kappa, a_omega, b_omega, l_lambda, u_lambda, a_theta, b_theta, iter = 5e+05, plot.step = 1000, save.step = 1000, filename )estgtp( X, kappa0 = exp(a_kappa + ((b_kappa^2)/2)), omega0 = exp(a_omega + ((b_omega^2)/2)), lambda0 = (l_lambda + u_lambda)/2, theta0 = exp(a_theta + ((b_theta^2)/2)), skappa, somega, dlambda, stheta, smove, a_kappa, b_kappa, a_omega, b_omega, l_lambda, u_lambda, a_theta, b_theta, iter = 5e+05, plot.step = 1000, save.step = 1000, filename )
X |
A point pattern dataset (object of class ppp) to which the model should be fitted. |
kappa0 |
Initial value for kappa, by default it will be set as expectation of prior for kappa. |
omega0 |
Initial value for omega, by default it will be set as expectation of prior for omega. |
lambda0 |
Initial value for lambda, by default it will be set as expectation of prior for lambda. |
theta0 |
Initial value for theta, by default it will be set as expectation of prior for theta. |
skappa |
variability of proposal for kappa: second parameter of log-normal distribution |
somega |
variability of proposal for omega: second parameter of log-normal distribution |
dlambda |
variability of proposal for lambda: half of range of uniform distribution |
stheta |
variability of proposal for theta: second parameter of log-normal distribution |
smove |
variability of proposal for moving center point: SD of normal distribution |
a_kappa |
First parameter of prior distribution for kappa, which is log-normal distribution. |
b_kappa |
Second parameter of prior distribution for kappa, which is log-normal distribution. |
a_omega |
First parameter of prior distribution for omega, which is log-normal distribution. |
b_omega |
Second parameter of prior distribution for omega, which is log-normal distribution. |
l_lambda |
First parameter of prior distribution for lambda, which is uniform distribution. |
u_lambda |
Second parameter of prior distribution for lambda, which is uniform distribution. |
a_theta |
First parameter of prior distribution for theta, which is log-normal distribution. |
b_theta |
Second parameter of prior distribution for theta, which is log-normal distribution. |
iter |
Number of iterations of MCMC. |
plot.step |
Step for the graph plotting. If the value is greater than iter parameter value, no plots will be visible. |
save.step |
Step for the parameters saving. The file must be specified or has to be set to larger than iter. |
filename |
The name of the output RDS file |
The output is an estimated MCMC chain of parameters, centers and connections.
library(spatstat) kappa = 10 omega = .1 lambda= .5 theta = 10 X = rgtp(kappa, omega, lambda, theta, win = owin(c(0, 1), c(0, 1))) plot(X$X) plot(X$C) a_kappa = 4 b_kappa = 1 x <- seq(0, 100, length = 100) hx <- dlnorm(x, a_kappa, b_kappa) plot(x, hx, type = "l", lty = 1, xlab = "x value", ylab = "Density", main = "Prior") a_omega = -3 b_omega = 1 x <- seq(0, 1, length = 100) hx <- dlnorm(x, a_omega, b_omega) plot(x, hx, type = "l", lty = 1, xlab = "x value", ylab = "Density", main = "Prior") l_lambda = -1 u_lambda = 0.99 x <- seq(-1, 1, length = 100) hx <- dunif(x, l_lambda, u_lambda) plot(x, hx, type = "l", lty = 1, xlab = "x value", ylab = "Density", main = "Prior") a_theta = 4 b_theta = 1 x <- seq(0, 100, length = 100) hx <- dlnorm(x, a_theta, b_theta) plot(x, hx, type = "l", lty = 1, xlab = "x value", ylab = "Density", main = "Prior") est = estgtp(X$X, skappa = exp(a_kappa + ((b_kappa ^ 2) / 2)) / 100, somega = exp(a_omega + ((b_omega ^ 2) / 2)) / 100, dlambda = 0.01, stheta = exp(a_theta + ((b_theta ^ 2) / 2)) / 100, smove = 0.1, a_kappa = a_kappa, b_kappa = b_kappa, a_omega = a_omega, b_omega = b_omega, l_lambda = l_lambda, u_lambda = u_lambda, a_theta = a_theta, b_theta = b_theta, iter = 50, plot.step = 50, save.step = 1e9, filename = "")library(spatstat) kappa = 10 omega = .1 lambda= .5 theta = 10 X = rgtp(kappa, omega, lambda, theta, win = owin(c(0, 1), c(0, 1))) plot(X$X) plot(X$C) a_kappa = 4 b_kappa = 1 x <- seq(0, 100, length = 100) hx <- dlnorm(x, a_kappa, b_kappa) plot(x, hx, type = "l", lty = 1, xlab = "x value", ylab = "Density", main = "Prior") a_omega = -3 b_omega = 1 x <- seq(0, 1, length = 100) hx <- dlnorm(x, a_omega, b_omega) plot(x, hx, type = "l", lty = 1, xlab = "x value", ylab = "Density", main = "Prior") l_lambda = -1 u_lambda = 0.99 x <- seq(-1, 1, length = 100) hx <- dunif(x, l_lambda, u_lambda) plot(x, hx, type = "l", lty = 1, xlab = "x value", ylab = "Density", main = "Prior") a_theta = 4 b_theta = 1 x <- seq(0, 100, length = 100) hx <- dlnorm(x, a_theta, b_theta) plot(x, hx, type = "l", lty = 1, xlab = "x value", ylab = "Density", main = "Prior") est = estgtp(X$X, skappa = exp(a_kappa + ((b_kappa ^ 2) / 2)) / 100, somega = exp(a_omega + ((b_omega ^ 2) / 2)) / 100, dlambda = 0.01, stheta = exp(a_theta + ((b_theta ^ 2) / 2)) / 100, smove = 0.1, a_kappa = a_kappa, b_kappa = b_kappa, a_omega = a_omega, b_omega = b_omega, l_lambda = l_lambda, u_lambda = u_lambda, a_theta = a_theta, b_theta = b_theta, iter = 50, plot.step = 50, save.step = 1e9, filename = "")
Calculates median values for kappa, omega, lambda, theta; calculates 2.5 and 97.5 quantile and draws trace plots.
estgtpr(est, discard = 100, step = 10)estgtpr(est, discard = 100, step = 10)
est |
Output from |
discard |
Number of iterations to be discarded as burn in for the estimation. |
step |
Every step iteration is taken in the parameter estimation. |
Median and quantile values and plots (kappa, omega, lambda, theta).
library(spatstat) kappa = 10 omega = .1 lambda= .5 theta = 10 X = rgtp(kappa, omega, lambda, theta, win = owin(c(0, 1), c(0, 1))) plot(X$X) plot(X$C) a_kappa = 4 b_kappa = 1 x <- seq(0, 100, length = 100) hx <- dlnorm(x, a_kappa, b_kappa) plot(x, hx, type = "l", lty = 1, xlab = "x value", ylab = "Density", main = "Prior") a_omega = -3 b_omega = 1 x <- seq(0, 1, length = 100) hx <- dlnorm(x, a_omega, b_omega) plot(x, hx, type = "l", lty = 1, xlab = "x value", ylab = "Density", main = "Prior") l_lambda = -1 u_lambda = 0.99 x <- seq(-1, 1, length = 100) hx <- dunif(x, l_lambda, u_lambda) plot(x, hx, type = "l", lty = 1, xlab = "x value", ylab = "Density", main = "Prior") a_theta = 4 b_theta = 1 x <- seq(0, 100, length = 100) hx <- dlnorm(x, a_theta, b_theta) plot(x, hx, type = "l", lty = 1, xlab = "x value", ylab = "Density", main = "Prior") est = estgtp(X$X, skappa = exp(a_kappa + ((b_kappa ^ 2) / 2)) / 100, somega = exp(a_omega + ((b_omega ^ 2) / 2)) / 100, dlambda = 0.01, stheta = exp(a_theta + ((b_theta ^ 2) / 2)) / 100, smove = 0.1, a_kappa = a_kappa, b_kappa = b_kappa, a_omega = a_omega, b_omega = b_omega, l_lambda = l_lambda, u_lambda = u_lambda, a_theta = a_theta, b_theta = b_theta, iter = 50, plot.step = 50, save.step = 1e9, filename = "") discard = 10 step = 10 result = estgtpr(est, discard, step)library(spatstat) kappa = 10 omega = .1 lambda= .5 theta = 10 X = rgtp(kappa, omega, lambda, theta, win = owin(c(0, 1), c(0, 1))) plot(X$X) plot(X$C) a_kappa = 4 b_kappa = 1 x <- seq(0, 100, length = 100) hx <- dlnorm(x, a_kappa, b_kappa) plot(x, hx, type = "l", lty = 1, xlab = "x value", ylab = "Density", main = "Prior") a_omega = -3 b_omega = 1 x <- seq(0, 1, length = 100) hx <- dlnorm(x, a_omega, b_omega) plot(x, hx, type = "l", lty = 1, xlab = "x value", ylab = "Density", main = "Prior") l_lambda = -1 u_lambda = 0.99 x <- seq(-1, 1, length = 100) hx <- dunif(x, l_lambda, u_lambda) plot(x, hx, type = "l", lty = 1, xlab = "x value", ylab = "Density", main = "Prior") a_theta = 4 b_theta = 1 x <- seq(0, 100, length = 100) hx <- dlnorm(x, a_theta, b_theta) plot(x, hx, type = "l", lty = 1, xlab = "x value", ylab = "Density", main = "Prior") est = estgtp(X$X, skappa = exp(a_kappa + ((b_kappa ^ 2) / 2)) / 100, somega = exp(a_omega + ((b_omega ^ 2) / 2)) / 100, dlambda = 0.01, stheta = exp(a_theta + ((b_theta ^ 2) / 2)) / 100, smove = 0.1, a_kappa = a_kappa, b_kappa = b_kappa, a_omega = a_omega, b_omega = b_omega, l_lambda = l_lambda, u_lambda = u_lambda, a_theta = a_theta, b_theta = b_theta, iter = 50, plot.step = 50, save.step = 1e9, filename = "") discard = 10 step = 10 result = estgtpr(est, discard, step)
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).
estintp( X, control, x_left, x_right, y_bottom, y_top, W_dil, z_beta = NULL, z_alpha = NULL, z_omega = NULL, verbose = TRUE )estintp( X, control, x_left, x_right, y_bottom, y_top, W_dil, z_beta = NULL, z_alpha = NULL, z_omega = NULL, verbose = TRUE )
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 |
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. |
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.
The model for the intensity function of the parent process is the following:
,
where is the vector of parameters
and is the list of covariates.
Note that choosing 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
. 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 is the following:
,
where is the vector of parameters and
is the list of covariates.
Note that choosing 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 is the following:
,
where is the vector of parameters and is the list of covariates.
Note that choosing is acceptable, resulting in a constant model.
In such a case z_omega must be an empty list or NULL.
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.
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.
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.
The prior distribution for alpha is normal with
mean = Prior_alpha_mean and SD = Prior_alpha_SD.
The prior distribution for the vector 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
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)).
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.
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)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)
For exploratory purposes it may be useful to perform the first step of the analysis only, to investigate the dependence of the intensity function of the parent process on given covariates, without running the MCMC chain.
first_step(X, z_beta, W_dil, plot = TRUE)first_step(X, z_beta, W_dil, plot = TRUE)
X |
observed point pattern in the |
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. |
W_dil |
the observation window dilated by the assumed maximal cluster radius. |
plot |
logical, should the estimates intensity function of the parent process be plotted? |
The calling the spatstat.model::ppm() function from the spatstat
package, with some additional computations useful when preparing
the run of the MCMC chain, is mainly performed in this function.
The function also contains a simple way to plot the estimated
intensity function of the parent process.
List containing the output of the spatstat.model::ppm()
function from the spatstat package, along with some auxiliary
objects useful for running the MCMC chain.
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) # Determine the union of rectangles: 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) # Estimating the intensity function of the parent process: aux = first_step(X, z_beta, W_dil, plot = TRUE)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) # Determine the union of rectangles: 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) # Estimating the intensity function of the parent process: aux = first_step(X, z_beta, W_dil, plot = TRUE)
A graphical representation of the posterior distributions in terms of histograms and trace plots.
plot_outputs(Output)plot_outputs(Output)
Output |
list, output of the main function |
If the covariate list z_beta was non-empty, the estimated
intensity function of the parent process is plotted. Then,
the estimated surface representing the location dependent
mean number of points in a cluster is plotted, and similarly,
the estimated surface representing the location dependent scale
of clusters is plotted.
After that, histograms of the sample posterior distributions
of the individual parameters are plotted, together with
the histograms of p-values giving significance of the individual
covariates in z_beta with respect to the population
of parent points.
Then, the trace plots for individual model parameters are plotted,
with highlighted sample median (full red line) and sample 2.5%
and 97.5% quantiles (dashed red lines), and similarly for
the p-values giving significance of the individual covariates
in z_beta with respect to the population of parent points.
Additionally, the following graphs are also plotted:
trace plot for the log-likelihood of the model,
trace plot for the number of parent points,
trace plot for the probability of accepting proposed updates of ,
trace plot for the fraction of accepted updates of in the last 1000 iterations,
trace plot for the probability of accepting proposed updates of ,
trace plot for the fraction of accepted updates of in the last 1000 iterations,
trace plot for the fraction of accepted updates of parent points in the last 1000 iterations.
Series of plots providing a graphical representation of the posterior distributions in terms of histograms and trace plots.
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 = 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) # Default parameters for prior distributions: control = list(NStep = 100, BurnIn = 20, SamplingFreq = 5) # 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)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 = 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) # Default parameters for prior distributions: control = list(NStep = 100, BurnIn = 20, SamplingFreq = 5) # 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)
The summaries of the posterior distributions in the text form are provided.
print_outputs(Output)print_outputs(Output)
Output |
list, output of the main function |
The parameter estimates (sample medians
from the empirical posterior distributions) and the 2.5%
and 97.5% quantiles from the empirical posterior
distributions are printed.
Additionally, during the run of the MCMC chain the significance
of the covariates in the list z_beta with respect to the
current population of parent points is repeatedly tested.
This function prints the medians of the series of p-values
obtained in this way, together with the corresponding
2.5% and 97.5% sample quantiles of the p-values for each covariate.
Text output summarizing the posterior distributions.
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 = 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) # Default parameters for prior distributions: control = list(NStep = 100, BurnIn = 20, SamplingFreq = 5) # 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 print_outputs(Output)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 = 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) # Default parameters for prior distributions: control = list(NStep = 100, BurnIn = 20, SamplingFreq = 5) # 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 print_outputs(Output)
After running the MCMC chain for the given number of steps, the trace plots may indicate that too small value of burn-in was used in the first place. This function enables re-estimating the posterior distributions with a different value of burn-in, without the need to run the MCMC chain again.
re_estimate(Output, BurnIn = 0)re_estimate(Output, BurnIn = 0)
Output |
list, output of the main function estintp. |
BurnIn |
new value of burn-in. |
The output of the main function binspp contains all
the intermediate states of the chain (sampled with the required
frequency) no matter what the original value of burn-in was.
This enables simple and quick re-estimation of the posterior
distributions with either higher or lower value of burn-in
than the one used originally. The output of this function has
the same structure as the output of the main function estintp().
List containing the parameter estimates along with the 2.5% and 97.5% quantiles of the posterior distributions, along with auxiliary objects needed for printing and plotting the outputs.
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 = 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) # Default parameters for prior distributions: control = list(NStep = 100, BurnIn = 50, SamplingFreq = 5) # 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) # Recompute the outputs when another value of burn-in is desired, # without running the chain again: Out2 <- re_estimate(Output, BurnIn = 80) print_outputs(Out2) plot_outputs(Out2)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 = 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) # Default parameters for prior distributions: control = list(NStep = 100, BurnIn = 50, SamplingFreq = 5) # 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) # Recompute the outputs when another value of burn-in is desired, # without running the chain again: Out2 <- re_estimate(Output, BurnIn = 80) print_outputs(Out2) plot_outputs(Out2)
Simulation of generalized Thomas process.
rgtp( kappa, omega, lambda, theta, win = owin(c(0, 1), c(0, 1)), nsim = 1, expand = 4 * omega )rgtp( kappa, omega, lambda, theta, win = owin(c(0, 1), c(0, 1)), nsim = 1, expand = 4 * omega )
kappa |
intensity of cluster centers. |
omega |
standard deviation of normal distribution specifying the clusters spread. |
lambda |
parameter of generalised Poisson distribution controlling over or under dispersion. |
theta |
parameter of generalised Poisson distribution controlling the mean number of points in a cluster. |
win |
window in which to simulate the pattern. An object in the |
nsim |
number of simulations. |
expand |
the size of expansion of window to simulate the centers of clusters. |
A list(X, C), where X is Generalized Thomas process, and C is Process of cluster centers for Generalized Thomas process.
library(spatstat) kappa = 10 omega = .1 lambda= .5 theta = 10 X = rgtp(kappa, omega, lambda, theta, win = owin(c(0, 1), c(0, 1))) plot(X$X) plot(X$C)library(spatstat) kappa = 10 omega = .1 lambda= .5 theta = 10 X = rgtp(kappa, omega, lambda, theta, win = owin(c(0, 1), c(0, 1))) plot(X$X) plot(X$C)
The means to simulate realizations from the Thomas-type cluster point process with complex inhomogeneities are provided.
rThomasInhom( kappa, alpha, omega, W, W_dil, betavec = NULL, alphavec = NULL, omegavec = NULL, z_beta = NULL, z_alpha = NULL, z_omega = NULL )rThomasInhom( kappa, alpha, omega, W, W_dil, betavec = NULL, alphavec = NULL, omegavec = NULL, z_beta = NULL, z_alpha = NULL, z_omega = NULL )
kappa |
intensity or intensity function of the parent process, scalar or pixel image object of class |
alpha |
scalar, influences the mean number of points in individual clusters, see Details. |
omega |
scalar, influences the spread of individual clusters, see Details. |
W |
the observation window where the realization is to be generated, in the |
W_dil |
the observation window dilated by the assumed maximal cluster radius, as a binary mask with the same resolution as the covariates. |
betavec |
vector of parameters describing the dependence of the intensity function of the parent process on covariates in the list z_beta. |
alphavec |
vector of parameters describing the dependence of the mean number of points in a cluster on covariates in the list z_alpha. |
omegavec |
vector of parameters describing the dependence of the spread of the clusters on covariates in the list z_omega. |
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. |
A realization of a Thomas-type cluster
point process model with possible inhomogeneity
(described by covariates) are produced in any or all of the following model
components: intensity function of the parent process, mean number
of points in a cluster, scale of the clusters.
Model parametrization is described in the documentation to the
function estintp(). The parent process is generated in the dilated
observation window W_dil to avoid edge-effects,
the resulting point pattern is eventually truncated to the smaller
observation window W.
A planar point pattern, object of the type spatstat.geom::ppp()
used in the spatstat package.
library(spatstat) # Unit square observation window: W <- owin() # Dilation of the observation window: W_dil <- dilation(W, 0.1) W_dil <- as.mask(W_dil) # Define covariates: f1 <- function(x, y) { x } f2 <- function(x, y) { y } f3 <- function(x, y) { 1 - (y - 0.5) ^ 2 } cov1 <- as.im(f1, W = W_dil) cov2 <- as.im(f2, W = W_dil) cov3 <- as.im(f3, W = W_dil) # Stationary Thomas process: X <- rThomasInhom(kappa = 50, alpha = log(10), omega = log(0.01), W = W, W_dil = W_dil) plot(X) # Thomas-type cluster process with inhomogeneity in all model components: X <- rThomasInhom(kappa = 10, betavec = c(1), z_beta = list(cov1), alpha = log(10), alphavec = c(1), z_alpha = list(cov2), omega = log(0.01), omegavec = c(1), z_omega = list(cov3), W = W, W_dil = W_dil) plot(X)library(spatstat) # Unit square observation window: W <- owin() # Dilation of the observation window: W_dil <- dilation(W, 0.1) W_dil <- as.mask(W_dil) # Define covariates: f1 <- function(x, y) { x } f2 <- function(x, y) { y } f3 <- function(x, y) { 1 - (y - 0.5) ^ 2 } cov1 <- as.im(f1, W = W_dil) cov2 <- as.im(f2, W = W_dil) cov3 <- as.im(f3, W = W_dil) # Stationary Thomas process: X <- rThomasInhom(kappa = 50, alpha = log(10), omega = log(0.01), W = W, W_dil = W_dil) plot(X) # Thomas-type cluster process with inhomogeneity in all model components: X <- rThomasInhom(kappa = 10, betavec = c(1), z_beta = list(cov1), alpha = log(10), alphavec = c(1), z_alpha = list(cov2), omega = log(0.01), omegavec = c(1), z_omega = list(cov3), W = W, W_dil = W_dil) plot(X)
The oak trees dataset sampled in 2009 in region consisting of 5 rectangles.
trees_N4trees_N4
A list with columns:
A list of region window definition.
Number of oak trees in the region.
Array of x coordinates of oak trees.
Array of y coordinates of oak trees.
The data contains point pattern of trees, 5 covariates (refor, reserve, slope, tdensity, tmi) and 4 vectors (x_left, x_right, y_bottom, y_top) of corners of rectangles forming the observation window.
Jesús Fernández-Habas, Pilar Fernández-Rebollo, Mónica Rivas Casado, Alma María García Moreno, Begoña Abellanas. Spatio-temporal analysis of oak decline process in open woodlands: A case study in SW Spain, Journal of Environmental Management, 248, 2019, 109308, doi:10.1016/j.jenvman.2019.109308.
plot(trees_N4)plot(trees_N4)
The vector of left horizontal corners of rectangles forming observation window for trees_N4 dataset.
x_left_N4x_left_N4
An object of class numeric of length 5.
The vector of right horizontal corners of rectangles forming observation window for trees_N4 dataset.
x_right_N4x_right_N4
An object of class numeric of length 5.
The vector of bottom vertical corners of rectangles forming observation window for trees_N4 dataset.
y_bottom_N4y_bottom_N4
An object of class numeric of length 5.
The vector of top vertical corners of rectangles forming observation window for trees_N4 dataset.
y_top_N4y_top_N4
An object of class numeric of length 5.