\documentclass[nojss]{jss}

% \VignetteIndexEntry{GET: Point patterns}
% \VignetteKeywords{global envelope test, goodness-of-fit, Monte Carlo test, R, spatial point pattern}
% \VignettePackage{GET}
% \VignetteEngine{R.rsp::rsp}

%% -- LaTeX packages and custom commands ---------------------------------------

%% recommended packages
\usepackage{thumbpdf,lmodern}

%% another package (only for this demo article)
\usepackage{framed}

%% new custom commands
\newcommand{\class}[1]{`\code{#1}'}
\newcommand{\fct}[1]{\code{#1()}}
\newcommand{\TT}{\mathbf{T}}
\newcommand{\lo}{\mathrm{low}}
\newcommand{\up}{\mathrm{upp}}
\providecommand{\1}{\mathbf{1}}
\newtheorem{definition}{Definition}[section]
\newtheorem{remark}{Remark}[section]

\usepackage{amsmath,amsfonts,amssymb}

\makeatletter
\setkeys{Gin}{width=\Gin@nat@width}
\makeatother

%% For Sweave-based articles about R packages:
%% need no \usepackage{Sweave}


%% -- Article metainformation (author, title, ...) -----------------------------

%% - \author{} with primary affiliation
%% - \Plainauthor{} without affiliations
%% - Separate authors by \And or \AND (in \author) or by comma (in \Plainauthor).
%% - \AND starts a new line, \And does not.
\author{Mari Myllym{\"a}ki\\Natural Resources Institute Finland (Luke)}
\Plainauthor{Mari Myllym{\"a}ki}

%% - \title{} in title case
%% - \Plaintitle{} without LaTeX markup (if any)
%% - \Shorttitle{} with LaTeX markup (if any), used as running title
\title{\pkg{GET}: Point patterns}
\Plaintitle{GET: Point patterns}
\Shorttitle{\pkg{GET}: Point patterns}

%% - \Abstract{} almost as usual
\Abstract{
This vignette gives examples of the use of global envelopes for point pattern analysis,
as implemented in the \proglang{R} package \pkg{GET}.
When citing the vignette and package please cite \citet{MyllymakiMrkvicka2023}
and references given by typing \code{citation("GET")} in \proglang{R}.
}

%% - \Keywords{} with LaTeX markup, at least one required
%% - \Plainkeywords{} without LaTeX markup (if necessary)
%% - Should be comma-separated and in sentence case.
\Keywords{global envelope test, goodness-of-fit, Monte Carlo test, \proglang{R}, spatial point pattern}
\Plainkeywords{global envelope test, goodness-of-fit, Monte Carlo test, R, spatial point pattern}

%% - \Address{} of at least one author
%% - May contain multiple affiliations for each author
%%   (in extra lines, separated by \emph{and}\\).
%% - May contain multiple authors for the same affiliation
%%   (in the same first line, separated by comma).
\Address{
  Mari Myllym{\"a}ki\\
  Natural Resources Institute Finland (Luke)\\
  Latokartanonkaari 9\\
  FI-00790 Helsinki, Finland\\
  E-mail: \email{mari.myllymaki@luke.fi}\\
  URL: \url{https://www.luke.fi/en/experts/mari-myllymaki/}
}

\begin{document}

\section{Introduction}

This vignette gives examples of the use of global envelopes for the analysis of spatial point patterns.
The examples utilize the \proglang{R} \citep{R2023} package \pkg{spatstat} \citep{spatstat2015} in addition to the \pkg{GET} package \citep{MyllymakiMrkvicka2023}. The envelope plots are produced by the use of the \pkg{ggplot2} package \citep{Wickham2016}, where we utilize the theme \code{theme_bw} for this document.
%
\begin{Schunk}
\begin{Sinput}
R> library("GET")
R> library("spatstat.model")
R> library("ggplot2")
R> theme_set(theme_bw(base_size = 9))
\end{Sinput}
\end{Schunk}
%

\section{General workflow of the tests}

\subsection{Utilizing the spatstat package}

In general, it is useful in the point pattern analysis utilize the \pkg{spatstat} package. The workflow utilizing spatstat with the \pkg{GET} package is typically the following:
Say we have a point pattern, for which we would like to test a hypothesis, as a \code{ppp} object of spatstat. E.g.
%
\begin{Schunk}
\begin{Sinput}
R> X <- spruces
R> X
\end{Sinput}
\begin{Soutput}
Marked planar point pattern: 134 points
marks are numeric, of storage type  'double'
window: rectangle = [0, 56] x [0, 38] metres
\end{Soutput}
\end{Schunk}
%

1. To test a simple hypothesis, e.g., complete spatial randomness (CSR):

\begin{itemize}
\item Use the function envelope of \pkg{spatstat} to create \code{nsim} simulations
under CSR and to calculate the functions you want.
Important: use the option \code{savefuns=TRUE} and
specify the number of simulations \code{nsim}.
See the help documentation in the \pkg{spatstat} package
for possible test functions (if the argument \code{fun} is not given, the function \fct{Kest} is used,
i.e. an estimator of the $K$-function).

Making 999 simulations of CSR and estimating $K$-function for each of them and data
(the argument \code{simulate} specifies how to perform simulations under CSR):
%
%
\begin{Schunk}
\begin{Sinput}
R> env <- envelope(X, nsim=1999, savefuns=TRUE,
+                  simulate=expression(runifpoint(ex=X)),
+                  verbose=FALSE)
\end{Sinput}
\end{Schunk}
%
\item Perform the test
%
\begin{Schunk}
\begin{Sinput}
R> res <- global_envelope_test(env)
\end{Sinput}
\end{Schunk}
%

\item Plot the result
%
\begin{Schunk}
\begin{Sinput}
R> plot(res)
\end{Sinput}
\end{Schunk}
\includegraphics{pointpatterns-workflow_res}
%
\end{itemize}

2. To test a goodness-of-fit of a parametric model (composite hypothesis case):

\begin{itemize}
\item Fit the model to your data by means of the function \fct{ppm} or \fct{kppm} of \pkg{spatstat}.
See the help documentation for possible models.
\item Use the function \fct{GET.composite} to create \code{nsim} simulations
from the fitted model, to calculate the functions you want,
and to make an adjusted global envelope test.
See the example below.
\item Plot the result.
\end{itemize}

More detailed examples are given below.

\subsection{The workflow when using your own programs for simulations}

\begin{itemize}
\item (Fit the model and) Create $s$ simulations from the (fitted) null model.
\item Calculate the functions $T_1(r), T_2(r), \dots, T_{s+1}(r)$.
\item Use \fct{create\_curve\_set} to create a \code{curve\_set} object from the functions $T_i(r)$, $i=1, \dots, s+1$.
\item Perform the test and plot the result
\end{itemize}

See example in the help file of the \fct{global\_envelope\_test} function.

\section{Testing simple hypotheses}

\subsection{Testing complete spatial randomness (CSR)}

Let us illustrate the CSR for the spruces data set from the R library \pkg{spatstat}.
%
\begin{Schunk}
\begin{Sinput}
R> X <- unmark(spruces)
R> par(mfrow = c(1,1), mgp = c(0, 0, 0), mar = c(0, 0, 0, 0))
R> plot(X, main = "")
\end{Sinput}
\end{Schunk}
\includegraphics{pointpatterns-spruces}
%

Below the function \fct{envelope} of the \pkg{spatstat} package is used to generate point patterns under CSR (specified in the argument \code{simulate}) and to calculate the centred $L$-functions (specified below by the arguments \code{fun}, \code{correction} and \code{transform}), which are used here as the test functions.
%
%
\begin{Schunk}
\begin{Sinput}
R> nsim <- 1999 # Number of simulations
R> env <- envelope(X, fun = "Lest", nsim = nsim,
+    savefuns = TRUE, # save the functions
+    correction = "translate", # edge correction for L
+    transform = expression(.-r), # centering
+    simulate = expression(runifpoint(ex = X)), # Simulate CSR
+    verbose = FALSE)
R> res <- global_envelope_test(env, type = "erl")
R> plot(res)
\end{Sinput}
\end{Schunk}
\includegraphics{pointpatterns-spruces_csr}
%

It is possible to cut the functions to an interval of distances $[r_{\min}, r_{\max}]$ (at the same time creating a \code{curve_set} from \code{env}) and perform the test on the functions on this interval only:
%
\begin{Schunk}
\begin{Sinput}
R> cset <- crop_curves(env, r_min = 1, r_max = 7)
R> # Do the rank envelope test (erl)
R> res <- global_envelope_test(cset, type = "erl")
R> plot(res) + ylab(expression(italic(hat(L)(r)-r)))
\end{Sinput}
\end{Schunk}
\includegraphics{pointpatterns-spruces_csr_crop}
%

\subsection{Testing random labeling of marks}

Let now the studied marked point pattern be the spruces data with marks:
%
\begin{Schunk}
\begin{Sinput}
R> mpp <- spruces
R> par(mfrow=c(1,1), mgp=c(0, 0, 0), mar=c(0, 0, 0, 0))
R> plot(mpp, main = "")
\end{Sinput}
\end{Schunk}
\includegraphics{pointpatterns-spruces_mpp}
%

As in the CSR test, to perform the test of random labelling hypothesis, first the \fct{envelope} function of \pkg{spatstat} can be used to generate simulations and calculate the test function $T(r)$ for the data pattern (mpp) and each simulation.
Below the estimator of the mark-weighted $L$-function, $L_{mm}(r)$, with translational edge correction is used as the test function. The argument \code{simulate} specifies the simulations under the random labeling, i.e., simple permutation of the marks.
%
%
\begin{Schunk}
\begin{Sinput}
R> nsim <- 1999 # Number of simulations
R> env <- envelope(mpp, fun = Kmark, nsim = nsim, f = function(m1, m2) { m1*m2 },
+    correction = "translate", returnL = TRUE,
+    simulate = expression(rlabel(mpp, permute = TRUE)), # Permute the marks
+    savefuns = TRUE, # Save the functions
+    verbose = FALSE)
\end{Sinput}
\end{Schunk}
Thereafter, the curves can be cropped to the desired $r$-interval and centered by the mean of the simulated functions for better visualization, before making the test.
%
\begin{Schunk}
\begin{Sinput}
R> # Crop curves to desired r-interval
R> curve_set <- crop_curves(env, r_min = 1.5, r_max = 9.5)
R> # Center the functions, i.e. take \hat{L}_mm(r)-the mean of simulated functions.
R> curve_set <- residual(curve_set)
R> # The global envelope test
R> res <- global_envelope_test(curve_set)
R> plot(res) + ylab(expression(italic(L[mm](r)-L(r))))
\end{Sinput}
\end{Schunk}
\includegraphics{pointpatterns-spruces_randomlabeling_test}
%

\subsection{A combined global envelope test}

Sometimes it may be desired to base the test on several test functions \citep{MrkvickaEtal2017}.
Below it is illustrated how the CSR can be tested simultaneously by means of $L$, $F$, $G$ and $J$ functions for the saplings data set available in the \pkg{GET} library.
First some setup:
%
\begin{Schunk}
\begin{Sinput}
R> data(saplings)
R> X <- as.ppp(saplings, W = square(75))
R> nsim <- 499 # Number of simulations
R> # Specify distances for different test functions
R> n <- 500 # the number of r-values
R> rmin <- 0; rmax <- 20; rstep <- (rmax-rmin)/n
R> rminJ <- 0; rmaxJ <- 8; rstepJ <- (rmaxJ-rminJ)/n
R> r <- seq(0, rmax, by = rstep)    # r-distances for Lest
R> rJ <- seq(0, rmaxJ, by = rstepJ) # r-distances for Fest, Gest, Jest
\end{Sinput}
\end{Schunk}
%
Then perform simulations of CSR and calculate the $L$-functions saving the simulated patterns and functions:
%
%
\begin{Schunk}
\begin{Sinput}
R> env_L <- envelope(X, nsim = nsim,
+    simulate = expression(runifpoint(ex = X)),
+    fun = "Lest", correction = "translate",
+    transform = expression(.-r), # Take the L(r)-r function instead of L(r)
+    r = r,                       # Specify the distance vector
+    savefuns = TRUE,             # Save the estimated functions
+    savepatterns = TRUE,         # Save the simulated patterns
+    verbose = FALSE)
R> # The simulations can be obtained from the returned object:
R> simulations <- attr(env_L, "simpatterns")
\end{Sinput}
\end{Schunk}
%
And then the other test functions $F$, $G$, $J$ should be calculated for each simulated pattern:
%
\begin{Schunk}
\begin{Sinput}
R> env_F <- envelope(X, nsim = nsim,
+    simulate = simulations,
+    fun = "Fest", correction = "Kaplan", r = rJ,
+    savefuns = TRUE, verbose = FALSE)
R> env_G <- envelope(X, nsim = nsim,
+    simulate = simulations,
+    fun = "Gest", correction = "km", r = rJ,
+    savefuns = TRUE, verbose = FALSE)
R> env_J <- envelope(X, nsim = nsim,
+    simulate = simulations,
+    fun = "Jest", correction = "none", r = rJ,
+    savefuns = TRUE, verbose = FALSE)
\end{Sinput}
\end{Schunk}
%
All the curves can then be cropped to the desired r-interval $I$, if needed,
%
\begin{Schunk}
\begin{Sinput}
R> curve_set_L <- crop_curves(env_L, r_min = rmin, r_max = rmax)
R> curve_set_F <- crop_curves(env_F, r_min = rminJ, r_max = rmaxJ)
R> curve_set_G <- crop_curves(env_G, r_min = rminJ, r_max = rmaxJ)
R> curve_set_J <- crop_curves(env_J, r_min = rminJ, r_max = rmaxJ)
\end{Sinput}
\end{Schunk}
%
and finally the combined global envelope calculated and plotted
%
\begin{Schunk}
\begin{Sinput}
R> res_combined <- global_envelope_test(curve_sets = list(curve_set_L, curve_set_F,
+                                            curve_set_G, curve_set_J))
R> plot(res_combined, labels = c("L(r)-r", "F(r)", "G(r)", "J(r)"))
\end{Sinput}
\end{Schunk}
\includegraphics{pointpatterns-combined_test}
%


\section{The problem of NA values in the curve set}

With some functions in spatial statistics, missing NA values often appear for some distances $r$, which typically are not of interest, e.g., for the $J$-function with too large distances and for the pair-correlation function at zero. There is currently no automatic way to compute the envelopes with NA values existing in the \code{curve\_set} object: when such distances are present, they should be removed before computation of the envelope. This can be done using the function \fct{crop\_curves}. The function allows 1) to crop the curves to a user-specified interval $[r_{\min}, r_{\max}]$ (arguments \code{r_min} and \code{r_max}), or 2) to crop away all $r$-distances with NA (or infinite) values.

Consider the example of adult trees from \pkg{GET} \citep[][Section 3.2]{MyllymakiMrkvicka2023}
replacing the $L$-function with the $J$-function.
The \fct{envelope} function returns some NAs for the range of $r$-values specified below.
Without cropping away these, the function \fct{global\_envelope\_test} returns an error.
A working example is as follows (instead of \code{allfinite} one could specify cropping through \code{r_max}):
%
\begin{Schunk}
\begin{Sinput}
R> data("adult_trees")
R> X <- as.ppp(adult_trees, W = square(75))
R> env <- envelope(X, nsim = 999, fun = "Jest", correction = "km",
+    simulate = expression(runifpoint(ex = X)),
+    savefuns = TRUE, verbose = FALSE, r = seq(0, 10, length = 512))
R> cset <- crop_curves(env, allfinite=TRUE)
R> res <- global_envelope_test(cset)
\end{Sinput}
\end{Schunk}
%


\section{A one-stage goodness-of-fit test (typically conservative!)}

It is possible to perform a one-stage goodness-of-fit test for point process models as follows, accepting that the test may be conservative (or liberal). In literature, it has been recommended that the tests may be used if the test function is not closely related to the estimation procedure that was used to fit the model. However, \fct{GET.composite} can be used for adjusted tests, see the help file of this function and an example below.
%
\begin{Schunk}
\begin{Sinput}
R> X <- unmark(spruces)
R> # Minimum distance between points in the pattern
R> min(nndist(X))
\end{Sinput}
\begin{Soutput}
[1] 1.044031
\end{Soutput}
\begin{Sinput}
R> # Fit a model
R> fittedmodel <- ppm(X, interaction = Hardcore(hc = 1)) # Hardcore process
\end{Sinput}
\end{Schunk}
%

Simulating Gibbs process by \fct{envelope} is slow, because it uses an MCMC algorithm
%
\begin{Schunk}
\begin{Sinput}
R> #env <- envelope(fittedmodel, fun = "Jest", nsim = 999, savefuns = TRUE,
R> #                correction = "none", r = seq(0, 4, length = 500))
\end{Sinput}
\end{Schunk}
%
Using direct algorihm can be faster, because the perfect simulation is used here.
Therefore, in the following we utilize the function \fct{rHardcore}:
%
%
\begin{Schunk}
\begin{Sinput}
R> simulations <- NULL
R> nsim <- 999 # Number of simulations
R> for(j in 1:nsim) {
+     simulations[[j]] <- rHardcore(beta = exp(fittedmodel$coef[1]),
+                           R = fittedmodel$interaction$par$hc,
+                           W = X$window)
+  }
R> env_HC <- envelope(X, simulate = simulations, fun = "Jest",
+    nsim = length(simulations),
+    savefuns = TRUE, correction = "none",
+    r = seq(0, 4, length = 500),
+    verbose = FALSE)
R> curve_set <- crop_curves(env_HC, r_min = 1, r_max = 3.5)
R> res_HC <- global_envelope_test(curve_set, type = "erl")
R> plot(res_HC) + ylab(expression(italic(J(r))))
\end{Sinput}
\end{Schunk}
\includegraphics{pointpatterns-HC_goodness-of-fit_test}
%

Note: Conditioning the Gibbs hard-core model on the number of points, the only parameter in
the hard-core model is the hard-core distance. It is possible to fix the hard-core distance, e.g.
to the minimum distance between two points in the data, and then the test of
the hard-core model with fixed hard-core distance is simple (no parameters involved), and thus exact.
This example is given in \citet{MyllymakiEtal2017} and it is possible to prepare simulations for this case utilizing
the function \fct{rmh} of the \pkg{spatstat} package.
In the above example conditioning on the number of points was not employed.

\section{Adjusted global envelope test for composite hypotheses}

Let us test the fit of a Matern cluster process for the sapling data as an example of a composite hypothesis test.
The adjusted test of the \pkg{GET} package is described in Section 2.3 of \citet{MyllymakiMrkvicka2023}. The procedure was suggested by \citet{BaddeleyEtal2017} and extended for global envelopes in \citet{MyllymakiMrkvicka2023}.
%
\begin{Schunk}
\begin{Sinput}
R> data(saplings)
R> saplings <- as.ppp(saplings, W = square(75))
R> par(mfrow = c(1,1), mgp = c(0, 0, 0), mar = c(0, 0, 0, 0))
R> plot(saplings, main = "")
\end{Sinput}
\end{Schunk}
\includegraphics{pointpatterns-saplings}
%

First define the $r$-distances and number of simulations. Here we set the number of simulations just to 19, for fast exploration of the code, but for serious analysis we recommend at least 499 simulations.
%
\begin{Schunk}
\begin{Sinput}
R> rmin <- 0.3; rmax <- 10; rstep <- (rmax-rmin)/500
R> r <- seq(0, rmax, by = rstep)
R> nsim <- 19 # Increase nsim for serious analysis!
\end{Sinput}
\end{Schunk}
%
The Matern cluster process can be fitted to the pattern using the \fct{ppm} function of \pkg{spatstat}. This utilizes minimum contrast estimation with the $K$-function.
Then the adjusted global area rank envelope test can be performed using the function \fct{GET.composite}. Below we use the centred $L(r)$ function as the test function. The argument \code{type} specifies the global envelope test, see the help file of the \fct{global\_envelope\_test} function.
%
%
\begin{Schunk}
\begin{Sinput}
R> M1 <- kppm(saplings~1, clusters = "MatClust", statistic = "K")
R> adjenvL <- GET.composite(X = M1, nsim = nsim,
+    testfuns = list(L = list(fun="Lest", correction = "translate",
+         transform = expression(.-r), r = r)), # passed to envelope
+    type = "area", r_min = rmin, r_max = rmax, verbose = FALSE)
R> plot(adjenvL)
\end{Sinput}
\end{Schunk}
\includegraphics{pointpatterns-saplings_adjusted_test}
%

\section{Testing global and local dependence on covariates}

The function \fct{GET.spatialF} performs the one-stage global envelope tests based on spatial $F$- and $S$-statistics \citep{MyllymakiEtal2020} to explore the effects of covariates in parametric point process models.

Let us look at a simple example of tropical rain forest trees.
%
\begin{Schunk}
\begin{Sinput}
R> data(bei)
\end{Sinput}
\end{Schunk}
%
Let us study the effect of gradient on the intensity of the trees. We define the full model including this interesting covariate and the reduced model, which is otherwise the same as the full model, but the interesting covariate is excluded. Further the function \fct{fitppm} defines how the (full or reduced) model can be fitted to the point pattern.
%
%
\begin{Schunk}
\begin{Sinput}
R> fullmodel <- ~ grad
R> reducedmodel <- ~ 1
R> fitppm <- function(X, model, covariates) {
+    ppm(X, model, covariates = covariates)
+  }
R> nsim <- 19 # Increase nsim for serious analysis!
R> res_sF <- GET.spatialF(bei, fullmodel, reducedmodel,
+                      fitppm, bei.extra, nsim)
\end{Sinput}
\begin{Soutput}
Generating 19 simulated patterns ...1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17,
18, 
19.
\end{Soutput}
\end{Schunk}
%
\begin{Schunk}
\begin{Sinput}
R> plot(res_sF$F, what = c("obs", "hi", "hi.sign"), sign.type = "col")
\end{Sinput}
\end{Schunk}
\includegraphics{pointpatterns-spatialF_bei_F}
%
\begin{Schunk}
\begin{Sinput}
R> plot(res_sF$S, what = c("obs", "hi", "hi.sign"), sign.type = "col")
\end{Sinput}
\end{Schunk}
\includegraphics{pointpatterns-spatialF_bei_S}
%

Another example is the test of the effect of elevation on the point pattern of lightnings.
%
%
\begin{Schunk}
\begin{Sinput}
R> # Example of forest fires
R> data("clmfires")
R> # Choose the locations of the lightnings in years 2004-2007:
R> pp.lightning <- unmark(subset(clmfires, cause == "lightning" &
+                   date >= "2004-01-01" & date < "2008-01-01"))
R> covariates <- clmfires.extra$clmcov100
R> covariates$forest <- 
+    covariates$landuse == "conifer" | covariates$landuse == "denseforest" |
+    covariates$landuse == "mixedforest"
R> fullmodel <- ~ elevation + landuse
R> reducedmodel <- ~ landuse
R> nsim <- 19 # Increase nsim for serious analysis!
R> res_sF2 <- GET.spatialF(pp.lightning, fullmodel, reducedmodel,
+                      fitppm, covariates, nsim)
\end{Sinput}
\begin{Soutput}
Generating 19 simulated patterns ...1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17,
18, 
19.
\end{Soutput}
\end{Schunk}
%
\begin{Schunk}
\begin{Sinput}
R> plot(res_sF2$F, what = c("obs", "hi", "hi.sign"), sign.type = "col")
\end{Sinput}
\end{Schunk}
\includegraphics{pointpatterns-spatialF_forestfires_F}
%
\begin{Schunk}
\begin{Sinput}
R> plot(res_sF2$S, what = c("obs", "hi", "hi.sign"), sign.type = "col")
\end{Sinput}
\end{Schunk}
\includegraphics{pointpatterns-spatialF_forestfires_S}
%

Above only the inhomogeneous Poisson process was used as the model. Examples of the \fct{fitfun} functions for clustered and regular processes are:
%
\begin{Schunk}
\begin{Sinput}
R> # fitfun for the log Gaussian Cox Process with exponential covariance function
R> fitLGCPexp <- function(X, model, covariates) {
+    kppm(X, model, clusters = "LGCP", model = "exponential", covariates = covariates)
+  }
R> # fitfun for the hardcore process with hardcore radius 0.01
R> fitHardcore <- function(X, model, covariates) {
+    ppm(X, model, interaction = Hardcore(0.01), covariates = covariates)
+  }
\end{Sinput}
\end{Schunk}
%

\section{An example analysis of the saplings data set}

This is the example of \citet[][Supplement S4]{MyllymakiEtal2017}.

The saplings data set is available at the \pkg{GET} package.
%
\begin{Schunk}
\begin{Sinput}
R> data(saplings)
R> saplings <- as.ppp(saplings, W = square(75))
\end{Sinput}
\end{Schunk}
%

First choose the r-distances for $L(r)$ and $J(r)$ functions, respectively.
%
\begin{Schunk}
\begin{Sinput}
R> nr <- 500
R> rmin <- 0.3; rminJ <- 0.3
R> rmax <- 10; rmaxJ <- 6
R> rstep <- (rmax-rmin)/nr; rstepJ <- (rmaxJ-rminJ)/nr
R> r <- seq(0, rmax, by = rstep)
R> rJ <- seq(0, rmaxJ, by = rstepJ)
\end{Sinput}
\end{Schunk}
%

\subsection{The CSR test based on the L(r)-r function}

Note: CSR is simulated by fixing the number of points and generating \code{nsim} simulations
from the binomial process, i.e. we deal with a simple hypothesis.

First, the envelope function of the spatstat package can be used to generate nsim simulations under CSR and
to calculate the centred $L$-function for the data and each simulation.
%
%
\begin{Schunk}
\begin{Sinput}
R> nsim <- 1999 # Number of simulations
R> env <- envelope(saplings, nsim = nsim,
+    simulate = expression(runifpoint(ex = saplings)), # Simulate CSR
+    fun = "Lest", correction = "translate", # estimator of L with transl. edge corr.
+    transform = expression(.-r), # Take the L(r)-r function instead of L(r)
+    r = r, # Specify the distance vector
+    savefuns = TRUE, # Save the estimated functions
+    verbose = FALSE)
\end{Sinput}
\end{Schunk}
%
Then the curves can be cropped to the desired interval of distances $[r_{\min}, r_{\max}]$
%
\begin{Schunk}
\begin{Sinput}
R> curve_set <- crop_curves(env, r_min = rmin, r_max = rmax)
\end{Sinput}
\end{Schunk}
%
And a global envelope test done by means of the \fct{global\_envelope\_test} function (\code{type="rank"} and larger \code{nsim} was used in \citet[][S4]{MyllymakiEtal2017}:
%
\begin{Schunk}
\begin{Sinput}
R> res_sapl <- global_envelope_test(curve_set, type = "erl")
R> plot(res_sapl) + ylab(expression(italic(hat(L)(r)-r)))
\end{Sinput}
\end{Schunk}
\includegraphics{pointpatterns-saplings_csr_L3}
%

The CSR hypothesis is clearly rejected and the rank envelope indicates clear
clustering of saplings. As a next step, we explore the Matern cluster process as a null model.
This is a composite hypothesis.

\subsection{Testing the fit of a Matern cluster process}

First we fit the Matern cluster process to the pattern. Here we use the minimum contrast estimation with the $K$-funcion (the pair correction function can be chosen by setting \code{statistic = "pcf"}).
%
\begin{Schunk}
\begin{Sinput}
R> fitted_model <- kppm(saplings~1, clusters = "MatClust", statistic = "K")
\end{Sinput}
\end{Schunk}
%
Next step is to perform the adjusted directional quantile global envelope test using the centred $L$-function.
(For the rank envelope test, choose \code{type = "rank"} instead and increase \code{nsim}.)
%
%
\begin{Schunk}
\begin{Sinput}
R> nsim <- 19 # 19 just for experimenting with the code!!
R> #nsim <- 499 # 499 is ok for type = 'qdir' (takes > 1 h)
R> adjenvL_sapl <- GET.composite(X = fitted_model,
+    fun = "Lest", correction = "translate",
+    transform = expression(.-r), r = r,
+    type = "qdir", nsim = nsim, nsimsub = nsim,
+    r_min = rmin, r_max = rmax, verbose = FALSE)
\end{Sinput}
\end{Schunk}
%
The result can then be plotted:
\begin{Schunk}
\begin{Sinput}
R> plot(adjenvL_sapl) + ylab(expression(italic(L(r)-r)))
\end{Sinput}
\end{Schunk}
\includegraphics{pointpatterns-saplings_matclust_L_plot}
%

From the test with the centred $L$-function, it appears that the Matern cluster model would be
a reasonable model for the saplings pattern.
To further explore the goodness-of-fit of the Matern cluster process, test the
model with the $J$-function:
This takes quite some time if \code{nsim} is reasonably large.
%
%
\begin{Schunk}
\begin{Sinput}
R> adjenvJ_sapl <- GET.composite(X = fitted_model,
+    fun = "Jest", correction = "none", r = rJ,
+    type = "qdir", nsim = nsim, nsimsub = nsim,
+    r_min = rminJ, r_max = rmaxJ, verbose = FALSE)
\end{Sinput}
\end{Schunk}
%
And, plot the result
%
\begin{Schunk}
\begin{Sinput}
R> plot(adjenvJ_sapl) + ylab(expression(italic(J(r))))
\end{Sinput}
\end{Schunk}
\includegraphics{pointpatterns-saplings_matclust_J_plot}
%

Thus, it appears that the Matern cluster process is not adequate for the saplings data.
The Matern cluster model might be interpreted as a regeneration process in circular gaps
between large trees. However, it is possible that the gap openings in the forest were not exactly
circular, thereby leading to the rejection of the model by the $J$-function.

It is also possible to test the fit of the Matern cluster process simultaneously by the two test functions:
%
%
\begin{Schunk}
\begin{Sinput}
R> adjenvLJ_sapl <- GET.composite(X = fitted_model,
+    testfuns = list(L = list(fun = "Lest", correction = "translate",
+                        transform = expression(.-r), r = r),
+                    J = list(fun = "Jest", correction = "none", r = rJ)),
+    type = "erl", nsim = nsim, nsimsub = nsim,
+    r_min = c(rmin, rminJ), r_max = c(rmax, rmaxJ),
+    save.cons.envelope = TRUE, verbose = FALSE)
R> plot(adjenvLJ_sapl)
\end{Sinput}
\end{Schunk}
\includegraphics{pointpatterns-saplings_matclust4}
%

%% -- Bibliography -------------------------------------------------------------

\bibliography{GETbibfile}


\end{document}