quadrat.test: Dispersion Test for Spatial Point Pattern Based on Quadrat Counts

Description

Performs a test of Complete Spatial Randomness for a given point pattern, based on quadrat counts. Alternatively performs a goodness-of-fit test of a fitted inhomogeneous Poisson model. By default performs chi-squared tests; can also perform Monte Carlo based tests.

Usage

quadrat.test(X, ...)
# S3 method for ppp
quadrat.test(X, nx=5, ny=nx,
                          alternative=c("two.sided", "regular", "clustered"),
                           method=c("Chisq", "MonteCarlo"),
                           conditional=TRUE, CR=1,
                           lambda=NULL, df.est=NULL,
                           ...,
                           xbreaks=NULL, ybreaks=NULL, tess=NULL,
                           nsim=1999)
# S3 method for ppm
quadrat.test(X, nx=5, ny=nx, 
                          alternative=c("two.sided", "regular", "clustered"),
                           method=c("Chisq", "MonteCarlo"),
                           conditional=TRUE, CR=1, df.est=NULL,
                           ..., 
                           xbreaks=NULL, ybreaks=NULL, tess=NULL,
                           nsim=1999)
# S3 method for quadratcount
quadrat.test(X,
                          alternative=c("two.sided", "regular", "clustered"),
                          method=c("Chisq", "MonteCarlo"),
                          conditional=TRUE, CR=1,
                          lambda=NULL, df.est=NULL,
                          ...,
                          nsim=1999)

Arguments

A point pattern (object of class "ppp") to be subjected to the goodness-of-fit test. Alternatively a fitted point process model (object of class "ppm") to be tested. Alternatively X can be the result of applying quadratcount to a point pattern.

nx,ny

Numbers of quadrats in the $x$ and $y$ directions. Incompatible with xbreaks and ybreaks.

alternative

Character string (partially matched) specifying the alternative hypothesis.

method

Character string (partially matched) specifying the test to use: either method="Chisq" for the chi-squared test (the default), or method="MonteCarlo" for a Monte Carlo test.

conditional

Logical. Should the Monte Carlo test be conducted conditionally upon the observed number of points of the pattern? Ignored if method="Chisq".

Optional. Numerical value. The exponent for the Cressie-Read test statistic. See Details.

lambda

Optional. Pixel image (object of class "im") or function (class "funxy") giving the predicted intensity of the point process.

df.est

Optional. Advanced use only. The number of fitted parameters, or the number of degrees of freedom lost by estimation of parameters.

…

Ignored.

xbreaks

Optional. Numeric vector giving the $x$ coordinates of the boundaries of the quadrats. Incompatible with nx.

ybreaks

Optional. Numeric vector giving the $y$ coordinates of the boundaries of the quadrats. Incompatible with ny.

tess

Tessellation (object of class "tess" or something acceptable to as.tess) determining the quadrats. Incompatible with nx, ny, xbreaks, ybreaks.

nsim

The number of simulated samples to generate when method="MonteCarlo".

Value

An object of class "htest". See chisq.test for explanation.

The return value is also an object of the special class "quadrattest", and there is a plot method for this class. See the examples.

Details

These functions perform $χ^{2}$ tests or Monte Carlo tests of goodness-of-fit for a point process model, based on quadrat counts.

The function quadrat.test is generic, with methods for point patterns (class "ppp"), split point patterns (class "splitppp"), point process models (class "ppm") and quadrat count tables (class "quadratcount").

if X is a point pattern, we test the null hypothesis that the data pattern is a realisation of Complete Spatial Randomness (the uniform Poisson point process). Marks in the point pattern are ignored. (If lambda is given then the null hypothesis is the Poisson process with intensity lambda.)
if X is a split point pattern, then for each of the component point patterns (taken separately) we test the null hypotheses of Complete Spatial Randomness. See quadrat.test.splitppp for documentation.
If X is a fitted point process model, then it should be a Poisson point process model. The data to which this model was fitted are extracted from the model object, and are treated as the data point pattern for the test. We test the null hypothesis that the data pattern is a realisation of the (inhomogeneous) Poisson point process specified by X.

In all cases, the window of observation is divided into tiles, and the number of data points in each tile is counted, as described in quadratcount. The quadrats are rectangular by default, or may be regions of arbitrary shape specified by the argument tess. The expected number of points in each quadrat is also calculated, as determined by CSR (in the first case) or by the fitted model (in the second case). Then the Pearson $X^{2}$ statistic $X^{2} = s u m ((o b s e r v e d - e x p e c t e d)^{2} / e x p e c t e d)$ is computed.

