rank_envelope: The rank envelope test

Description Usage Arguments Details Value References See Also Examples

Description

The rank envelope test, p-value and global envelope

Usage

1
2
rank_envelope(curve_set, alpha = 0.05, savedevs = FALSE,
  alternative = c("two.sided", "less", "greater"), lexo = FALSE, ties)

Arguments

curve_set

A curve_set (see create_curve_set) or an envelope object. If an envelope object is given, it must contain the summary functions from the simulated patterns which can be achieved by setting savefuns = TRUE when calling envelope().

alpha

The significance level. The 100(1-alpha)% global envelope will be calculated.

savedevs

Logical. Should the global rank values k_i, i=1,...,nsim+1 be returned? Default: FALSE.

alternative

A character string specifying the alternative hypothesis. Must be one of the following: "two.sided" (default), "less" or "greater".

lexo

Logical, whether or not to use rank count ordering for calculation of the p-value. See details.

ties

Ties method to be passed to estimate_p_value. Used to obtain a point estimate for the p-value. The default point estimate is the mid-rank p-value.

Details

The rank envelope test is a completely non-parametric test, which provides a p-value interval given by the most liberal and the most conservative p-value estimate and the 100(1-alpha)% global envelope for the chosen test function T(r) on the chosen interval of distances.

Given a curve_set (or an envelope) object, the test is carried out as follows.

For each curve in the curve_set, both the data curve and the simulations, the global rank measure R is determined. If savedevs = TRUE, then the global rank values R_1, R_2, ..., R_(s+1) are returned in the component 'k', where k[1] is the value for the data.

Based on R_i, i=1, ..., s+1, the p-interval is calculated. This interval is by default plotted for the object returned by the rank_envelope function. Also a single p-value is calculated and returned in component 'p'. By default this p-value is the mid-rank p-value, but another option can be used by specifying ties argument which is passed to estimate_p_value. For options see estimate_p_value.

The 100(1-alpha)% global envelope is given by the 'k_alpha'th lower and upper envelope. For details see Myllymäki et al. (2013).

The above holds for p-value calculation if lexo == FALSE and then the test corresponds to the rank envelope test by Myllymaki et. al (2013). If lexo == TRUE, then all the pointwise ranks are used to rank the curves by rank count ordering (Myllymäki et al., 2015) and the single p-value in p is the p-value based on the rank count ordering.

The rank count ordering test allows in principle a lower number of simulations to be used, but then the test may no longer be usable as a graphical test.

Value

An "envelope_test" object containing the following fields:

References

Myllymäki, M., Mrkvička, T., Seijo, H., Grabarnik, P. (2013). Global envelope tests for spatial point patterns. arXiv:1307.0239 [stat.ME]

Myllymäki, M., Mrkvička, T., Grabarnik, P., Seijo, H. and Hahn, U. (2015). Global envelope tests for spatial point patterns. arXiv:1307.0239v4 [stat.ME]

See Also

random_labelling, plot.envelope_test

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
## Testing complete spatial randomness (CSR)
#-------------------------------------------
require(spatstat)
pp <- unmark(spruces)
# Generate nsim simulations under CSR, calculate L-function for the data and simulations
env <- envelope(pp, fun="Lest", nsim=2499, savefuns=TRUE, correction="translate")
# The rank envelope test
res <- rank_envelope(env)
# Plot the result.
# - The central curve is now obtained from env[['theo']], which is the
# value of the L-function under the null hypothesis (L(r) = r).
plot(res)
# or (requires R library ggplot2)
plot(res, use_ggplot2=TRUE)

## Advanced use:
# Choose the interval of distances [r_min, r_max] (at the same time create a curve_set from 'env')
curve_set <- crop_curves(env, r_min = 1, r_max = 7)
# For better visualisation, take the L(r)-r function
curve_set <- residual(curve_set, use_theo = TRUE)
# Do the rank envelope test
res <- rank_envelope(curve_set); plot(res, use_ggplot2=TRUE)

## Random labeling test
#----------------------
# requires library 'marksummary'
mpp <- spruces
# 1) Perform simulations under the random labelling hypothesis and calculate
# the test function T(r) for the data pattern (mpp) and each simulation.
# The command below specifies that the test function is T(r) = \hat{L}_m(r),
# which is an estimator of the mark-weighted L function, L_m(r),
# with translational edge correction (default).
# The random_labelling function returns the centred functions \hat{L}_m(r)-T_0(r),
# where T_0(r) = \hat{L}(r) is the unmarked L function.
curve_set <- random_labelling(mpp, mtf_name = 'm', nsim=2499, r_min=1.5, r_max=9.5)
# 2) Do the rank envelope test
res <- rank_envelope(curve_set)
# 3) Plot the test result
plot(res, use_ggplot2=TRUE, ylab=expression(italic(L[m](r)-L(r))))

# Make the test using instead the test function T(r) = \hat{L}_mm(r);
# which is an estimator of the mark-weighted L function, L_mm(r),
# with translational edge correction (default).
curve_set <- random_labelling(mpp, mtf_name = 'mm', nsim=2499, r_min=1.5, r_max=9.5)
res <- rank_envelope(curve_set)
plot(res, use_ggplot2=TRUE, ylab=expression(italic(L[mm](r)-L(r))))

## Goodness-of-fit test (typically conservative)
#-----------------------------------------------
pp <- unmark(spruces)
# Minimum distance between points in the pattern
min(nndist(pp))
# Fit a model
fittedmodel <- ppm(pp, interaction=Hardcore(hc=1)) # Hardcore process

## Not run: 
# Simulating Gibbs process by 'envelope' is slow, because it uses the MCMC algorithm
#env <- envelope(fittedmodel, fun="Jest", nsim=999, savefuns=TRUE,
                 correction="none", r=seq(0, 4, length=500))

# Using direct algorihm can be faster, because the perfect simulation is used here.
simulations <- NULL
for(j in 1:2499) {
   simulations[[j]] <- rHardcore(beta=exp(fittedmodel$coef[1]),
                                 R = fittedmodel$interaction$par$hc,
                                 W = pp$window);
   if(j%%10==0) cat(j, "...", sep="")
}
env <- envelope(pp, simulate=simulations, fun="Jest", nsim=length(simulations),
                savefuns=TRUE, correction="none", r=seq(0, 4, length=500))
curve_set <- crop_curves(env, r_min = 1, r_max = 3.5)
res <- rank_envelope(curve_set); plot(res, use_ggplot2=TRUE)

## End(Not run)

myllym/spptest documentation built on May 23, 2019, noon