envelope.R
#
# envelope.R
#
# computes simulation envelopes
#
# $Revision: 2.4 $ $Date: 2010/03/18 14:37:02 $
#
envelope <- function(Y, fun, ...) {
UseMethod("envelope")
}
# .................................................................
# A 'simulation recipe' contains the following variables
#
# type = Type of simulation
# "csr": uniform Poisson process
# "rmh": simulated realisation of fitted Gibbs model
# "kppm": simulated realisation of fitted cluster model
# "expr": result of evaluating a user-supplied expression
# "list": user-supplied list of point patterns
#
# expr = expression that is repeatedly evaluated to generate simulations
#
# envir = environment in which to evaluate the expression `expr'
#
# 'csr' = TRUE iff the model is (known to be) uniform Poisson
#
# ...................................................................
simulrecipe <- function(type, expr, envir, csr) {
out <- list(type=type,
expr=expr,
envir=envir,
csr=csr)
class(out) <- "simulrecipe"
out
}
envelope.ppp <-
function(Y, fun=Kest, nsim=99, nrank=1, ...,
simulate=NULL, verbose=TRUE, clipdata=TRUE,
transform=NULL, global=FALSE, ginterval=NULL,
savefuns=FALSE, savepatterns=FALSE, nsim2=nsim,
VARIANCE=FALSE, nSD=2,
Yname=NULL) {
cl <- match.call()
if(is.null(Yname)) Yname <- deparse(substitute(Y))
envir.user <- parent.frame()
envir.here <- sys.frame(sys.nframe())
if(is.null(simulate)) {
# ...................................................
# Realisations of complete spatial randomness
# will be generated by rpoispp
# Data pattern X is argument Y
# Data pattern determines intensity of Poisson process
X <- Y
sY <- summary(Y, checkdup=FALSE)
Yintens <- sY$intensity
Ywin <- Y$window
# expression that will be evaluated
simexpr <-
if(!is.marked(Y)) {
# unmarked point pattern
expression(rpoispp(Yintens, win=Ywin))
} else {
# marked point pattern
Ymarx <- marks(Y, dfok=FALSE)
expression({A <- rpoispp(Yintens, win=Ywin);
A %mark% sample(Ymarx, A$n, replace=TRUE)})
}
# evaluate in THIS environment
simrecipe <- simulrecipe(type = "csr",
expr = simexpr,
envir = envir.here,
csr = TRUE)
} else {
# ...................................................
# Simulations are determined by 'simulate' argument
# Processing is deferred to envelopeEngine
simrecipe <- simulate
# Data pattern is argument Y
X <- Y
}
envelopeEngine(X=X, fun=fun, simul=simrecipe,
nsim=nsim, nrank=nrank, ...,
verbose=verbose, clipdata=clipdata,
transform=transform, global=global, ginterval=ginterval,
savefuns=savefuns, savepatterns=savepatterns, nsim2=nsim2,
VARIANCE=VARIANCE, nSD=nSD,
Yname=Yname, cl=cl,
envir.user=envir.user)
}
envelope.ppm <-
function(Y, fun=Kest, nsim=99, nrank=1, ...,
simulate=NULL, verbose=TRUE, clipdata=TRUE,
start=NULL, control=list(nrep=1e5, expand=1.5),
transform=NULL, global=FALSE, ginterval=NULL,
savefuns=FALSE, savepatterns=FALSE, nsim2=nsim,
VARIANCE=FALSE, nSD=2,
Yname=NULL) {
cl <- match.call()
if(is.null(Yname)) Yname <- deparse(substitute(Y))
envir.user <- parent.frame()
envir.here <- sys.frame(sys.nframe())
# Extract data pattern X from fitted model Y
X <- data.ppm(Y)
if(is.null(simulate)) {
# ...................................................
# Simulated realisations of the fitted model Y
# will be generated using rmh
# Set up parameters for rmh
rmodel <- rmhmodel(Y, verbose=FALSE)
if(is.null(start))
start <- list(n.start=X$n)
rstart <- rmhstart(start)
rcontr <- rmhcontrol(control)
# pre-digest arguments
rmhinfolist <- rmh(rmodel, rstart, rcontr, preponly=TRUE, verbose=FALSE)
# expression that will be evaluated
simexpr <- expression(rmhEngine(rmhinfolist, verbose=FALSE))
envir <- envir.here
# evaluate in THIS environment
simrecipe <- simulrecipe(type = "rmh",
expr = simexpr,
envir = envir.here,
csr = FALSE)
} else {
# ...................................................
# Simulations are determined by 'simulate' argument
# Processing is deferred to envelopeEngine
simrecipe <- simulate
}
envelopeEngine(X=X, fun=fun, simul=simrecipe,
nsim=nsim, nrank=nrank, ...,
verbose=verbose, clipdata=clipdata,
transform=transform, global=global, ginterval=ginterval,
savefuns=savefuns, savepatterns=savepatterns, nsim2=nsim2,
VARIANCE=VARIANCE, nSD=nSD,
Yname=Yname, cl=cl,
envir.user=envir.user)
}
envelope.kppm <-
function(Y, fun=Kest, nsim=99, nrank=1, ...,
simulate=NULL, verbose=TRUE, clipdata=TRUE,
transform=NULL, global=FALSE, ginterval=NULL,
savefuns=FALSE, savepatterns=FALSE, nsim2=nsim,
VARIANCE=FALSE, nSD=2, Yname=NULL)
{
cl <- match.call()
if(is.null(Yname)) Yname <- deparse(substitute(Y))
envir.user <- parent.frame()
envir.here <- sys.frame(sys.nframe())
# Extract data pattern X from fitted model Y
X <- Y$X
if(is.null(simulate)) {
# Simulated realisations of the fitted model Y
# will be generated using simulate.kppm
kmodel <- Y
# expression that will be evaluated
simexpr <- expression(simulate(kmodel)[[1]])
# evaluate in THIS environment
simrecipe <- simulrecipe(type = "kppm",
expr = simexpr,
envir = envir.here,
csr = FALSE)
} else {
# ...................................................
# Simulations are determined by 'simulate' argument
# Processing is deferred to envelopeEngine
simrecipe <- simulate
}
envelopeEngine(X=X, fun=fun, simul=simrecipe,
nsim=nsim, nrank=nrank, ...,
verbose=verbose, clipdata=clipdata,
transform=transform, global=global, ginterval=ginterval,
savefuns=savefuns, savepatterns=savepatterns, nsim2=nsim2,
VARIANCE=VARIANCE, nSD=nSD,
Yname=Yname, cl=cl,
envir.user=envir.user)
}
## .................................................................
## Engine for simulating and computing envelopes
## .................................................................
#
# X is the data point pattern, which could be ppp, pp3, ppx etc
# X determines the class of pattern expected from the simulations
#
envelopeEngine <-
function(X, fun, simul,
nsim=99, nrank=1, ...,
verbose=TRUE, clipdata=TRUE,
transform=NULL, global=FALSE, ginterval=NULL,
savefuns=FALSE, savepatterns=FALSE, nsim2=nsim,
VARIANCE=FALSE, nSD=2,
Yname=NULL, internal=NULL, cl=NULL,
envir.user=envir.user) {
#
envir.here <- sys.frame(sys.nframe())
# Identify class of patterns to be simulated, from data pattern X
Xclass <- if(is.ppp(X)) "ppp" else
if(is.pp3(X)) "pp3" else
if(is.ppx(X)) "ppx" else
stop("Unrecognised class of point pattern")
Xobjectname <- paste("point pattern of class", sQuote(Xclass))
# Identify type of simulation from argument 'simul'
if(inherits(simul, "simulrecipe")) {
# ..................................................
# simulation recipe is given
simtype <- simul$type
simexpr <- simul$expr
envir <- simul$envir
csr <- simul$csr
} else {
# ...................................................
# simulation is specified by argument `simulate' to envelope()
simulate <- simul
# which should be an expression, or a list of point patterns,
# or an envelope object.
csr <- FALSE
# override
if(!is.null(icsr <- internal$csr)) csr <- icsr
model <- NULL
if(inherits(simulate, "envelope")) {
# envelope object: see if it contains stored point patterns
simpat <- attr(simulate, "simpatterns")
if(!is.null(simpat))
simulate <- simpat
else
stop(paste("The argument", sQuote("simulate"),
"is an envelope object but does not contain",
"any saved point patterns."))
}
if(is.expression(simulate)) {
# The user-supplied expression 'simulate' will be evaluated repeatedly
simtype <- "expr"
simexpr <- simulate
envir <- envir.user
} else if(is.list(simulate) &&
( (is.ppp(X) && all(unlist(lapply(simulate, is.ppp))))
|| (is.pp3(X) && all(unlist(lapply(simulate, is.pp3))))
|| (is.ppx(X) && all(unlist(lapply(simulate, is.ppx)))))) {
# The user-supplied list of point patterns will be used
simtype <- "list"
SimDataList <- simulate
# expression that will be evaluated
simexpr <- expression(SimDataList[[i]])
envir <- envir.here
# ensure that `i' is defined
i <- 1
# any messages?
if(!is.null(mess <- attr(simulate, "internal"))) {
# determine whether these point patterns are realisations of CSR
csr <- !is.null(mc <- mess$csr) && mc
}
} else stop(paste(sQuote("simulate"),
"should be an expression, or a list of point patterns"))
}
#--------------------------------------------------------------------
# Determine number of simulations
#
#
## determine whether dual simulations are required
## (one set of simulations to calculate the theoretical mean,
## another independent set of simulations to obtain the critical point.)
dual <- (global && !csr && !VARIANCE)
Nsim <- if(!dual) nsim else (nsim + nsim2)
# if taking data from a list of point patterns,
# check there are enough of them
if(simtype == "list" && Nsim > length(SimDataList))
stop(paste("Number of simulations",
paren(if(!dual)
paste(nsim) else
paste(nsim, "+", nsim2, "=", Nsim)
),
"exceeds number of point pattern datasets supplied"))
# Undocumented secret exit
# ------------------------------------------------------------------
if(!is.null(eject <- internal$eject) && eject == "patterns") {
# generate simulated realisations and return only these patterns
if(verbose) {
action <- if(simtype == "list") "Extracting" else "Generating"
descrip <- switch(simtype,
csr = "simulations of CSR",
rmh = "simulated realisations of fitted Gibbs model",
kppm = "simulated realisations of fitted cluster model",
expr = "simulations by evaluating expression",
list = "point patterns from list",
"simulated realisations")
cat(paste(action, Nsim, descrip, "...\n"))
}
XsimList <- list()
# start simulation loop
for(i in 1:Nsim) {
if(verbose) progressreport(i, Nsim)
Xsim <- eval(simexpr, envir=envir)
if(!inherits(Xsim, Xclass))
switch(simtype,
csr={
stop(paste("Internal error:", Xobjectname, "not generated"))
},
rmh={
stop(paste("Internal error: rmh did not return an",
Xobjectname))
},
kppm={
stop(paste("Internal error: simulate.kppm did not return an",
Xobjectname))
},
expr={
stop(paste("Evaluating the expression", sQuote("simulate"),
"did not yield an", Xobjectname))
},
list={
stop(paste("Internal error: list entry was not an",
Xobjectname))
},
stop(paste("Internal error:", Xobjectname, "not generated"))
)
XsimList[[i]] <- Xsim
}
if(verbose) {
cat(paste("Done.\n"))
flush.console()
}
attr(XsimList, "internal") <- list(csr=csr)
return(XsimList)
}
# ------------------------------------------------------------------
# Summary function to be applied
# Name of function, for error messages
fname <- if(is.name(substitute(fun))) deparse(substitute(fun))
else if(is.character(fun)) fun else "fun"
fname <- sQuote(fname)
# R function to apply
if(is.character(fun))
fun <- get(fun)
if(!is.function(fun))
stop(paste("unrecognised format for function", fname))
fargs <- names(formals(fun))
if(!any(c("r", "...") %in% fargs))
stop(paste(fname, "should have an argument",
sQuote("r"), "or", sQuote("...")))
usecorrection <- ("correction" %in% fargs)
# ---------------------------------------------------------------------
# validate other arguments
if((nrank %% 1) != 0)
stop("nrank must be an integer")
stopifnot(nrank > 0 && nrank < nsim/2)
rgiven <- "r" %in% names(list(...))
if(tran <- !is.null(transform)) {
stopifnot(is.expression(transform))
# 'transform' is an expression
aa <- substitute(substitute(tt, list(.=as.name("funX"))),
list(tt=transform))
# 'aa' is a language expression invoking 'substitute'
bb <- eval(parse(text=deparse(aa)))
# 'bb' is an expression obtained by replacing "." by "funX"
transform.funX <- as.call(bb)
transform.funX[[1]] <- as.name("eval.fv")
# 'transform.funX' is an unevaluated call to eval.fv
aa <- substitute(substitute(tt, list(.=as.name("funXsim"))),
list(tt=transform))
bb <- eval(parse(text=deparse(aa)))
transform.funXsim <- as.call(bb)
transform.funXsim[[1]] <- as.name("eval.fv")
}
if(!is.null(ginterval))
stopifnot(is.numeric(ginterval) && length(ginterval) == 2)
# -------------------------------------------------------------------
# Determine clipping window
if(clipdata) {
# Generate one realisation
Xsim <- eval(simexpr, envir=envir)
if(!inherits(Xsim, Xclass))
switch(simtype,
csr=stop(paste("Internal error:", Xobjectname, "not generated")),
rmh=stop(paste("Internal error: rmh did not return an",
Xobjectname)),
kppm=stop(paste("Internal error: simulate.kppm did not return an",
Xobjectname)),
expr=stop(paste("Evaluating the expression", sQuote("simulate"),
"did not yield an", Xobjectname)),
list=stop(paste("Internal error: list entry was not an",
Xobjectname)),
stop(paste("Internal error:", Xobjectname, "not generated"))
)
# Extract window
clipwin <- Xsim$window
if(!is.subset.owin(clipwin, X$window))
warning("Window containing simulated patterns is not a subset of data window")
}
# ---------------------------------------------------------------------
# Evaluate function for data pattern X
Xarg <- if(!clipdata) X else X[clipwin]
funX <- do.call(fun,
resolve.defaults(list(Xarg),
list(...),
if(usecorrection) list(correction="best") else NULL))
if(!inherits(funX, "fv"))
stop(paste("The function", fname,
"must return an object of class", sQuote("fv")))
argname <- fvnames(funX, ".x")
valname <- fvnames(funX, ".y")
if(tran) {
# extract only the recommended value
if(csr)
funX <- funX[, c(argname, valname, "theo")]
else
funX <- funX[, c(argname, valname)]
# apply the transformation to it
funX <- eval(transform.funX)
}
rvals <- funX[[argname]]
fX <- funX[[valname]]
# default domain over which to maximise
alim <- attr(funX, "alim")
if(global && is.null(ginterval))
ginterval <- if(rgiven) range(rvals) else alim
# capture main decision parameters
EnvelopeInfo <- list(call=cl,
Yname=Yname,
valname=valname,
csr=csr,
simtype=simtype,
nrank=nrank,
nsim=nsim,
global=global,
dual=dual,
nsim2=nsim2,
VARIANCE=VARIANCE,
nSD=nSD)
######### SIMULATE #######################
if(verbose) {
action <- if(simtype == "list") "Extracting" else "Generating"
descrip <- switch(simtype,
csr = "simulations of CSR",
rmh = "simulated realisations of fitted Gibbs model",
kppm = "simulated realisations of fitted cluster model",
expr = "simulations by evaluating expression",
list = "point patterns from list",
"simulated patterns")
cat(paste(action, Nsim, descrip, "...\n"))
}
# determine whether simulated point patterns should be saved
catchpatterns <- savepatterns && simtype != "list"
Caughtpatterns <- list()
# allocate space for computed function values
nrvals <- length(rvals)
simvals <- matrix(, nrow=nrvals, ncol=Nsim)
# start simulation loop
for(i in 1:Nsim) {
Xsim <- eval(simexpr, envir=envir)
if(!inherits(Xsim, Xclass))
switch(simtype,
csr=stop(paste("Internal error:", Xobjectname, "not generated")),
rmh=stop(paste("Internal error: rmh did not return an",
Xobjectname)),
kppm=stop(paste("Internal error: simulate.kppm did not return an",
Xobjectname)),
expr=stop(paste("Evaluating the expression", sQuote("simulate"),
"did not yield an", Xobjectname)),
list=stop(paste("Internal error: list entry was not an",
Xobjectname)),
stop(paste("Internal error:", Xobjectname, "not generated"))
)
if(catchpatterns)
Caughtpatterns[[i]] <- Xsim
# apply function
funXsim <- do.call(fun,
resolve.defaults(list(Xsim),
if(rgiven) NULL else list(r=rvals),
list(...),
if(usecorrection) list(correction="best") else NULL))
# sanity checks
if(!inherits(funXsim, "fv"))
stop(paste("When applied to a simulated pattern, the function",
fname, "did not return an object of class",
sQuote("fv")))
argname.sim <- fvnames(funXsim, ".x")
valname.sim <- fvnames(funXsim, ".y")
if(argname.sim != argname)
stop(paste("The objects returned by", fname,
"when applied to a simulated pattern",
"and to the data pattern",
"are incompatible. They have different argument names",
sQuote(argname.sim), "and", sQuote(argname),
"respectively"))
if(valname.sim != valname)
stop(paste("When", fname, "is applied to a simulated pattern",
"it provides an estimate named", sQuote(valname.sim),
"whereas the estimate for the data pattern is named",
sQuote(valname),
". Try using the argument", sQuote("correction"),
"to make them compatible"))
if(tran) {
# extract only the recommended value
if(csr)
funXsim <- funXsim[, c(argname, valname, "theo")]
else
funXsim <- funXsim[, c(argname, valname)]
# apply the transformation to it
funXsim <- eval(transform.funXsim)
}
# extract the values for simulation i
simvals.i <- funXsim[[valname]]
if(length(simvals.i) != nrvals)
stop("Vectors of function values have incompatible lengths")
simvals[ , i] <- funXsim[[valname]]
if(verbose)
progressreport(i, Nsim)
}
## end simulation loop
if(verbose) {
cat("\nDone.\n")
flush.console()
}
# ...........................................................
# save functions and/or patterns if so commanded
if(savefuns) {
alldata <- cbind(rvals, simvals)
simnames <- paste("sim", 1:nsim, sep="")
colnames(alldata) <- c("r", simnames)
alldata <- as.data.frame(alldata)
SimFuns <- fv(alldata,
argu="r",
ylab=attr(funX, "ylab"),
valu="sim1",
fmla= deparse(. ~ r),
alim=attr(funX, "alim"),
labl=names(alldata),
desc=c("distance argument r",
paste("Simulation ", 1:nsim, sep="")))
fvnames(SimFuns, ".") <- simnames
}
if(savepatterns)
SimPats <- if(simtype == "list") SimDataList else Caughtpatterns
######### COMPUTE ENVELOPES #######################
orderstat <- function(x, n) { sort(x)[n] }
if(VARIANCE) {
# pointwise mean, variance etc
simvals[is.infinite(simvals)] <- NA
Ef <- apply(simvals, 1, mean, na.rm=TRUE)
varf <- apply(simvals, 1, var, na.rm=TRUE)
# derived quantities
sd <- sqrt(varf)
stdres <- (fX-Ef)/sd
stdres[!is.finite(stdres)] <- NA
# critical limits
lo <- Ef - nSD * sd
hi <- Ef + nSD * sd
lo.name <- paste("lower", nSD, "sigma critical limit for %s")
hi.name <- paste("upper", nSD, "sigma critical limit for %s")
# put together
if(csr) {
results <- data.frame(r=rvals,
obs=fX,
theo=funX[["theo"]],
lo=lo,
hi=hi)
morestuff <- data.frame(mmean=Ef,
var=varf,
res=fX-Ef,
stdres=stdres)
mslabl <- c("mean(r)", "var(r)", "res(r)", "stdres(r)")
msdesc <- c("sample mean of %s from simulations",
"sample variance of %s from simulations",
"raw residual", "standardised residual")
} else {
results <- data.frame(r=rvals,
obs=fX,
mmean=Ef,
lo=lo,
hi=hi)
morestuff <- data.frame(var=varf,
res=fX-Ef,
stdres=stdres)
mslabl <- c("var(r)", "res(r)", "stdres(r)")
msdesc <- c("sample variance of %s from simulations",
"raw residual", "standardised residual")
}
} else if(!global) {
# POINTWISE ENVELOPES
simvals[is.infinite(simvals)] <- NA
lo <- apply(simvals, 1, orderstat, n=nrank)
hi <- apply(simvals, 1, orderstat, n=nsim-nrank+1)
lo.name <- paste("lower pointwise envelope of %s from simulations")
hi.name <- paste("upper pointwise envelope of %s from simulations")
#
if(csr) {
results <- data.frame(r=rvals,
obs=fX,
theo=funX[["theo"]],
lo=lo,
hi=hi)
} else {
m <- apply(simvals, 1, mean, na.rm=TRUE)
results <- data.frame(r=rvals,
obs=fX,
mmean=m,
lo=lo,
hi=hi)
}
} else {
# SIMULTANEOUS ENVELOPES
domain <- (rvals >= ginterval[1]) & (rvals <= ginterval[2])
funX <- funX[domain, ]
simvals <- simvals[domain, ]
simvals[is.infinite(simvals)] <- NA
if(csr)
theory <- funX[["theo"]]
else {
# use the first 'nsim' results to estimate the mean under H0
theory <- m <- apply(simvals[, 1:nsim], 1, mean, na.rm=TRUE)
# then discard them and use the remaining 'nsim2' simulations
# to construct the envelopes.
simvals <- simvals[, nsim + 1:nsim2]
}
# compute max absolute deviations
deviations <- sweep(simvals, 1, theory)
suprema <- apply(abs(deviations), 2, max, na.rm=TRUE)
# ranked deviations
dmax <- orderstat(suprema, n=nsim-nrank+1)
# simultaneous bands
lo <- theory - dmax
hi <- theory + dmax
lo.name <- "lower critical boundary for %s"
hi.name <- "upper critical boundary for %s"
if(csr)
results <- data.frame(r=rvals[domain],
obs=fX[domain],
theo=theory,
lo=lo,
hi=hi)
else
results <- data.frame(r=rvals[domain],
obs=fX[domain],
mmean=m,
lo=lo,
hi=hi)
}
############ WRAP UP #########################
result <- fv(results,
argu="r",
ylab=attr(funX, "ylab"),
valu="obs",
fmla= deparse(. ~ r),
alim=attr(funX, "alim"),
labl=c("r", "obs(r)", if(csr) "theo(r)" else "mean(r)",
"lo(r)", "hi(r)"),
desc=c("distance argument r",
"observed value of %s for data pattern",
if(csr) "theoretical value of %s for CSR"
else "sample mean of %s from simulations",
lo.name, hi.name),
fname=attr(funX, "fname"),
yexp =attr(funX, "yexp"))
# columns to be plotted by default
dotty <- c("obs", if(csr) "theo" else "mmean", "hi", "lo")
if(VARIANCE) {
# add more stuff
result <- bind.fv(result, morestuff, mslabl, msdesc)
if(csr) dotty <- c(dotty, "mmean")
}
fvnames(result, ".") <- dotty
unitname(result) <- unitname(funX)
class(result) <- c("envelope", class(result))
# tack on envelope information
attr(result, "einfo") <- EnvelopeInfo
# tack on functions and/or patterns if so commanded
if(savefuns)
attr(result, "simfuns") <- SimFuns
if(savepatterns)
attr(result, "simpatterns") <- SimPats
return(result)
}
print.envelope <- function(x, ...) {
e <- attr(x, "einfo")
g <- e$global
csr <- e$csr
simtype <- e$simtype
nr <- e$nrank
nsim <- e$nsim
V <- e$VARIANCE
fname <- deparse(attr(x, "ylab"))
type <- if(V) paste("Pointwise", e$nSD, "sigma") else
if(g) "Simultaneous" else "Pointwise"
cat(paste(type, "critical envelopes for", fname, "\n"))
if(!is.null(valname <- e$valname))
cat(paste("Edge correction:",
dQuote(valname), "\n"))
# determine *actual* type of simulation
descrip <-
if(csr) "simulations of CSR"
else if(!is.null(simtype)) {
switch(simtype,
csr="simulations of CSR",
rmh="simulations of fitted Gibbs model",
kppm="simulations of fitted cluster model",
expr="evaluations of user-supplied expression",
list="point pattern datasets in user-supplied list")
} else "simulations of fitted model"
#
cat(paste("Obtained from", nsim, descrip, "\n"))
#
if(!is.null(e$dual) && e$dual)
cat(paste("Theoretical (i.e. null) mean value of", fname,
"estimated from a separate set of",
e$nsim2, "simulations\n"))
if(!is.null(attr(x, "simfuns")))
cat("(All simulated function values are stored)\n")
if(!is.null(attr(x, "simpatterns")))
cat("(All simulated point patterns are stored)\n")
alpha <- if(g) { nr/(nsim+1) } else { 2 * nr/(nsim+1) }
if(!V) {
# significance interpretation!
cat(paste("Significance level of",
if(!g) "pointwise",
"Monte Carlo test:",
paste(if(g) nr else 2 * nr,
"/", nsim+1, sep=""),
"=", alpha, "\n"))
}
cat(paste("Data:", e$Yname, "\n"))
NextMethod("print")
}
summary.envelope <- function(object, ...) {
e <- attr(object, "einfo")
g <- e$global
nr <- e$nrank
nsim <- e$nsim
csr <- e$csr
simtype <- e$simtype
fname <- deparse(attr(object, "ylab"))
V <- e$VARIANCE
type <- if(V) paste("Pointwise", e$nSD, "sigma") else
if(g) "Simultaneous" else "Pointwise"
cat(paste(type, "critical envelopes for", fname, "\n"))
# determine *actual* type of simulation
descrip <-
if(csr) "simulations of CSR"
else if(!is.null(simtype)) {
switch(simtype,
csr="simulations of CSR",
rmh="simulations of fitted Gibbs model",
kppm="simulations of fitted cluster model",
expr="evaluations of user-supplied expression",
list="point pattern datasets in user-supplied list",
"simulated point patterns")
} else "simulations of fitted model"
#
cat(paste("Obtained from", nsim, descrip, "\n"))
#
if(!is.null(attr(object, "simfuns")))
cat(paste("(All", nsim, "simulated function values",
"are stored in attr(,", dQuote("simfuns"), ") )\n"))
if(!is.null(attr(object, "simpatterns")))
cat(paste("(All", nsim, "simulated point patterns",
"are stored in attr(,", dQuote("simpatterns"), ") )\n"))
#
if(V) {
# nSD envelopes
cat(paste("Envelopes computed as sample mean plus/minus",
e$nSD, "sample standard deviations\n"))
} else {
# critical envelopes
lo.ord <- if(nr == 1) "minimum" else paste(ordinal(nr), "smallest")
hi.ord <- if(nr == 1) "maximum" else paste(ordinal(nsim-nr+1), "largest")
if(g)
cat(paste("Envelopes computed as",
if(csr) "theoretical curve" else "mean of simulations",
"plus/minus", hi.ord,
"simulated value of maximum absolute deviation\n"))
else {
cat(paste("Upper envelope: pointwise", hi.ord,"of simulated curves\n"))
cat(paste("Lower envelope: pointwise", lo.ord,"of simulated curves\n"))
}
alpha <- if(g) { nr/(nsim+1) } else { 2 * nr/(nsim+1) }
cat(paste("Significance level of Monte Carlo test:",
paste(if(g) nr else 2 * nr,
"/", nsim+1, sep=""),
"=", alpha, "\n"))
}
cat(paste("Data:", e$Yname, "\n"))
return(invisible(NULL))
}
