Description Usage Arguments Details Value References See Also Examples
The rank envelope test, p-value and global envelope
1 2 | rank_envelope(curve_set, alpha = 0.05, savedevs = FALSE,
alternative = c("two.sided", "less", "greater"), lexo = FALSE, ties)
|
curve_set |
A curve_set (see |
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 |
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.
An "envelope_test" object containing the following fields:
r = Distances for which the test was made.
method = The name of the envelope test.
alternative = The alternative specified in the function call.
p = A point estimate for the p-value (default is the mid-rank p-value).
ties = As the argument ties
.
p_interval = The p-value interval [p_liberal, p_conservative].
k_alpha = The value of k corresponding to the 100(1-alpha)% global envelope.
k = Global rank values (k[1] is the value for the data pattern). Returned only if savedevs = TRUE.
central_curve = If the curve_set (or envelope object) contains a component 'theo', then this function is used as the central curve and returned in this component. Otherwise, the central_curve is the mean of the test functions T_i(r), i=2, ..., s+1. Used for visualization only.
data_curve = The test function for the data.
lower = The lower envelope.
upper = The upper envelope.
call = The call of the function.
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]
random_labelling
, plot.envelope_test
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)
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.