If method="Chisq" then a $χ^{2}$ test of goodness-of-fit is performed by comparing the test statistic to the $χ^{2}$ distribution with $m - k$ degrees of freedom, where m is the number of quadrats and $k$ is the number of fitted parameters (equal to 1 for quadrat.test.ppp). The default is to compute the two-sided $p$ -value, so that the test will be declared significant if $X^{2}$ is either very large or very small. One-sided $p$ -values can be obtained by specifying the alternative. An important requirement of the $χ^{2}$ test is that the expected counts in each quadrat be greater than 5.

If method="MonteCarlo" then a Monte Carlo test is performed, obviating the need for all expected counts to be at least 5. In the Monte Carlo test, nsim random point patterns are generated from the null hypothesis (either CSR or the fitted point process model). The Pearson $X^{2}$ statistic is computed as above. The $p$ -value is determined by comparing the $X^{2}$ statistic for the observed point pattern, with the values obtained from the simulations. Again the default is to compute the two-sided $p$ -value.

If conditional is TRUE then the simulated samples are generated from the multinomial distribution with the number of “trials” equal to the number of observed points and the vector of probabilities equal to the expected counts divided by the sum of the expected counts. Otherwise the simulated samples are independent Poisson counts, with means equal to the expected counts.

If the argument CR is given, then instead of the Pearson $X^{2}$ statistic, the Cressie-Read (1984) power divergence test statistic $2 n I = \frac{2}{C R (C R + 1)} \sum_{i} [{(\frac{X_{i}}{E_{i}})}^{C} R - 1]$ is computed, where $X_{i}$ is the $i$ th observed count and $E_{i}$ is the corresponding expected count. The value CR=1 gives the Pearson $X^{2}$ statistic; CR=0 gives the likelihood ratio test statistic $G^{2}$ ; CR=-1/2 gives the Freeman-Tukey statistic $T^{2}$ ; CR=-1 gives the modified likelihood ratio test statistic $G M^{2}$ ; and CR=-2 gives Neyman's modified statistic $N M^{2}$ . In all cases the asymptotic distribution of this test statistic is the same $χ^{2}$ distribution as above.

The return value is an object of class "htest". Printing the object gives comprehensible output about the outcome of the test.

The return value also belongs to the special class "quadrat.test". Plotting the object will display the quadrats, annotated by their observed and expected counts and the Pearson residuals. See the examples.

References

Cressie, N. and Read, T.R.C. (1984) Multinomial goodness-of-fit tests. Journal of the Royal Statistical Society, Series B 46, 440--464.

Examples

Run this code

# NOT RUN {
  quadrat.test(simdat)
  quadrat.test(simdat, 4, 3)

  quadrat.test(simdat, alternative="regular")
  quadrat.test(simdat, alternative="clustered")

  ## Likelihood ratio test
  quadrat.test(simdat, CR=0)
  ## Power divergence tests
  quadrat.test(simdat, CR=-1)$p.value
  quadrat.test(simdat, CR=-2)$p.value

  # Using Monte Carlo p-values
  quadrat.test(swedishpines) # Get warning, small expected values.
  # quadrat.test(swedishpines, method="M", nsim=4999)
  #  quadrat.test(swedishpines, method="M", nsim=4999, conditional=FALSE)

  
# }
# NOT RUN {
  # quadrat counts
  qS <- quadratcount(simdat, 4, 3)
  quadrat.test(qS)

  # fitted model: inhomogeneous Poisson
  fitx <- ppm(simdat ~ x)
  quadrat.test(fitx)

  # an equivalent test (results differ due to discretisation effects):
  quadrat.test(simdat, lambda=predict(fitx), df.est=length(coef(fitx)))

  te <- quadrat.test(simdat, 4)
  residuals(te)  # Pearson residuals

  plot(te)

  plot(simdat, pch="+", cols="green", lwd=2)
  plot(te, add=TRUE, col="red", cex=1.4, lty=2, lwd=3)

  sublab <- eval(substitute(expression(p[chi^2]==z),
                       list(z=signif(te$p.value,3))))
  title(sub=sublab, cex.sub=3)

  # quadrats of irregular shape
  B <- dirichlet(runifpoint(6, Window(simdat)))
  qB <- quadrat.test(simdat, tess=B)
  plot(simdat, main="quadrat.test(simdat, tess=B)", pch="+")
  plot(qB, add=TRUE, col="red", lwd=2, cex=1.2)

# }

Run the code above in your browser using DataLab

Get 50% off unlimited learning