Nothing
R/envelope.R
In spatstat.explore: Exploratory Data Analysis for the 'spatstat' Family
Defines functions
as.data.frame.envelope
resolveEinfo
envelope.envelope
envelope.matrix
summary.envelope
print.envelope
plot.envelope
envelopeEngine
envelope.ppp
simulrecipe
envelope
Documented in
as.data.frame.envelope
envelope
envelopeEngine
envelope.envelope
envelope.matrix
envelope.ppp
plot.envelope
print.envelope
resolveEinfo
simulrecipe
summary.envelope
#
# envelope.R
#
# computes simulation envelopes
#
# $Revision: 2.125 $ $Date: 2023/08/15 08:07:52 $
#
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 or Poisson model
# "kppm": simulated realisation of fitted cluster model
# "slrm": simulated realisation of spatial logistic regression
# "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
#
# pois = TRUE if model is known to be Poisson
#
# constraints = additional information about simulation (e.g. 'with fixed n')
#
# ...................................................................
simulrecipe <- function(type, expr, envir, csr, pois=csr, constraints="") {
if(csr && !pois) warning("Internal error: csr=TRUE but pois=FALSE")
out <- list(type=type,
expr=expr,
envir=envir,
csr=csr,
pois=pois,
constraints=constraints)
class(out) <- "simulrecipe"
out
}
envelope.ppp <-
function(Y, fun=Kest, nsim=99, nrank=1, ...,
funargs=list(), funYargs=funargs,
simulate=NULL, fix.n=FALSE, fix.marks=FALSE,
verbose=TRUE, clipdata=TRUE,
transform=NULL, global=FALSE, ginterval=NULL, use.theory=NULL,
alternative=c("two.sided", "less", "greater"),
scale=NULL, clamp=FALSE,
savefuns=FALSE, savepatterns=FALSE, nsim2=nsim,
VARIANCE=FALSE, nSD=2,
Yname=NULL,
maxnerr=nsim, rejectNA=FALSE, silent=FALSE,
do.pwrong=FALSE,
envir.simul=NULL) {
cl <- short.deparse(sys.call())
if(is.null(Yname)) Yname <- short.deparse(substitute(Y))
if(is.null(fun)) fun <- Kest
envir.user <- if(!is.null(envir.simul)) envir.simul else parent.frame()
envir.here <- sys.frame(sys.nframe())
ismarked <- is.marked(Y)
ismulti <- is.multitype(Y)
fix.marks <- fix.marks && ismarked
if(!is.null(simulate)) {
# ...................................................
# Simulations are determined by 'simulate' argument
if(fix.n || fix.marks)
warning("fix.n and fix.marks were ignored, because 'simulate' was given")
# Processing is deferred to envelopeEngine
simrecipe <- simulate
# Data pattern is argument Y
X <- Y
} else if(!fix.n && !fix.marks) {
# ...................................................
# 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
nY <- npoints(Y)
Ywin <- Y$window
Ymarx <- marks(Y)
# expression that will be evaluated
simexpr <- if(is.null(Ymarx)) {
# unmarked point pattern
expression(rpoispp(Yintens, win=Ywin))
} else if(is.null(dim(Ymarx))) {
# single column of marks
expression({
A <- rpoispp(Yintens, win=Ywin);
j <- sample(nY, npoints(A), replace=TRUE);
A %mark% Ymarx[j]
})
} else {
# multiple columns of marks
expression({
A <- rpoispp(Yintens, win=Ywin);
j <- sample(nY, npoints(A), replace=TRUE);
A %mark% Ymarx[j, , drop=FALSE]
})
}
dont.complain.about(Yintens, Ywin)
# evaluate in THIS environment
simrecipe <- simulrecipe(type = "csr",
expr = simexpr,
envir = envir.here,
csr = TRUE,
pois = TRUE)
} else if(fix.marks) {
# ...................................................
# Data pattern is argument Y
X <- Y
# Realisations of binomial process
# with fixed number of points and fixed marks
# will be generated by runifpoint
nY <- npoints(Y)
Ywin <- as.owin(Y)
Ymarx <- marks(Y)
# expression that will be evaluated
simexpr <- expression(runifpoint(nY, Ywin) %mark% Ymarx)
# suppress warnings from code checkers
dont.complain.about(nY, Ywin, Ymarx)
# simulation constraints (explanatory string)
constraints <- if(ismulti) "with fixed number of points of each type" else
"with fixed number of points and fixed marks"
# evaluate in THIS environment
simrecipe <- simulrecipe(type = "csr",
expr = simexpr,
envir = envir.here,
csr = TRUE,
pois = TRUE,
constraints = constraints)
} else {
# ...................................................
# Data pattern is argument Y
X <- Y
# Realisations of binomial process
# will be generated by runifpoint
nY <- npoints(Y)
Ywin <- as.owin(Y)
Ymarx <- marks(Y)
# expression that will be evaluated
simexpr <- if(is.null(Ymarx)) {
## unmarked
expression(runifpoint(nY, Ywin))
} else if(is.null(dim(Ymarx))) {
## single column of marks
expression({
A <- runifpoint(nY, Ywin);
j <- sample(nY, npoints(A), replace=TRUE);
A %mark% Ymarx[j]
})
} else {
## multiple columns of marks
expression({
A <- runifpoint(nY, Ywin);
j <- sample(nY, npoints(A), replace=TRUE);
A %mark% Ymarx[j, ,drop=FALSE]
})
}
dont.complain.about(nY, Ywin)
# evaluate in THIS environment
simrecipe <- simulrecipe(type = "csr",
expr = simexpr,
envir = envir.here,
csr = TRUE,
pois = TRUE,
constraints = "with fixed number of points")
}
envelopeEngine(X=X, fun=fun, simul=simrecipe,
nsim=nsim, nrank=nrank, ...,
funargs=funargs, funYargs=funYargs,
verbose=verbose, clipdata=clipdata,
transform=transform,
global=global, ginterval=ginterval, use.theory=use.theory,
alternative=alternative, scale=scale, clamp=clamp,
savefuns=savefuns, savepatterns=savepatterns, nsim2=nsim2,
VARIANCE=VARIANCE, nSD=nSD,
Yname=Yname,
maxnerr=maxnerr, rejectNA=rejectNA, silent=silent,
cl=cl, envir.user=envir.user, do.pwrong=do.pwrong)
}
## Code for envelope.ppm, envelope.kppm, envelope.slrm
## is moved to spatstat.model
## .................................................................
## 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, ..., funargs=list(), funYargs=funargs,
verbose=TRUE, clipdata=TRUE,
transform=NULL, global=FALSE, ginterval=NULL, use.theory=NULL,
alternative=c("two.sided", "less", "greater"),
scale=NULL, clamp=FALSE,
savefuns=FALSE, savepatterns=FALSE,
saveresultof=NULL,
weights=NULL,
nsim2=nsim,
VARIANCE=FALSE, nSD=2,
Yname=NULL,
maxnerr=nsim,
rejectNA=FALSE,
silent=FALSE,
maxerr.action=c("fatal", "warn", "null"),
internal=NULL, cl=NULL,
envir.user=envir.user,
expected.arg="r",
do.pwrong=FALSE,
foreignclass=NULL,
collectrubbish=FALSE)
{
#
envir.here <- sys.frame(sys.nframe())
alternative <- match.arg(alternative)
maxerr.action <- match.arg(maxerr.action)
foreignclass <- as.character(foreignclass)
if(length(foreignclass) != 0 && clipdata) {
warning(paste("Ignoring clipdata=TRUE:",
"I don't know how to clip objects of class",
sQuote(paste(foreignclass, collapse=","))))
clipdata <- FALSE
}
# ----------------------------------------------------------
# Determine Simulation
# ----------------------------------------------------------
check.1.integer(nsim)
stopifnot(nsim > 1)
# 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
if(inherits(X, foreignclass)) foreignclass else
stop("Unrecognised class of point pattern")
Xobjectname <- paste("point pattern of class", sQuote(Xclass))
# Option to use weighted average
if(use.weights <- !is.null(weights)) {
# weight can be either a numeric vector or a function
if(is.numeric(weights)) {
compute.weights <- FALSE
weightfun <- NULL
} else if(is.function(weights)) {
compute.weights <- TRUE
weightfun <- weights
weights <- NULL
} else stop("weights should be either a function or a numeric vector")
} else compute.weights <- FALSE
# Undocumented option to generate patterns only.
patterns.only <- identical(internal$eject, "patterns")
# Undocumented option to evaluate 'something' for each simulation
if(savevalues <- !is.null(saveresultof)) {
stopifnot(is.function(saveresultof))
SavedValues <- list()
## might be a function of the pattern only, or both pattern and summary fun
result.depends.both <- (length(formals(saveresultof)) >= 2)
}
# 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
pois <- simul$pois
constraints <- simul$constraints
} 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, or a function to be applied to the data
csr <- FALSE
# override
if(!is.null(icsr <- internal$csr)) csr <- icsr
pois <- csr
constraints <- ""
# 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.function(simulate)) {
## User-supplied function 'simulate' will be repeatedly evaluated on X
simtype <- "func"
simexpr <- expression(simulate(X))
envir <- envir.here
} else if(is.list(simulate) &&
all(sapply(simulate, inherits, what=Xclass))) {
#' The user-supplied list of point patterns will be used
simtype <- "list"
SimDataList <- simulate
#' expression that will be evaluated
simexpr <- expression(SimDataList[[i+nerr]])
dont.complain.about(SimDataList)
envir <- envir.here
#' ensure that `i' is defined
i <- 1L
nerr <- 0L
maxnerr <- min(length(SimDataList)-nsim, maxnerr)
#' any messages?
if(!is.null(mess <- attr(simulate, "internal"))) {
# determine whether these point patterns are realisations of CSR
csr <- isTRUE(mess$csr)
}
} else if(is.list(simulate) &&
all(sapply(simulate, is.list)) &&
all(lengths(simulate) == 1) &&
all(sapply((elements <- lapply(simulate, "[[", i=1)), inherits, what=Xclass))) {
#' malformed argument: list(list(ppp), list(ppp), ....)
SimDataList <- elements
simtype <- "list"
#' expression that will be evaluated
simexpr <- expression(SimDataList[[i+nerr]])
dont.complain.about(SimDataList)
envir <- envir.here
#' ensure that `i' is defined
i <- 1L
nerr <- 0L
maxnerr <- min(length(SimDataList)-nsim, maxnerr)
#' any messages?
if(!is.null(mess <- attr(simulate, "internal"))) {
# determine whether these point patterns are realisations of CSR
csr <- isTRUE(mess$csr)
}
} else stop(paste(sQuote("simulate"),
"should be an expression,",
"or a list of point patterns of the same kind as X"))
}
# -------------------------------------------------------------------
# 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)),
slrm=stop(paste("Internal error: simulate.slrm did not return an",
Xobjectname)),
expr=stop(paste("Evaluating the expression", sQuote("simulate"),
"did not yield an", Xobjectname)),
func=stop(paste("Evaluating the function", 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")
}
# ------------------------------------------------------------------
# Summary function to be applied
# ------------------------------------------------------------------
if(is.null(fun))
stop("Internal error: fun is NULL")
# Name of function, for error messages
fname <- if(is.name(substitute(fun))) short.deparse(substitute(fun)) else
if(is.character(fun)) fun else "fun"
fname <- sQuote(fname)
# R function to apply
if(is.character(fun)) {
gotfun <- try(get(fun, mode="function"))
if(inherits(gotfun, "try-error"))
stop(paste("Could not find a function named", sQuote(fun)))
fun <- gotfun
} else if(!is.function(fun))
stop(paste("unrecognised format for function", fname))
fargs <- names(formals(fun))
if(!any(c(expected.arg, "...") %in% fargs))
stop(paste(fname, "should have",
ngettext(length(expected.arg), "an argument", "arguments"),
"named", commasep(sQuote(expected.arg)),
"or a", sQuote("..."), "argument"))
usecorrection <- any(c("correction", "...") %in% fargs)
# ---------------------------------------------------------------------
# validate other arguments
if((nrank %% 1) != 0)
stop("nrank must be an integer")
if((nsim %% 1) != 0)
stop("nsim must be an integer")
stopifnot(nrank > 0 && nrank < nsim/2)
arg.given <- any(expected.arg %in% names(list(...)))
if(tran <- !is.null(transform)) {
stopifnot(is.expression(transform))
# prepare expressions to be evaluated each time
transform.funX <- inject.expr("with(funX,.)", transform)
transform.funXsim <- inject.expr("with(funXsim,.)", transform)
# .... old code using 'eval.fv' ......
# transform.funX <- dotexpr.to.call(transform, "funX", "eval.fv")
# transform.funXsim <- dotexpr.to.call(transform, "funXsim", "eval.fv")
# 'transform.funX' and 'transform.funXsim' are unevaluated calls to eval.fv
}
if(!is.null(ginterval))
stopifnot(is.numeric(ginterval) && length(ginterval) == 2)
# ---------------------------------------------------------------------
# Evaluate function for data pattern X
# ------------------------------------------------------------------
Xarg <- if(!clipdata) X else X[clipwin]
corrx <- if(usecorrection) list(correction="best") else NULL
dont.complain.about(Xarg)
funX <- do.call(fun,
resolve.defaults(list(quote(Xarg)),
list(...),
funYargs,
corrx))
if(!inherits(funX, "fv"))
stop(paste("The function", fname,
"must return an object of class", sQuote("fv")))
## catch 'conservation' parameters
conserveargs <- attr(funX, "conserve")
if(!is.null(conserveargs) && !any(c("conserve", "...") %in% fargs))
stop(paste("In this usage, the function", fname,
"should have an argument named 'conserve' or '...'"))
## warn about 'dangerous' arguments
if(!is.null(dang <- attr(funX, "dangerous")) &&
any(uhoh <- dang %in% names(list(...)))) {
nuh <- sum(uhoh)
warning(paste("Envelope may be invalid;",
ngettext(nuh, "argument", "arguments"),
commasep(sQuote(dang[uhoh])),
ngettext(nuh, "appears", "appear"),
"to have been fixed."),
call.=FALSE)
}
argname <- fvnames(funX, ".x")
valname <- fvnames(funX, ".y")
has.theo <- "theo" %in% fvnames(funX, "*")
csr.theo <- csr && has.theo
use.theory <- if(is.null(use.theory)) csr.theo else (use.theory && has.theo)
if(tran) {
# extract only the recommended value
if(use.theory)
funX <- funX[, c(argname, valname, "theo")]
else
funX <- funX[, c(argname, valname)]
# apply the transformation to it
funX <- eval(transform.funX)
}
argvals <- funX[[argname]]
# fX <- funX[[valname]]
arg.desc <- attr(funX, "desc")[match(argname, colnames(funX))]
# default domain over which to maximise
alim <- attr(funX, "alim")
if(global && is.null(ginterval))
ginterval <- if(arg.given || is.null(alim)) range(argvals) else alim
#--------------------------------------------------------------------
# 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 && !use.theory && !VARIANCE)
if(dual) {
check.1.integer(nsim2)
stopifnot(nsim2 >= 1)
}
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(patterns.only) {
# 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 = paste("simulated realisations of fitted",
if(pois) "Poisson" else "Gibbs",
"model"),
kppm = "simulated realisations of fitted cluster model",
slrm = "simulated realisations of spatial logistic regression model",
expr = "simulations by evaluating expression",
func = "simulations by evaluating function",
list = "point patterns from list",
"simulated realisations")
if(!is.null(constraints) && nzchar(constraints))
descrip <- paste(descrip, constraints)
explan <- if(dual) paren(paste(nsim2, "to estimate the mean and",
nsim, "to calculate envelopes")) else ""
splat(action, Nsim, descrip, explan, "...")
}
XsimList <- list()
# start simulation loop
sstate <- list()
for(i in 1:Nsim) {
if(verbose) sstate <- progressreport(i, Nsim, state=sstate)
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))
},
slrm={
stop(paste("Internal error: simulate.slrm did not return an",
Xobjectname))
},
expr={
stop(paste("Evaluating the expression", sQuote("simulate"),
"did not yield an", Xobjectname))
},
func={
stop(paste("Evaluating the function", 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)
}
# capture main decision parameters
envelopeInfo <- list(call=cl,
Yname=Yname,
valname=valname,
csr=csr,
csr.theo=csr.theo,
use.theory=use.theory,
pois=pois,
simtype=simtype,
constraints=constraints,
nrank=nrank,
nsim=nsim,
Nsim=Nsim,
global=global,
ginterval=ginterval,
dual=dual,
nsim2=nsim2,
VARIANCE=VARIANCE,
nSD=nSD,
alternative=alternative,
scale=scale,
clamp=clamp,
use.weights=use.weights,
do.pwrong=do.pwrong)
# ----------------------------------------
######### SIMULATE #######################
# ----------------------------------------
if(verbose) {
action <- if(simtype == "list") "Extracting" else "Generating"
descrip <- switch(simtype,
csr = "simulations of CSR",
rmh = paste("simulated realisations of fitted",
if(pois) "Poisson" else "Gibbs",
"model"),
kppm = "simulated realisations of fitted cluster model",
slrm = "simulated realisations of fitted spatial logistic regression model",
expr = "simulations by evaluating expression",
func = "simulations by evaluating function",
list = "point patterns from list",
"simulated patterns")
if(!is.null(constraints) && nzchar(constraints))
descrip <- paste(descrip, constraints)
explan <- if(dual) paren(paste(nsim2, "to estimate the mean and",
nsim, "to calculate envelopes")) else ""
splat(action, Nsim, descrip, explan, "...")
}
# determine whether simulated point patterns should be saved
catchpatterns <- savepatterns && simtype != "list"
Caughtpatterns <- list()
# allocate space for computed function values
nargvals <- length(argvals)
simvals <- matrix(, nrow=nargvals, ncol=Nsim)
# allocate space for weights to be computed
if(compute.weights)
weights <- numeric(Nsim)
# inferred values of function argument 'r' or equivalent parameters
if(identical(expected.arg, "r") && identical(argname, "r")) {
# Kest, etc
inferred.r.args <- list(r=argvals)
} else if(identical(expected.arg, c("rmax", "nrval"))) {
# K3est, etc
inferred.r.args <- list(rmax=max(argvals), nrval=length(argvals))
} else if(any(c(argname, "...") %in% fargs)) {
## assume it accepts the vector of argument values
inferred.r.args <- structure(list(argvals), names=argname)
} else {
stop(paste("Don't know how to infer values of",
if(length(expected.arg))
commasep(sQuote(expected.arg)) else "function argument"))
}
# arguments for function when applied to simulated patterns
funargs <-
resolve.defaults(funargs,
inferred.r.args,
list(...),
conserveargs,
if(usecorrection) list(correction="best") else NULL)
# reject simulated pattern if function values are all NA (etc)
rejectNA <- isTRUE(rejectNA)
# start simulation loop
nerr <- 0
gaveup <- FALSE
if(verbose) pstate <- list()
for(i in 1:Nsim) {
## safely generate a random pattern and apply function
success <- FALSE
while(!success && !gaveup) {
Xsim <- eval(simexpr, envir=envir)
## check valid point pattern
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)),
slrm=stop(paste("Internal error:",
"simulate.slrm did not return an",
Xobjectname)),
expr=stop(paste("Evaluating the expression", sQuote("simulate"),
"did not yield an", Xobjectname)),
func=stop(paste("Evaluating the function", 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
if(savevalues && !result.depends.both)
SavedValues[[i]] <- saveresultof(Xsim)
if(compute.weights) {
wti <- weightfun(Xsim)
if(!is.numeric(wti))
stop("weightfun did not return a numeric value")
if(length(wti) != 1L)
stop("weightfun should return a single numeric value")
weights[i] <- wti
}
## apply function safely
funXsim <- try(do.call(fun, c(list(Xsim), funargs)), silent=silent)
success <-
!inherits(funXsim, "try-error") &&
inherits(funXsim, "fv") &&
(!rejectNA || any(is.finite(funXsim[[valname]])))
if(!success) {
#' error in computing summary function
nerr <- nerr + 1L
if(nerr > maxnerr) {
gaveup <- TRUE
errtype <- if(rejectNA) "fatal errors or NA function values"
if(simtype == "list") {
whinge <- paste("Exceeded maximum possible number of errors",
"when evaluating summary function:",
length(SimDataList), "patterns provided,",
nsim, "patterns required,",
nerr, ngettext(nerr, "pattern", "pattern"),
"rejected due to", errtype)
} else {
whinge <- paste("Exceeded maximum permissible number of",
errtype,
paren(paste("maxnerr =", maxnerr)),
"when evaluating summary function",
"for simulated point patterns")
}
switch(maxerr.action,
fatal = stop(whinge, call.=FALSE),
warn = warning(whinge, call.=FALSE),
null = {})
} else if(!silent) cat("[retrying]\n")
}
#' ..... end of while(!success) ................
}
if(gaveup) break; # exit loop now
## sanity checks
if(i == 1L) {
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"))
rfunX <- with(funX, ".x")
rfunXsim <- with(funXsim, ".x")
if(!identical(rfunX, rfunXsim))
stop(paste("When", fname, "is applied to a simulated pattern,",
"the values of the argument", sQuote(argname.sim),
"are different from those used for the data."))
}
if(savevalues && result.depends.both)
SavedValues[[i]] <- saveresultof(Xsim, funXsim)
if(tran) {
# extract only the recommended value
if(use.theory)
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) != nargvals)
stop("Vectors of function values have incompatible lengths")
simvals[ , i] <- funXsim[[valname]]
if(verbose)
pstate <- progressreport(i, Nsim, state=pstate)
if(collectrubbish) {
rm(Xsim)
rm(funXsim)
gc()
}
}
## end simulation loop
if(verbose) {
cat("\nDone.\n")
flush.console()
}
# ...........................................................
# save functions and/or patterns if so commanded
if(!gaveup) {
if(savefuns) {
alldata <- cbind(argvals, simvals)
simnames <- paste("sim", 1:Nsim, sep="")
colnames(alldata) <- c(argname, simnames)
alldata <- as.data.frame(alldata)
SimFuns <- fv(alldata,
argu=argname,
ylab=attr(funX, "ylab"),
valu="sim1",
fmla= paste(". ~", argname),
alim=attr(funX, "alim"),
labl=names(alldata),
desc=c(arg.desc,
paste("Simulation ", 1:Nsim, sep="")),
fname=attr(funX, "fname"),
yexp=attr(funX, "yexp"),
unitname=unitname(funX))
fvnames(SimFuns, ".") <- simnames
}
if(savepatterns)
SimPats <- if(simtype == "list") SimDataList else Caughtpatterns
}
######### COMPUTE ENVELOPES #######################
etype <- if(global) "global" else if(VARIANCE) "variance" else "pointwise"
if(dual) {
jsim <- 1:nsim
jsim.mean <- nsim + 1:nsim2
} else {
jsim <- jsim.mean <- NULL
}
result <- envelope.matrix(simvals, funX=funX,
jsim=jsim, jsim.mean=jsim.mean,
type=etype, alternative=alternative,
scale=scale, clamp=clamp,
csr=csr, use.theory=use.theory,
nrank=nrank, ginterval=ginterval, nSD=nSD,
Yname=Yname, do.pwrong=do.pwrong,
weights=weights, gaveup=gaveup)
## tack on envelope information
attr(result, "einfo") <- resolve.defaults(envelopeInfo,
attr(result, "einfo"))
if(!gaveup) {
## tack on functions and/or patterns if so commanded
if(savefuns)
attr(result, "simfuns") <- SimFuns
if(savepatterns) {
attr(result, "simpatterns") <- SimPats
attr(result, "datapattern") <- X
}
## undocumented - tack on values of some other quantity
if(savevalues) {
attr(result, "simvalues") <- SavedValues
attr(result, "datavalue") <-
if(result.depends.both) saveresultof(X, funX) else saveresultof(X)
}
}
## save function weights
if(use.weights)
attr(result, "weights") <- weights
return(result)
}
plot.envelope <- function(x, ..., main) {
if(missing(main)) main <- short.deparse(substitute(x))
shade.given <- ("shade" %in% names(list(...)))
shade.implied <- !is.null(fvnames(x, ".s"))
if(!(shade.given || shade.implied)) {
# ensure x has default 'shade' attribute
# (in case x was produced by an older version of spatstat)
if(all(c("lo", "hi") %in% colnames(x)))
fvnames(x, ".s") <- c("lo", "hi")
else warning("Unable to determine shading for envelope")
}
NextMethod("plot", main=main)
}
print.envelope <- function(x, ...) {
e <- attr(x, "einfo")
g <- e$global
csr <- e$csr
pois <- e$pois
if(is.null(pois)) pois <- csr
simtype <- e$simtype
constraints <- e$constraints
nr <- e$nrank
nsim <- e$nsim
V <- e$VARIANCE
fname <- flat.deparse(attr(x, "ylab"))
type <- if(V) paste("Pointwise", e$nSD, "sigma") else
if(g) "Simultaneous" else "Pointwise"
splat(type, "critical envelopes for", fname,
"\nand observed value for", sQuote(e$Yname))
if(!is.null(valname <- e$valname) && waxlyrical('extras'))
splat("Edge correction:", dQuote(valname))
## determine *actual* type of simulation
descrip <-
if(csr) "simulations of CSR"
else if(!is.null(simtype)) {
switch(simtype,
csr="simulations of CSR",
rmh=paste("simulations of fitted",
if(pois) "Poisson" else "Gibbs",
"model"),
kppm="simulations of fitted cluster model",
slrm="simulations of fitted spatial logistic regression model",
expr="evaluations of user-supplied expression",
func="evaluations of user-supplied function",
list="point pattern datasets in user-supplied list",
funs="columns of user-supplied data")
} else "simulations of fitted model"
if(!is.null(constraints) && nzchar(constraints))
descrip <- paste(descrip, constraints)
#
splat("Obtained from", nsim, descrip)
#
if(waxlyrical('extras')) {
dual <- isTRUE(e$dual)
usetheory <- isTRUE(e$use.theory)
hownull <- if(usetheory) {
"(known exactly)"
} else if(dual) {
paste("(estimated from a separate set of",
e$nsim2, "simulations)")
} else NULL
formodel <- if(csr) "for CSR" else NULL
if(g) {
splat("Envelope based on maximum deviation of", fname,
"from null value", formodel, hownull)
} else if(dual) {
splat("Null value of", fname, formodel, hownull)
}
if(!is.null(attr(x, "simfuns")))
splat("(All simulated function values are stored)")
if(!is.null(attr(x, "simpatterns")))
splat("(All simulated point patterns are stored)")
}
splat("Alternative:", e$alternative)
if(!V && waxlyrical('extras')) {
## significance interpretation!
alpha <- if(g) { nr/(nsim+1) } else { 2 * nr/(nsim+1) }
splat("Significance level of",
if(g) "simultaneous" else "pointwise",
"Monte Carlo test:",
paste0(if(g) nr else 2 * nr,
"/", nsim+1),
"=", signif(alpha, 3))
}
if(waxlyrical('gory') && !is.null(pwrong <- attr(x, "pwrong"))) {
splat("\t[Estimated significance level of pointwise excursions:",
paste0("pwrong=", signif(pwrong, 3), "]"))
}
NextMethod("print")
}
summary.envelope <- function(object, ...) {
e <- attr(object, "einfo")
g <- e$global
V <- e$VARIANCE
nr <- e$nrank
nsim <- e$nsim
csr <- e$csr
pois <- e$pois
if(is.null(pois)) pois <- csr
simtype <- e$simtype
constraints <- e$constraints
alternative <- e$alternative
use.theory <- e$use.theory
has.theo <- "theo" %in% fvnames(object, "*")
csr.theo <- csr && has.theo
use.theory <- if(is.null(use.theory)) csr.theo else (use.theory && has.theo)
fname <- short.deparse(attr(object, "ylab"))
type <- if(V) paste("Pointwise", e$nSD, "sigma") else
if(g) "Simultaneous" else "Pointwise"
splat(type, "critical envelopes for", fname,
"\nand observed value for", sQuote(e$Yname))
# determine *actual* type of simulation
descrip <-
if(csr) "simulations of CSR"
else if(!is.null(simtype)) {
switch(simtype,
csr="simulations of CSR",
rmh=paste("simulations of fitted",
if(pois) "Poisson" else "Gibbs",
"model"),
kppm="simulations of fitted cluster model",
slrm="simulations of fitted spatial logistic regression model",
expr="evaluations of user-supplied expression",
func="evaluations of user-supplied function",
list="point pattern datasets in user-supplied list",
funs="columns of user-supplied data",
"simulated point patterns")
} else "simulations of fitted model"
if(!is.null(constraints) && nzchar(constraints))
descrip <- paste(descrip, constraints)
#
splat("Obtained from", nsim, descrip)
#
if(waxlyrical('extras')) {
if(!is.null(e$dual) && e$dual)
splat("Theoretical (i.e. null) mean value of", fname,
"estimated from a separate set of",
e$nsim2, "simulations")
if(!is.null(attr(object, "simfuns")))
splat("(All", nsim, "simulated function values",
"are stored in attr(,", dQuote("simfuns"), ") )")
if(!is.null(attr(object, "simpatterns")))
splat("(All", nsim, "simulated point patterns",
"are stored in attr(,", dQuote("simpatterns"), ") )")
}
#
splat("Alternative:", alternative)
if(V) {
# nSD envelopes
splat(switch(alternative,
two.sided = "Envelopes",
"Critical boundary"),
"computed as sample mean",
switch(alternative,
two.sided="plus/minus",
less="minus",
greater="plus"),
e$nSD, "sample standard deviations")
} else {
# critical envelopes
lo.ord <- if(nr == 1L) "minimum" else paste(ordinal(nr), "smallest")
hi.ord <- if(nr == 1L) "maximum" else paste(ordinal(nr), "largest")
if(g)
splat(switch(alternative,
two.sided = "Envelopes",
"Critical boundary"),
"computed as",
if(use.theory) "theoretical curve" else "mean of simulations",
switch(alternative,
two.sided="plus/minus",
less="minus",
greater="plus"),
hi.ord,
"simulated value of maximum",
switch(alternative,
two.sided="absolute",
less="negative",
greater="positive"),
"deviation")
else {
if(alternative != "less")
splat("Upper envelope: pointwise", hi.ord, "of simulated curves")
if(alternative != "greater")
splat("Lower envelope: pointwise", lo.ord, "of simulated curves")
}
symmetric <- (alternative == "two.sided") && !g
alpha <- if(!symmetric) { nr/(nsim+1) } else { 2 * nr/(nsim+1) }
splat("Significance level of Monte Carlo test:",
paste0(if(!symmetric) nr else 2 * nr,
"/", nsim+1),
"=", alpha)
}
splat("Data:", e$Yname)
return(invisible(NULL))
}
# envelope.matrix
# core functionality to compute envelope values
# theory = funX[["theo"]]
# observed = fX
envelope.matrix <- function(Y, ...,
argvals=rvals, rvals=NULL, ## rvals is old name
observed=NULL, theory=NULL,
funX=NULL,
nsim=NULL, nsim2=NULL,
jsim=NULL, jsim.mean=NULL,
type=c("pointwise", "global", "variance"),
alternative=c("two.sided", "less", "greater"),
scale = NULL, clamp=FALSE,
csr=FALSE, use.theory = csr,
nrank=1, ginterval=NULL, nSD=2,
savefuns=FALSE,
check=TRUE,
Yname=NULL,
argname=NULL,
arg.desc=NULL,
do.pwrong=FALSE,
weights=NULL,
precomputed=NULL,
gaveup=FALSE) {
if(is.null(Yname))
Yname <- short.deparse(substitute(Y))
type <- match.arg(type)
alternative <- match.arg(alternative)
if(!is.null(funX))
stopifnot(is.fv(funX))
pwrong <- NULL
use.weights <- !is.null(weights)
cheat <- !is.null(precomputed)
if(is.null(argvals) && is.null(observed) && !is.null(funX)) {
## assume funX is summary function for observed data
argvals <- with(funX, .x)
observed <- with(funX, .y)
argname.orig <- fvnames(funX, ".x")
if(is.null(argname)) argname <- argname.orig
if(is.null(arg.desc))
arg.desc <- attr(funX, "desc")[match(argname.orig, colnames(funX))]
theory <- if(use.theory) (theory %orifnull% funX[["theo"]]) else NULL
if(check) stopifnot(nrow(funX) == nrow(Y))
} else {
## construct envelope from raw data
if(is.null(argname)) argname <- "r"
if(is.null(arg.desc)) arg.desc <- "distance argument r"
if(check) {
## validate vectors of data
if(is.null(argvals)) stop("argvals must be supplied")
if(is.null(observed)) stop("observed must be supplied")
stopifnot(length(argvals) == nrow(Y))
stopifnot(length(observed) == length(argvals))
}
}
use.theory <- use.theory && !is.null(theory)
if(use.theory && check) stopifnot(length(theory) == length(argvals))
simvals <- Y
fX <- observed
if(!is.null(funX)) {
atr <- attributes(funX)
fname <- atr$fname
yexp <- atr$yexp
} else {
fname <- "f"
yexp <- substitute(f(r), list(r=as.name(argname)))
atr <- list(alim=range(argvals),
ylab=yexp,
yexp=yexp,
fname="f")
}
NAvector <- rep(NA_real_, length(argvals))
if(!cheat) {
## ................ standard calculation .....................
## validate weights
if(use.weights && !gaveup)
check.nvector(weights, ncol(simvals),
things="simulated functions", naok=TRUE, vname="weights")
## determine numbers of columns used
Ncol <- if(!gaveup) ncol(simvals) else Inf
if(Ncol < 2)
stop("Need at least 2 columns of function values")
## all columns are used unless 'nsim' or 'jsim' given.
if(!(is.null(nsim) && is.null(jsim))) {
if(is.null(jsim)) {
jsim <- 1:nsim
} else if(is.null(nsim)) {
nsim <- length(jsim)
} else stopifnot(length(jsim) == nsim)
if(nsim > Ncol)
stop(paste(nsim, "simulations are not available; only",
Ncol, "columns provided"))
}
## nsim2 or jsim.mean may be given, and imply dual calculation
if(!(is.null(nsim2) && is.null(jsim.mean))) {
if(is.null(jsim.mean)) {
jsim.mean <- setdiff(seq_len(Ncol), jsim)[1:nsim2]
} else if(is.null(nsim2)) {
nsim2 <- length(jsim.mean)
} else stopifnot(length(jsim.mean) == nsim2)
if(nsim + nsim2 > Ncol)
stop(paste(nsim, "+", nsim2, "=", nsim+nsim2,
"simulations are not available; only",
Ncol, "columns provided"))
if(length(intersect(jsim, jsim.mean)))
warning("Internal warning: Indices in jsim and jsim.mean overlap")
}
restrict.columns <- !is.null(jsim)
dual <- !is.null(jsim.mean)
} else {
## ................ precomputed values ..................
## validate weights
if(use.weights)
check.nvector(weights, nsim,
things="simulations", naok=TRUE, vname="weights")
restrict.columns <- FALSE
dual <- FALSE
}
shadenames <- NULL
nsim.mean <- NULL
switch(type,
pointwise = {
## ....... POINTWISE ENVELOPES ...............................
if(gaveup) {
lo <- hi <- NAvector
} else if(cheat) {
stopifnot(checkfields(precomputed, c("lo", "hi")))
lo <- precomputed$lo
hi <- precomputed$hi
} else {
simvals[is.infinite(simvals)] <- NA
if(restrict.columns) {
simvals <- simvals[,jsim]
if(use.weights) weights <- weights[jsim]
}
nsim <- ncol(simvals)
if(nrank == 1L) {
lohi <- apply(simvals, 1L, range)
} else {
lohi <- apply(simvals, 1L,
# function(x, n) { sort(x)[n] },
orderstats,
k=c(nrank, nsim-nrank+1L))
}
lo <- lohi[1L,]
hi <- lohi[2L,]
}
lo.name <- "lower pointwise envelope of %s from simulations"
hi.name <- "upper pointwise envelope of %s from simulations"
##
if(!gaveup)
switch(alternative,
two.sided = { },
less = {
hi <- rep.int(Inf, length(hi))
hi.name <- "infinite upper limit"
},
greater = {
lo <- rep.int(-Inf, length(lo))
lo.name <- "infinite lower limit"
})
if(use.theory) {
results <- data.frame(r=argvals,
obs=fX,
theo=theory,
lo=lo,
hi=hi)
} else {
m <- if(gaveup) NAvector else
if(cheat) precomputed$mmean else
if(!use.weights) apply(simvals, 1L, mean, na.rm=TRUE) else
apply(simvals, 1L, weighted.mean, w=weights, na.rm=TRUE)
results <- data.frame(r=argvals,
obs=fX,
mmean=m,
lo=lo,
hi=hi)
}
colnames(results)[1] <- argname
shadenames <- c("lo", "hi")
if(do.pwrong) {
## estimate the p-value for the 'wrong test'
if(gaveup) {
pwrong <- NA_real_
} else if(cheat) {
pwrong <- precomputed$pwrong
do.pwrong <- !is.null(pwrong) && !badprobability(pwrong, FALSE)
} else {
dataranks <- t(apply(simvals, 1, rank, ties.method="random"))
upper.signif <- (dataranks <= nrank)
lower.signif <- (dataranks >= nsim-nrank+1L)
is.signif <- switch(alternative,
less = lower.signif,
greater = upper.signif,
two.sided = lower.signif | upper.signif)
is.signif.somewhere <- matcolany(is.signif)
pwrong <- sum(is.signif.somewhere)/nsim
}
}
},
global = {
## ..... SIMULTANEOUS ENVELOPES ..........................
if(gaveup) {
lo <- hi <- reference <- NAvector
} else if(cheat) {
## ... use precomputed values ..
stopifnot(checkfields(precomputed, c("lo", "hi")))
lo <- precomputed$lo
hi <- precomputed$hi
if(use.theory) {
reference <- theory
} else {
stopifnot(checkfields(precomputed, "mmean"))
reference <- precomputed$mmean
}
domain <- rep.int(TRUE, length(argvals))
} else {
## ... normal case: compute envelopes from simulations
if(!is.null(ginterval)) {
domain <- (argvals >= ginterval[1L]) & (argvals <= ginterval[2L])
funX <- funX[domain, ]
simvals <- simvals[domain, ]
} else domain <- rep.int(TRUE, length(argvals))
simvals[is.infinite(simvals)] <- NA
if(use.theory) {
reference <- theory[domain]
if(restrict.columns) {
simvals <- simvals[, jsim]
if(use.weights) weights <- weights[jsim]
}
} else if(dual) {
# Estimate the mean from one set of columns
# Form envelopes from another set of columns
simvals.mean <- simvals[, jsim.mean]
# mmean <-
reference <-
if(!use.weights) apply(simvals.mean, 1L, mean, na.rm=TRUE) else
apply(simvals.mean, 1L, weighted.mean, w=weights[jsim.mean],
na.rm=TRUE)
nsim.mean <- ncol(simvals.mean)
# retain only columns used for envelope
simvals <- simvals[, jsim]
} else {
# Compute the mean and envelopes using the same data
if(restrict.columns) {
simvals <- simvals[, jsim]
if(use.weights) weights <- weights[jsim]
}
# mmean <-
reference <-
if(!use.weights) apply(simvals, 1L, mean, na.rm=TRUE) else
apply(simvals, 1L, weighted.mean, w=weights, na.rm=TRUE)
}
nsim <- ncol(simvals)
# compute deviations
deviations <- sweep(simvals, 1L, reference)
deviations <-
switch(alternative,
two.sided = abs(deviations),
greater = if(clamp) pmax(0, deviations) else deviations,
less = if(clamp) pmax(0, -deviations) else (-deviations))
deviations <- matrix(deviations,
nrow=nrow(simvals), ncol=ncol(simvals))
## rescale ?
sc <- 1
if(!is.null(scale)) {
stopifnot(is.function(scale))
sc <- scale(argvals)
sname <- paste0("scale", paren(argname))
tname <- paste("values of", argname)
ans <- check.nvector(sc, length(argvals), things=tname,
fatal=FALSE, vname=sname)
if(!ans)
stop(attr(ans, "whinge"), call.=FALSE)
if(any(bad <- (sc <= 0))) {
## issue a warning unless this only happens at r=0
if(any(bad[argvals > 0]))
warning(paste("Some values of", sname,
"were negative or zero:",
"scale was reset to 1 for these values"),
call.=FALSE)
sc[bad] <- 1
}
deviations <- sweep(deviations, 1L, sc, "/")
}
## compute max (scaled) deviations
suprema <- apply(deviations, 2L, max, na.rm=TRUE)
# ranked deviations
dmax <- sort(suprema)[nsim-nrank+1L]
# simultaneous bands
lo <- reference - sc * dmax
hi <- reference + sc * dmax
}
lo.name <- "lower critical boundary for %s"
hi.name <- "upper critical boundary for %s"
if(!gaveup)
switch(alternative,
two.sided = { },
less = {
hi <- rep.int(Inf, length(hi))
hi.name <- "infinite upper boundary"
},
greater = {
lo <- rep.int(-Inf, length(lo))
lo.name <- "infinite lower boundary"
})
if(use.theory) {
results <- data.frame(r=argvals[domain],
obs=fX[domain],
theo=reference,
lo=lo,
hi=hi)
} else {
results <- data.frame(r=argvals[domain],
obs=fX[domain],
mmean=reference,
lo=lo,
hi=hi)
}
colnames(results)[1] <- argname
shadenames <- c("lo", "hi")
if(do.pwrong)
warning(paste("Argument", sQuote("do.pwrong=TRUE"), "ignored;",
"it is not relevant to global envelopes"))
},
variance={
## ....... POINTWISE MEAN, VARIANCE etc ......................
if(gaveup) {
Ef <- varf <- NAvector
} else if(cheat) {
# .... use precomputed values ....
stopifnot(checkfields(precomputed, c("Ef", "varf")))
Ef <- precomputed$Ef
varf <- precomputed$varf
} else {
## .... normal case: compute from simulations
simvals[is.infinite(simvals)] <- NA
if(restrict.columns) {
simvals <- simvals[, jsim]
if(use.weights) weights <- weights[jsim]
}
nsim <- ncol(simvals)
if(!use.weights) {
Ef <- apply(simvals, 1L, mean, na.rm=TRUE)
varf <- apply(simvals, 1L, var, na.rm=TRUE)
} else {
Ef <- apply(simvals, 1L, weighted.mean, w=weights, na.rm=TRUE)
varf <- apply(simvals, 1L, weighted.var, w=weights, na.rm=TRUE)
}
}
if(gaveup) {
sd <- stdres <- lo <- hi <- loCI <- hiCI <- NAvector
} else {
## derived quantities
sd <- sqrt(varf)
stdres <- (fX-Ef)/sd
stdres[!is.finite(stdres)] <- NA
## critical limits
lo <- Ef - nSD * sd
hi <- Ef + nSD * sd
## confidence interval
loCI <- Ef - nSD * sd/sqrt(nsim)
hiCI <- Ef + nSD * sd/sqrt(nsim)
}
lo.name <- paste("lower", nSD, "sigma critical limit for %s")
hi.name <- paste("upper", nSD, "sigma critical limit for %s")
loCI.name <- paste("lower", nSD, "sigma confidence bound",
"for mean of simulated %s")
hiCI.name <- paste("upper", nSD, "sigma confidence bound",
"for mean of simulated %s")
##
if(!gaveup)
switch(alternative,
two.sided = { },
less = {
hi <- hiCI <- rep.int(Inf, length(hi))
hi.name <- "infinite upper boundary"
hiCI.name <- "infinite upper confidence limit"
},
greater = {
lo <- loCI <- rep.int(-Inf, length(lo))
lo.name <- "infinite lower boundary"
loCI.name <- "infinite lower confidence limit"
})
## put together
if(use.theory) {
results <- data.frame(r=argvals,
obs=fX,
theo=theory,
lo=lo,
hi=hi)
colnames(results)[1] <- argname
shadenames <- c("lo", "hi")
morestuff <- data.frame(mmean=Ef,
var=varf,
res=fX-Ef,
stdres=stdres,
loCI=loCI,
hiCI=hiCI)
loCIlabel <-
if(alternative == "greater" && !gaveup) "-infinity" else
makefvlabel(NULL, NULL, fname, "loCI", argname=argname)
hiCIlabel <-
if(alternative == "less" && !gaveup) "infinity" else
makefvlabel(NULL, NULL, fname, "hiCI", argname=argname)
mslabl <- c(makefvlabel(NULL, "bar", fname, argname=argname),
makefvlabel("var", "hat", fname, argname=argname),
makefvlabel("res", "hat", fname, argname=argname),
makefvlabel("stdres", "hat", fname, argname=argname),
loCIlabel,
hiCIlabel)
wted <- if(use.weights) "weighted " else NULL
msdesc <- c(paste0(wted, "sample mean of %s from simulations"),
paste0(wted, "sample variance of %s from simulations"),
"raw residual",
"standardised residual",
loCI.name, hiCI.name)
} else {
results <- data.frame(r=argvals,
obs=fX,
mmean=Ef,
lo=lo,
hi=hi)
colnames(results)[1] <- argname
shadenames <- c("lo", "hi")
morestuff <- data.frame(var=varf,
res=fX-Ef,
stdres=stdres,
loCI=loCI,
hiCI=hiCI)
loCIlabel <-
if(alternative == "greater" && !gaveup) "-infinity" else
makefvlabel(NULL, NULL, fname, "loCI", argname=argname)
hiCIlabel <-
if(alternative == "less" && !gaveup) "infinity" else
makefvlabel(NULL, NULL, fname, "hiCI", argname=argname)
mslabl <- c(makefvlabel("var", "hat", fname, argname=argname),
makefvlabel("res", "hat", fname, argname=argname),
makefvlabel("stdres", "hat", fname, argname=argname),
loCIlabel,
hiCIlabel)
msdesc <- c(paste0(if(use.weights) "weighted " else NULL,
"sample variance of %s from simulations"),
"raw residual",
"standardised residual",
loCI.name, hiCI.name)
}
if(do.pwrong) {
## estimate the p-value for the 'wrong test'
if(gaveup) {
pwrong <- NA_real_
} else if(cheat) {
pwrong <- precomputed$pwrong
do.pwrong <- !is.null(pwrong) && !badprobability(pwrong, FALSE)
} else {
upper.signif <- (simvals > hi)
lower.signif <- (simvals < lo)
is.signif <- switch(alternative,
less = lower.signif,
greater = upper.signif,
two.sided = lower.signif | upper.signif)
# is.signif.somewhere <- apply(is.signif, 2, any)
is.signif.somewhere <- matcolany(is.signif)
pwrong <- sum(is.signif.somewhere)/nsim
}
}
}
)
############ WRAP UP #########################
if(use.theory) {
# reference is computed curve `theo'
reflabl <- makefvlabel(NULL, NULL, fname, "theo", argname=argname)
refdesc <- paste0("theoretical value of %s", if(csr) " for CSR" else NULL)
} else {
# reference is sample mean of simulations
reflabl <- makefvlabel(NULL, "bar", fname, argname=argname)
refdesc <- paste0(if(use.weights) "weighted " else NULL,
"sample mean of %s from simulations")
}
lolabl <- if(alternative == "greater" && !gaveup) "-infinity" else
makefvlabel(NULL, "hat", fname, "lo", argname=argname)
hilabl <- if(alternative == "less"&& !gaveup) "infinity" else
makefvlabel(NULL, "hat", fname, "hi", argname=argname)
result <- fv(results,
argu=argname,
ylab=atr$ylab,
valu="obs",
fmla= paste(". ~", argname),
alim=intersect.ranges(atr$alim, range(results[[argname]])),
labl=c(argname,
makefvlabel(NULL, "hat", fname, "obs", argname=argname),
reflabl,
lolabl,
hilabl),
desc=c(arg.desc,
"observed value of %s for data pattern",
refdesc, lo.name, hi.name),
fname=atr$fname,
yexp =atr$yexp)
# columns to be plotted by default
dotty <- c("obs", if(use.theory) "theo" else "mmean", "hi", "lo")
if(type == "variance") {
# add more stuff
result <- bind.fv(result, morestuff, mslabl, msdesc)
if(use.theory) dotty <- c(dotty, "mmean")
}
fvnames(result, ".") <- dotty
fvnames(result, ".s") <- shadenames
unitname(result) <- unitname(funX)
class(result) <- c("envelope", class(result))
# tack on envelope information
attr(result, "einfo") <- list(global = (type =="global"),
ginterval = ginterval,
alternative=alternative,
scale = scale,
clamp = clamp,
csr = csr,
use.theory = use.theory,
csr.theo = csr && use.theory,
simtype = "funs",
constraints = "",
nrank = nrank,
nsim = nsim,
VARIANCE = (type == "variance"),
nSD = nSD,
valname = NULL,
dual = dual,
nsim = nsim,
nsim2 = nsim.mean,
Yname = Yname,
do.pwrong=do.pwrong,
use.weights=use.weights,
gaveup = gaveup)
# tack on saved functions
if(savefuns && !gaveup) {
nSim <- ncol(Y)
alldata <- cbind(argvals, Y)
simnames <- paste("sim", 1:nSim, sep="")
colnames(alldata) <- c(argname, simnames)
alldata <- as.data.frame(alldata)
SimFuns <- fv(alldata,
argu=argname,
ylab=atr$ylab,
valu="sim1",
fmla= paste(". ~", argname),
alim=atr$alim,
labl=names(alldata),
desc=c(arg.desc,
paste("Simulation ", 1:nSim, sep="")),
unitname=unitname(funX))
fvnames(SimFuns, ".") <- simnames
attr(result, "simfuns") <- SimFuns
}
if(do.pwrong)
attr(result, "pwrong") <- pwrong
if(use.weights)
attr(result, "weights") <- weights
return(result)
}
envelope.envelope <- function(Y, fun=NULL, ...,
transform=NULL, global=FALSE, VARIANCE=FALSE) {
Yname <- short.deparse(substitute(Y))
stopifnot(inherits(Y, "envelope"))
Yorig <- Y
aargh <- list(...)
X <- attr(Y, "datapattern")
sf <- attr(Y, "simfuns")
sp <- attr(Y, "simpatterns")
wt <- attr(Y, "weights")
einfo <- attr(Y, "einfo")
csr <- aargh$internal$csr %orifnull% einfo$csr
if(is.null(fun) && is.null(sf)) {
# No simulated functions - must compute them from simulated patterns
if(is.null(sp))
stop(paste("Cannot compute envelope:",
"Y does not contain simulated functions",
"(was not generated with savefuns=TRUE)",
"and does not contain simulated patterns",
"(was not generated with savepatterns=TRUE)"))
# set default fun=Kest
fun <- Kest
}
if(!is.null(fun)) {
# apply new function
# point patterns are required
if(is.null(sp))
stop(paste("Object Y does not contain simulated point patterns",
"(attribute", dQuote("simpatterns"), ");",
"cannot apply a new", sQuote("fun")))
if(is.null(X))
stop(paste("Cannot apply a new", sQuote("fun"),
"; object Y generated by an older version of spatstat"))
## send signal if simulations were CSR
internal <- aargh$internal
if(csr) {
if(is.null(internal)) internal <- list()
internal$csr <- TRUE
}
## compute new envelope
result <- do.call(envelope,
resolve.defaults(list(Y=quote(X), fun=fun, simulate=sp),
aargh,
list(transform=transform,
global=global,
VARIANCE=VARIANCE,
internal=internal,
Yname=Yname,
nsim=einfo$nsim,
nsim2=einfo$nsim2,
weights=wt),
.StripNull=TRUE))
} else {
#' compute new envelope with existing simulated functions
if(is.null(sf))
stop(paste("Y does not contain a", dQuote("simfuns"), "attribute",
"(it was not generated with savefuns=TRUE)"))
if(!is.null(transform)) {
# Apply transformation to Y and sf
stopifnot(is.expression(transform))
## cc <- dotexpr.to.call(transform, "Y", "eval.fv")
cc <- inject.expr("with(Y, .)", transform)
Y <- eval(cc)
## cc <- dotexpr.to.call(transform, "sf", "eval.fv")
cc <- inject.expr("with(sf, .)", transform)
sf <- eval(cc)
}
#' catch discrepancy between domains of observed and simulated functions
if(nrow(sf) != nrow(Y)) {
rrsim <- sf[[fvnames(sf, ".x")]]
rrobs <- Y[[fvnames(Y, ".x")]]
ra <- intersect.ranges(range(rrsim), range(rrobs))
delta <- min(mean(diff(rrsim)), mean(diff(rrobs)))/2
oksim <- (rrsim >= ra[1] - delta) & (rrsim <= ra[2] + delta)
okobs <- (rrobs >= ra[1] - delta) & (rrobs <= ra[2] + delta)
if(sum(oksim) != sum(okobs))
stop("Internal error: Unable to reconcile the domains",
"of the observed and simulated functions", call.=FALSE)
if(mean(abs(rrsim[oksim] - rrobs[okobs])) > delta)
stop("Internal error: Unable to reconcile the r values",
"of the observed and simulated functions", call.=FALSE)
sf <- sf[oksim, ,drop=FALSE]
Y <- Y[okobs, ,drop=FALSE]
}
# extract simulated function values
df <- as.data.frame(sf)
rname <- fvnames(sf, ".x")
df <- df[, (names(df) != rname)]
# interface with 'envelope.matrix'
etype <- if(global) "global" else if(VARIANCE) "variance" else "pointwise"
dfm <- as.matrix(df)
dont.complain.about(dfm)
result <- do.call(envelope.matrix,
resolve.defaults(list(Y=quote(dfm)),
aargh,
list(type=etype,
csr=csr,
funX=Y,
Yname=Yname,
weights=wt),
.StripNull=TRUE))
}
if(!is.null(transform)) {
# post-process labels
labl <- attr(result, "labl")
dnames <- colnames(result)
dnames <- dnames[dnames %in% fvnames(result, ".")]
# expand "."
ud <- as.call(lapply(c("cbind", dnames), as.name))
dont.complain.about(ud)
expandtransform <- eval(substitute(substitute(tr, list(.=ud)),
list(tr=transform[[1L]])))
# compute new labels
attr(result, "fname") <- attr(Yorig, "fname")
mathlabl <- as.character(fvlegend(result, expandtransform))
# match labels to columns
evars <- all.vars(expandtransform)
used.dotnames <- evars[evars %in% dnames]
mathmap <- match(colnames(result), used.dotnames)
okmath <- !is.na(mathmap)
# update appropriate labels
labl[okmath] <- mathlabl[mathmap[okmath]]
attr(result, "labl") <- labl
}
# Tack on envelope info
copyacross <- c("Yname", "csr.theo", "use.theory", "simtype", "constraints")
attr(result, "einfo")[copyacross] <- attr(Yorig, "einfo")[copyacross]
attr(result, "einfo")$csr <- csr
# Save data
return(result)
}
pool.envelope <- local({
pool.envelope <- function(..., savefuns=FALSE, savepatterns=FALSE) {
Yname <- short.deparse(sys.call())
if(nchar(Yname) > 60) Yname <- paste(substr(Yname, 1L, 40L), "[..]")
Elist <- unname(list(...))
nE <- length(Elist)
if(nE == 0) return(NULL)
#' ........ validate envelopes .....................
#' All arguments must be envelopes
notenv <- !unlist(lapply(Elist, inherits, what="envelope"))
if(any(notenv)) {
n <- sum(notenv)
why <- paste(ngettext(n, "Argument", "Arguments"),
commasep(which(notenv)),
ngettext(n, "does not", "do not"),
"belong to the class",
dQuote("envelope"))
stop(why)
}
E1 <- Elist[[1L]]
## Only one envelope?
if(nE == 1)
return(E1)
## envelopes must be compatible
ok <- do.call(compatible, Elist)
if(!ok)
stop("Envelopes are not compatible")
## find name of function argument
argname <- fvnames(E1, ".x")
arg.desc <- attr(E1, "desc")[match(argname, colnames(E1))]
## ... reconcile parameters in different envelopes .......
eilist <- lapply(Elist, attr, which="einfo")
nrank <- resolveEinfo(eilist, "nrank", 1)
global <- resolveEinfo(eilist, "global", FALSE)
ginterval <- resolveEinfo(eilist, "ginterval", NULL, atomic=FALSE)
VARIANCE <- resolveEinfo(eilist, "VARIANCE", FALSE)
alternative <- resolveEinfo(eilist, "alternative", FALSE)
scale <- resolveEinfo(eilist, "scale", NULL, atomic=FALSE)
clamp <- resolveEinfo(eilist, "clamp", FALSE)
resolveEinfo(eilist, "simtype", "funs",
"Envelopes were generated using different types of simulation")
resolveEinfo(eilist, "constraints", "",
"Envelopes were generated using different types of conditioning")
resolveEinfo(eilist, "csr.theo", FALSE, NULL)
csr <- resolveEinfo(eilist, "csr", FALSE, NULL)
use.weights <- resolveEinfo(eilist, "use.weights" , FALSE,
"Weights were used in some, but not all, envelopes: they will be ignored")
use.theory <- resolveEinfo(eilist, "use.theory", csr, NULL)
##
weights <-
if(use.weights) unlist(lapply(Elist, attr, which="weights")) else NULL
type <- if(global) "global" else if(VARIANCE) "variance" else "pointwise"
## ........ validate saved functions .....................
if(savefuns || !VARIANCE) {
## Individual simulated functions are required
SFlist <- lapply(Elist, attr, which="simfuns")
isnul <- unlist(lapply(SFlist, is.null))
if(any(isnul)) {
n <- sum(isnul)
comply <- if(!VARIANCE) "compute the envelope:" else
"save the simulated functions:"
why <- paste("Cannot", comply,
ngettext(n, "argument", "arguments"),
commasep(which(isnul)),
ngettext(n, "does not", "do not"),
"contain a", dQuote("simfuns"), "attribute",
"(not generated with savefuns=TRUE)")
stop(why)
}
## Simulated functions must be the same function
fnames <- unique(lapply(SFlist, attr, which="fname"))
if(length(fnames) > 1L) {
fnames <- unlist(lapply(fnames, flatfname))
stop(paste("Envelope objects contain values",
"of different functions:",
commasep(sQuote(fnames))))
}
## vectors of r values must be identical
rlist <- lapply(SFlist, getargvals)
argvals <- rlist[[1L]]
samer <- unlist(lapply(rlist, identical, y=argvals))
if(!all(samer))
stop(paste("Simulated function values are not compatible",
"(different values of function argument)"))
## Extract function values and assemble into one matrix
matlist <- lapply(SFlist, getdotvals)
SFmatrix <- do.call(cbind, matlist)
}
## compute pooled envelope
switch(type,
pointwise = {
result <- envelope(SFmatrix, funX=E1,
type=type, alternative=alternative,
clamp=clamp,
nrank=nrank,
csr=csr, use.theory=use.theory,
Yname=Yname, weights=weights,
savefuns=savefuns)
},
global = {
simfunmatrix <- if(is.null(ginterval)) SFmatrix else {
## savefuns have not yet been clipped to ginterval
## while envelope data have been clipped.
domain <- (argvals >= ginterval[1L]) & (argvals <= ginterval[2L])
SFmatrix[domain, , drop=FALSE]
}
result <- envelope(simfunmatrix, funX=E1,
type=type, alternative=alternative,
scale=scale, clamp=clamp,
csr=csr, use.theory=use.theory,
nrank=nrank,
ginterval=ginterval,
Yname=Yname, weights=weights,
savefuns=savefuns)
},
variance = {
## Pool sample means and variances
nsims <- unlist(lapply(eilist, getElement, name="nsim"))
mmeans <- lapply(Elist, getElement, name="mmean")
vars <- lapply(Elist, getElement, name="var")
mmeans <- matrix(unlist(mmeans), ncol=nE)
vars <- matrix(unlist(vars), ncol=nE)
if(!use.weights) {
w.mean <- nsims
d.mean <- sum(nsims)
w.var <- nsims - 1
d.var <- sum(nsims) - 1
} else {
weightlist <- lapply(Elist, attr, which="weights")
w.mean <- unlist(lapply(weightlist, sum))
d.mean <- sum(w.mean)
ssw <- unlist(lapply(weightlist, meansqfrac))
## meansqfrac : function(x) {sum((x/sum(x))^2)}))
w.var <- w.mean * (1 - ssw)
d.var <- d.mean * (1 - sum(ssw))
}
poolmmean <- as.numeric(mmeans %*% matrix(w.mean/d.mean, ncol=1L))
within <- vars %*% matrix(w.var, ncol=1L)
between <- ((mmeans - poolmmean[])^2) %*% matrix(w.mean, ncol=1L)
poolvar <- as.numeric((within + between)/d.var)
## feed precomputed data to envelope.matrix
pc <- list(Ef=poolmmean[],
varf=poolvar[])
nsim <- sum(nsims)
result <- envelope.matrix(NULL, funX=E1,
type=type, alternative=alternative,
csr=csr, Yname=Yname,
weights=weights,
savefuns=savefuns,
nsim=nsim,
precomputed=pc)
})
## Copy envelope info that is not handled by envelope.matrix
copyacross <- c("Yname", "csr.theo", "use.theory", "simtype", "constraints")
attr(result, "einfo")[copyacross] <- attr(E1, "einfo")[copyacross]
## ..............saved patterns .....................
if(savepatterns) {
SPlist <- lapply(Elist, attr, which="simpatterns")
isnul <- unlist(lapply(SPlist, is.null))
if(any(isnul)) {
n <- sum(isnul)
why <- paste("Cannot save the simulated patterns:",
ngettext(n, "argument", "arguments"),
commasep(which(isnul)),
ngettext(n, "does not", "do not"),
"contain a", dQuote("simpatterns"), "attribute",
"(not generated with savepatterns=TRUE)")
warning(why)
} else {
attr(result, "simpatterns") <- Reduce(append, SPlist)
}
}
## ..............saved summary functions ................
if(savefuns) {
alldata <- cbind(argvals, SFmatrix)
Nsim <- ncol(SFmatrix)
simnames <- paste0("sim", 1:Nsim)
colnames(alldata) <- c(argname, simnames)
alldata <- as.data.frame(alldata)
SFtemplate <- SFlist[[1L]]
SimFuns <- fv(alldata,
argu=argname,
ylab=attr(SFtemplate, "ylab"),
valu="sim1",
fmla=paste(". ~", argname),
alim=attr(SFtemplate, "alim"),
labl=names(alldata),
desc=c(arg.desc,
paste("Simulation ", 1:Nsim, sep="")),
fname=attr(SFtemplate, "fname"),
yexp=attr(SFtemplate, "yexp"),
unitname=unitname(SFtemplate))
fvnames(SimFuns, ".") <- simnames
attr(result, "simfuns") <- SimFuns
}
dotnames <- lapply(Elist, fvnames, a=".")
dn <- dotnames[[1L]]
if(all(unlist(lapply(dotnames, identical, y=dn))))
fvnames(result, ".") <- dn
shadenames <- lapply(Elist, fvnames, a=".s")
sh <- shadenames[[1L]]
if(all(unlist(lapply(shadenames, identical, y=sh))))
fvnames(result, ".s") <- sh
return(result)
}
getargvals <- function(z) { as.matrix(z)[, fvnames(z, ".x")] }
getdotvals <- function(z) { as.matrix(z)[, fvnames(z, "."), drop=FALSE] }
meansqfrac <- function(x) {sum((x/sum(x))^2)}
pool.envelope
})
# resolve matching entries in different envelope objects
# x is a list of envelope info objects
resolveEinfo <- function(x, what, fallback, warn, atomic=TRUE) {
if(atomic) {
y <- unique(unlist(lapply(x, getElement, name=what)))
if(length(y) == 1L)
return(y)
} else {
y <- unique(lapply(x, getElement, name=what))
if(length(y) == 1L)
return(y[[1L]])
}
if(missing(warn))
warn <- paste("Envelopes were generated using different values",
"of argument", paste(sQuote(what), ";", sep=""),
"reverting to default value")
if(!is.null(warn))
warning(warn, call.=FALSE)
return(fallback)
}
as.data.frame.envelope <- function(x, ..., simfuns=FALSE) {
if(simfuns && !is.null(sf <- attr(x, "simfuns"))) {
# tack on the simulated functions as well
y <- as.data.frame(bind.fv(x, sf, clip=TRUE))
return(y)
}
NextMethod("as.data.frame")
}
Try the spatstat.explore package in your browser
Any scripts or data that you put into this service are public.
spatstat.explore documentation built on Oct. 23, 2023, 1:07 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions as.data.frame.envelope resolveEinfo envelope.envelope envelope.matrix summary.envelope print.envelope plot.envelope envelopeEngine envelope.ppp simulrecipe envelope
Documented in as.data.frame.envelope envelope envelopeEngine envelope.envelope envelope.matrix envelope.ppp plot.envelope print.envelope resolveEinfo simulrecipe summary.envelope
#
# envelope.R
#
# computes simulation envelopes
#
# $Revision: 2.125 $ $Date: 2023/08/15 08:07:52 $
#
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 or Poisson model
# "kppm": simulated realisation of fitted cluster model
# "slrm": simulated realisation of spatial logistic regression
# "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
#
# pois = TRUE if model is known to be Poisson
#
# constraints = additional information about simulation (e.g. 'with fixed n')
#
# ...................................................................
simulrecipe <- function(type, expr, envir, csr, pois=csr, constraints="") {
if(csr && !pois) warning("Internal error: csr=TRUE but pois=FALSE")
out <- list(type=type,
expr=expr,
envir=envir,
csr=csr,
pois=pois,
constraints=constraints)
class(out) <- "simulrecipe"
out
}
envelope.ppp <-
function(Y, fun=Kest, nsim=99, nrank=1, ...,
funargs=list(), funYargs=funargs,
simulate=NULL, fix.n=FALSE, fix.marks=FALSE,
verbose=TRUE, clipdata=TRUE,
transform=NULL, global=FALSE, ginterval=NULL, use.theory=NULL,
alternative=c("two.sided", "less", "greater"),
scale=NULL, clamp=FALSE,
savefuns=FALSE, savepatterns=FALSE, nsim2=nsim,
VARIANCE=FALSE, nSD=2,
Yname=NULL,
maxnerr=nsim, rejectNA=FALSE, silent=FALSE,
do.pwrong=FALSE,
envir.simul=NULL) {
cl <- short.deparse(sys.call())
if(is.null(Yname)) Yname <- short.deparse(substitute(Y))
if(is.null(fun)) fun <- Kest
envir.user <- if(!is.null(envir.simul)) envir.simul else parent.frame()
envir.here <- sys.frame(sys.nframe())
ismarked <- is.marked(Y)
ismulti <- is.multitype(Y)
fix.marks <- fix.marks && ismarked
if(!is.null(simulate)) {
# ...................................................
# Simulations are determined by 'simulate' argument
if(fix.n || fix.marks)
warning("fix.n and fix.marks were ignored, because 'simulate' was given")
# Processing is deferred to envelopeEngine
simrecipe <- simulate
# Data pattern is argument Y
X <- Y
} else if(!fix.n && !fix.marks) {
# ...................................................
# 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
nY <- npoints(Y)
Ywin <- Y$window
Ymarx <- marks(Y)
# expression that will be evaluated
simexpr <- if(is.null(Ymarx)) {
# unmarked point pattern
expression(rpoispp(Yintens, win=Ywin))
} else if(is.null(dim(Ymarx))) {
# single column of marks
expression({
A <- rpoispp(Yintens, win=Ywin);
j <- sample(nY, npoints(A), replace=TRUE);
A %mark% Ymarx[j]
})
} else {
# multiple columns of marks
expression({
A <- rpoispp(Yintens, win=Ywin);
j <- sample(nY, npoints(A), replace=TRUE);
A %mark% Ymarx[j, , drop=FALSE]
})
}
dont.complain.about(Yintens, Ywin)
# evaluate in THIS environment
simrecipe <- simulrecipe(type = "csr",
expr = simexpr,
envir = envir.here,
csr = TRUE,
pois = TRUE)
} else if(fix.marks) {
# ...................................................
# Data pattern is argument Y
X <- Y
# Realisations of binomial process
# with fixed number of points and fixed marks
# will be generated by runifpoint
nY <- npoints(Y)
Ywin <- as.owin(Y)
Ymarx <- marks(Y)
# expression that will be evaluated
simexpr <- expression(runifpoint(nY, Ywin) %mark% Ymarx)
# suppress warnings from code checkers
dont.complain.about(nY, Ywin, Ymarx)
# simulation constraints (explanatory string)
constraints <- if(ismulti) "with fixed number of points of each type" else
"with fixed number of points and fixed marks"
# evaluate in THIS environment
simrecipe <- simulrecipe(type = "csr",
expr = simexpr,
envir = envir.here,
csr = TRUE,
pois = TRUE,
constraints = constraints)
} else {
# ...................................................
# Data pattern is argument Y
X <- Y
# Realisations of binomial process
# will be generated by runifpoint
nY <- npoints(Y)
Ywin <- as.owin(Y)
Ymarx <- marks(Y)
# expression that will be evaluated
simexpr <- if(is.null(Ymarx)) {
## unmarked
expression(runifpoint(nY, Ywin))
} else if(is.null(dim(Ymarx))) {
## single column of marks
expression({
A <- runifpoint(nY, Ywin);
j <- sample(nY, npoints(A), replace=TRUE);
A %mark% Ymarx[j]
})
} else {
## multiple columns of marks
expression({
A <- runifpoint(nY, Ywin);
j <- sample(nY, npoints(A), replace=TRUE);
A %mark% Ymarx[j, ,drop=FALSE]
})
}
dont.complain.about(nY, Ywin)
# evaluate in THIS environment
simrecipe <- simulrecipe(type = "csr",
expr = simexpr,
envir = envir.here,
csr = TRUE,
pois = TRUE,
constraints = "with fixed number of points")
}
envelopeEngine(X=X, fun=fun, simul=simrecipe,
nsim=nsim, nrank=nrank, ...,
funargs=funargs, funYargs=funYargs,
verbose=verbose, clipdata=clipdata,
transform=transform,
global=global, ginterval=ginterval, use.theory=use.theory,
alternative=alternative, scale=scale, clamp=clamp,
savefuns=savefuns, savepatterns=savepatterns, nsim2=nsim2,
VARIANCE=VARIANCE, nSD=nSD,
Yname=Yname,
maxnerr=maxnerr, rejectNA=rejectNA, silent=silent,
cl=cl, envir.user=envir.user, do.pwrong=do.pwrong)
}
## Code for envelope.ppm, envelope.kppm, envelope.slrm
## is moved to spatstat.model
## .................................................................
## 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, ..., funargs=list(), funYargs=funargs,
verbose=TRUE, clipdata=TRUE,
transform=NULL, global=FALSE, ginterval=NULL, use.theory=NULL,
alternative=c("two.sided", "less", "greater"),
scale=NULL, clamp=FALSE,
savefuns=FALSE, savepatterns=FALSE,
saveresultof=NULL,
weights=NULL,
nsim2=nsim,
VARIANCE=FALSE, nSD=2,
Yname=NULL,
maxnerr=nsim,
rejectNA=FALSE,
silent=FALSE,
maxerr.action=c("fatal", "warn", "null"),
internal=NULL, cl=NULL,
envir.user=envir.user,
expected.arg="r",
do.pwrong=FALSE,
foreignclass=NULL,
collectrubbish=FALSE)
{
#
envir.here <- sys.frame(sys.nframe())
alternative <- match.arg(alternative)
maxerr.action <- match.arg(maxerr.action)
foreignclass <- as.character(foreignclass)
if(length(foreignclass) != 0 && clipdata) {
warning(paste("Ignoring clipdata=TRUE:",
"I don't know how to clip objects of class",
sQuote(paste(foreignclass, collapse=","))))
clipdata <- FALSE
}
# ----------------------------------------------------------
# Determine Simulation
# ----------------------------------------------------------
check.1.integer(nsim)
stopifnot(nsim > 1)
# 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
if(inherits(X, foreignclass)) foreignclass else
stop("Unrecognised class of point pattern")
Xobjectname <- paste("point pattern of class", sQuote(Xclass))
# Option to use weighted average
if(use.weights <- !is.null(weights)) {
# weight can be either a numeric vector or a function
if(is.numeric(weights)) {
compute.weights <- FALSE
weightfun <- NULL
} else if(is.function(weights)) {
compute.weights <- TRUE
weightfun <- weights
weights <- NULL
} else stop("weights should be either a function or a numeric vector")
} else compute.weights <- FALSE
# Undocumented option to generate patterns only.
patterns.only <- identical(internal$eject, "patterns")
# Undocumented option to evaluate 'something' for each simulation
if(savevalues <- !is.null(saveresultof)) {
stopifnot(is.function(saveresultof))
SavedValues <- list()
## might be a function of the pattern only, or both pattern and summary fun
result.depends.both <- (length(formals(saveresultof)) >= 2)
}
# 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
pois <- simul$pois
constraints <- simul$constraints
} 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, or a function to be applied to the data
csr <- FALSE
# override
if(!is.null(icsr <- internal$csr)) csr <- icsr
pois <- csr
constraints <- ""
# 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.function(simulate)) {
## User-supplied function 'simulate' will be repeatedly evaluated on X
simtype <- "func"
simexpr <- expression(simulate(X))
envir <- envir.here
} else if(is.list(simulate) &&
all(sapply(simulate, inherits, what=Xclass))) {
#' The user-supplied list of point patterns will be used
simtype <- "list"
SimDataList <- simulate
#' expression that will be evaluated
simexpr <- expression(SimDataList[[i+nerr]])
dont.complain.about(SimDataList)
envir <- envir.here
#' ensure that `i' is defined
i <- 1L
nerr <- 0L
maxnerr <- min(length(SimDataList)-nsim, maxnerr)
#' any messages?
if(!is.null(mess <- attr(simulate, "internal"))) {
# determine whether these point patterns are realisations of CSR
csr <- isTRUE(mess$csr)
}
} else if(is.list(simulate) &&
all(sapply(simulate, is.list)) &&
all(lengths(simulate) == 1) &&
all(sapply((elements <- lapply(simulate, "[[", i=1)), inherits, what=Xclass))) {
#' malformed argument: list(list(ppp), list(ppp), ....)
SimDataList <- elements
simtype <- "list"
#' expression that will be evaluated
simexpr <- expression(SimDataList[[i+nerr]])
dont.complain.about(SimDataList)
envir <- envir.here
#' ensure that `i' is defined
i <- 1L
nerr <- 0L
maxnerr <- min(length(SimDataList)-nsim, maxnerr)
#' any messages?
if(!is.null(mess <- attr(simulate, "internal"))) {
# determine whether these point patterns are realisations of CSR
csr <- isTRUE(mess$csr)
}
} else stop(paste(sQuote("simulate"),
"should be an expression,",
"or a list of point patterns of the same kind as X"))
}
# -------------------------------------------------------------------
# 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)),
slrm=stop(paste("Internal error: simulate.slrm did not return an",
Xobjectname)),
expr=stop(paste("Evaluating the expression", sQuote("simulate"),
"did not yield an", Xobjectname)),
func=stop(paste("Evaluating the function", 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")
}
# ------------------------------------------------------------------
# Summary function to be applied
# ------------------------------------------------------------------
if(is.null(fun))
stop("Internal error: fun is NULL")
# Name of function, for error messages
fname <- if(is.name(substitute(fun))) short.deparse(substitute(fun)) else
if(is.character(fun)) fun else "fun"
fname <- sQuote(fname)
# R function to apply
if(is.character(fun)) {
gotfun <- try(get(fun, mode="function"))
if(inherits(gotfun, "try-error"))
stop(paste("Could not find a function named", sQuote(fun)))
fun <- gotfun
} else if(!is.function(fun))
stop(paste("unrecognised format for function", fname))
fargs <- names(formals(fun))
if(!any(c(expected.arg, "...") %in% fargs))
stop(paste(fname, "should have",
ngettext(length(expected.arg), "an argument", "arguments"),
"named", commasep(sQuote(expected.arg)),
"or a", sQuote("..."), "argument"))
usecorrection <- any(c("correction", "...") %in% fargs)
# ---------------------------------------------------------------------
# validate other arguments
if((nrank %% 1) != 0)
stop("nrank must be an integer")
if((nsim %% 1) != 0)
stop("nsim must be an integer")
stopifnot(nrank > 0 && nrank < nsim/2)
arg.given <- any(expected.arg %in% names(list(...)))
if(tran <- !is.null(transform)) {
stopifnot(is.expression(transform))
# prepare expressions to be evaluated each time
transform.funX <- inject.expr("with(funX,.)", transform)
transform.funXsim <- inject.expr("with(funXsim,.)", transform)
# .... old code using 'eval.fv' ......
# transform.funX <- dotexpr.to.call(transform, "funX", "eval.fv")
# transform.funXsim <- dotexpr.to.call(transform, "funXsim", "eval.fv")
# 'transform.funX' and 'transform.funXsim' are unevaluated calls to eval.fv
}
if(!is.null(ginterval))
stopifnot(is.numeric(ginterval) && length(ginterval) == 2)
# ---------------------------------------------------------------------
# Evaluate function for data pattern X
# ------------------------------------------------------------------
Xarg <- if(!clipdata) X else X[clipwin]
corrx <- if(usecorrection) list(correction="best") else NULL
dont.complain.about(Xarg)
funX <- do.call(fun,
resolve.defaults(list(quote(Xarg)),
list(...),
funYargs,
corrx))
if(!inherits(funX, "fv"))
stop(paste("The function", fname,
"must return an object of class", sQuote("fv")))
## catch 'conservation' parameters
conserveargs <- attr(funX, "conserve")
if(!is.null(conserveargs) && !any(c("conserve", "...") %in% fargs))
stop(paste("In this usage, the function", fname,
"should have an argument named 'conserve' or '...'"))
## warn about 'dangerous' arguments
if(!is.null(dang <- attr(funX, "dangerous")) &&
any(uhoh <- dang %in% names(list(...)))) {
nuh <- sum(uhoh)
warning(paste("Envelope may be invalid;",
ngettext(nuh, "argument", "arguments"),
commasep(sQuote(dang[uhoh])),
ngettext(nuh, "appears", "appear"),
"to have been fixed."),
call.=FALSE)
}
argname <- fvnames(funX, ".x")
valname <- fvnames(funX, ".y")
has.theo <- "theo" %in% fvnames(funX, "*")
csr.theo <- csr && has.theo
use.theory <- if(is.null(use.theory)) csr.theo else (use.theory && has.theo)
if(tran) {
# extract only the recommended value
if(use.theory)
funX <- funX[, c(argname, valname, "theo")]
else
funX <- funX[, c(argname, valname)]
# apply the transformation to it
funX <- eval(transform.funX)
}
argvals <- funX[[argname]]
# fX <- funX[[valname]]
arg.desc <- attr(funX, "desc")[match(argname, colnames(funX))]
# default domain over which to maximise
alim <- attr(funX, "alim")
if(global && is.null(ginterval))
ginterval <- if(arg.given || is.null(alim)) range(argvals) else alim
#--------------------------------------------------------------------
# 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 && !use.theory && !VARIANCE)
if(dual) {
check.1.integer(nsim2)
stopifnot(nsim2 >= 1)
}
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(patterns.only) {
# 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 = paste("simulated realisations of fitted",
if(pois) "Poisson" else "Gibbs",
"model"),
kppm = "simulated realisations of fitted cluster model",
slrm = "simulated realisations of spatial logistic regression model",
expr = "simulations by evaluating expression",
func = "simulations by evaluating function",
list = "point patterns from list",
"simulated realisations")
if(!is.null(constraints) && nzchar(constraints))
descrip <- paste(descrip, constraints)
explan <- if(dual) paren(paste(nsim2, "to estimate the mean and",
nsim, "to calculate envelopes")) else ""
splat(action, Nsim, descrip, explan, "...")
}
XsimList <- list()
# start simulation loop
sstate <- list()
for(i in 1:Nsim) {
if(verbose) sstate <- progressreport(i, Nsim, state=sstate)
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))
},
slrm={
stop(paste("Internal error: simulate.slrm did not return an",
Xobjectname))
},
expr={
stop(paste("Evaluating the expression", sQuote("simulate"),
"did not yield an", Xobjectname))
},
func={
stop(paste("Evaluating the function", 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)
}
# capture main decision parameters
envelopeInfo <- list(call=cl,
Yname=Yname,
valname=valname,
csr=csr,
csr.theo=csr.theo,
use.theory=use.theory,
pois=pois,
simtype=simtype,
constraints=constraints,
nrank=nrank,
nsim=nsim,
Nsim=Nsim,
global=global,
ginterval=ginterval,
dual=dual,
nsim2=nsim2,
VARIANCE=VARIANCE,
nSD=nSD,
alternative=alternative,
scale=scale,
clamp=clamp,
use.weights=use.weights,
do.pwrong=do.pwrong)
# ----------------------------------------
######### SIMULATE #######################
# ----------------------------------------
if(verbose) {
action <- if(simtype == "list") "Extracting" else "Generating"
descrip <- switch(simtype,
csr = "simulations of CSR",
rmh = paste("simulated realisations of fitted",
if(pois) "Poisson" else "Gibbs",
"model"),
kppm = "simulated realisations of fitted cluster model",
slrm = "simulated realisations of fitted spatial logistic regression model",
expr = "simulations by evaluating expression",
func = "simulations by evaluating function",
list = "point patterns from list",
"simulated patterns")
if(!is.null(constraints) && nzchar(constraints))
descrip <- paste(descrip, constraints)
explan <- if(dual) paren(paste(nsim2, "to estimate the mean and",
nsim, "to calculate envelopes")) else ""
splat(action, Nsim, descrip, explan, "...")
}
# determine whether simulated point patterns should be saved
catchpatterns <- savepatterns && simtype != "list"
Caughtpatterns <- list()
# allocate space for computed function values
nargvals <- length(argvals)
simvals <- matrix(, nrow=nargvals, ncol=Nsim)
# allocate space for weights to be computed
if(compute.weights)
weights <- numeric(Nsim)
# inferred values of function argument 'r' or equivalent parameters
if(identical(expected.arg, "r") && identical(argname, "r")) {
# Kest, etc
inferred.r.args <- list(r=argvals)
} else if(identical(expected.arg, c("rmax", "nrval"))) {
# K3est, etc
inferred.r.args <- list(rmax=max(argvals), nrval=length(argvals))
} else if(any(c(argname, "...") %in% fargs)) {
## assume it accepts the vector of argument values
inferred.r.args <- structure(list(argvals), names=argname)
} else {
stop(paste("Don't know how to infer values of",
if(length(expected.arg))
commasep(sQuote(expected.arg)) else "function argument"))
}
# arguments for function when applied to simulated patterns
funargs <-
resolve.defaults(funargs,
inferred.r.args,
list(...),
conserveargs,
if(usecorrection) list(correction="best") else NULL)
# reject simulated pattern if function values are all NA (etc)
rejectNA <- isTRUE(rejectNA)
# start simulation loop
nerr <- 0
gaveup <- FALSE
if(verbose) pstate <- list()
for(i in 1:Nsim) {
## safely generate a random pattern and apply function
success <- FALSE
while(!success && !gaveup) {
Xsim <- eval(simexpr, envir=envir)
## check valid point pattern
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)),
slrm=stop(paste("Internal error:",
"simulate.slrm did not return an",
Xobjectname)),
expr=stop(paste("Evaluating the expression", sQuote("simulate"),
"did not yield an", Xobjectname)),
func=stop(paste("Evaluating the function", 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
if(savevalues && !result.depends.both)
SavedValues[[i]] <- saveresultof(Xsim)
if(compute.weights) {
wti <- weightfun(Xsim)
if(!is.numeric(wti))
stop("weightfun did not return a numeric value")
if(length(wti) != 1L)
stop("weightfun should return a single numeric value")
weights[i] <- wti
}
## apply function safely
funXsim <- try(do.call(fun, c(list(Xsim), funargs)), silent=silent)
success <-
!inherits(funXsim, "try-error") &&
inherits(funXsim, "fv") &&
(!rejectNA || any(is.finite(funXsim[[valname]])))
if(!success) {
#' error in computing summary function
nerr <- nerr + 1L
if(nerr > maxnerr) {
gaveup <- TRUE
errtype <- if(rejectNA) "fatal errors or NA function values"
if(simtype == "list") {
whinge <- paste("Exceeded maximum possible number of errors",
"when evaluating summary function:",
length(SimDataList), "patterns provided,",
nsim, "patterns required,",
nerr, ngettext(nerr, "pattern", "pattern"),
"rejected due to", errtype)
} else {
whinge <- paste("Exceeded maximum permissible number of",
errtype,
paren(paste("maxnerr =", maxnerr)),
"when evaluating summary function",
"for simulated point patterns")
}
switch(maxerr.action,
fatal = stop(whinge, call.=FALSE),
warn = warning(whinge, call.=FALSE),
null = {})
} else if(!silent) cat("[retrying]\n")
}
#' ..... end of while(!success) ................
}
if(gaveup) break; # exit loop now
## sanity checks
if(i == 1L) {
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"))
rfunX <- with(funX, ".x")
rfunXsim <- with(funXsim, ".x")
if(!identical(rfunX, rfunXsim))
stop(paste("When", fname, "is applied to a simulated pattern,",
"the values of the argument", sQuote(argname.sim),
"are different from those used for the data."))
}
if(savevalues && result.depends.both)
SavedValues[[i]] <- saveresultof(Xsim, funXsim)
if(tran) {
# extract only the recommended value
if(use.theory)
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) != nargvals)
stop("Vectors of function values have incompatible lengths")
simvals[ , i] <- funXsim[[valname]]
if(verbose)
pstate <- progressreport(i, Nsim, state=pstate)
if(collectrubbish) {
rm(Xsim)
rm(funXsim)
gc()
}
}
## end simulation loop
if(verbose) {
cat("\nDone.\n")
flush.console()
}
# ...........................................................
# save functions and/or patterns if so commanded
if(!gaveup) {
if(savefuns) {
alldata <- cbind(argvals, simvals)
simnames <- paste("sim", 1:Nsim, sep="")
colnames(alldata) <- c(argname, simnames)
alldata <- as.data.frame(alldata)
SimFuns <- fv(alldata,
argu=argname,
ylab=attr(funX, "ylab"),
valu="sim1",
fmla= paste(". ~", argname),
alim=attr(funX, "alim"),
labl=names(alldata),
desc=c(arg.desc,
paste("Simulation ", 1:Nsim, sep="")),
fname=attr(funX, "fname"),
yexp=attr(funX, "yexp"),
unitname=unitname(funX))
fvnames(SimFuns, ".") <- simnames
}
if(savepatterns)
SimPats <- if(simtype == "list") SimDataList else Caughtpatterns
}
######### COMPUTE ENVELOPES #######################
etype <- if(global) "global" else if(VARIANCE) "variance" else "pointwise"
if(dual) {
jsim <- 1:nsim
jsim.mean <- nsim + 1:nsim2
} else {
jsim <- jsim.mean <- NULL
}
result <- envelope.matrix(simvals, funX=funX,
jsim=jsim, jsim.mean=jsim.mean,
type=etype, alternative=alternative,
scale=scale, clamp=clamp,
csr=csr, use.theory=use.theory,
nrank=nrank, ginterval=ginterval, nSD=nSD,
Yname=Yname, do.pwrong=do.pwrong,
weights=weights, gaveup=gaveup)
## tack on envelope information
attr(result, "einfo") <- resolve.defaults(envelopeInfo,
attr(result, "einfo"))
if(!gaveup) {
## tack on functions and/or patterns if so commanded
if(savefuns)
attr(result, "simfuns") <- SimFuns
if(savepatterns) {
attr(result, "simpatterns") <- SimPats
attr(result, "datapattern") <- X
}
## undocumented - tack on values of some other quantity
if(savevalues) {
attr(result, "simvalues") <- SavedValues
attr(result, "datavalue") <-
if(result.depends.both) saveresultof(X, funX) else saveresultof(X)
}
}
## save function weights
if(use.weights)
attr(result, "weights") <- weights
return(result)
}
plot.envelope <- function(x, ..., main) {
if(missing(main)) main <- short.deparse(substitute(x))
shade.given <- ("shade" %in% names(list(...)))
shade.implied <- !is.null(fvnames(x, ".s"))
if(!(shade.given || shade.implied)) {
# ensure x has default 'shade' attribute
# (in case x was produced by an older version of spatstat)
if(all(c("lo", "hi") %in% colnames(x)))
fvnames(x, ".s") <- c("lo", "hi")
else warning("Unable to determine shading for envelope")
}
NextMethod("plot", main=main)
}
print.envelope <- function(x, ...) {
e <- attr(x, "einfo")
g <- e$global
csr <- e$csr
pois <- e$pois
if(is.null(pois)) pois <- csr
simtype <- e$simtype
constraints <- e$constraints
nr <- e$nrank
nsim <- e$nsim
V <- e$VARIANCE
fname <- flat.deparse(attr(x, "ylab"))
type <- if(V) paste("Pointwise", e$nSD, "sigma") else
if(g) "Simultaneous" else "Pointwise"
splat(type, "critical envelopes for", fname,
"\nand observed value for", sQuote(e$Yname))
if(!is.null(valname <- e$valname) && waxlyrical('extras'))
splat("Edge correction:", dQuote(valname))
## determine *actual* type of simulation
descrip <-
if(csr) "simulations of CSR"
else if(!is.null(simtype)) {
switch(simtype,
csr="simulations of CSR",
rmh=paste("simulations of fitted",
if(pois) "Poisson" else "Gibbs",
"model"),
kppm="simulations of fitted cluster model",
slrm="simulations of fitted spatial logistic regression model",
expr="evaluations of user-supplied expression",
func="evaluations of user-supplied function",
list="point pattern datasets in user-supplied list",
funs="columns of user-supplied data")
} else "simulations of fitted model"
if(!is.null(constraints) && nzchar(constraints))
descrip <- paste(descrip, constraints)
#
splat("Obtained from", nsim, descrip)
#
if(waxlyrical('extras')) {
dual <- isTRUE(e$dual)
usetheory <- isTRUE(e$use.theory)
hownull <- if(usetheory) {
"(known exactly)"
} else if(dual) {
paste("(estimated from a separate set of",
e$nsim2, "simulations)")
} else NULL
formodel <- if(csr) "for CSR" else NULL
if(g) {
splat("Envelope based on maximum deviation of", fname,
"from null value", formodel, hownull)
} else if(dual) {
splat("Null value of", fname, formodel, hownull)
}
if(!is.null(attr(x, "simfuns")))
splat("(All simulated function values are stored)")
if(!is.null(attr(x, "simpatterns")))
splat("(All simulated point patterns are stored)")
}
splat("Alternative:", e$alternative)
if(!V && waxlyrical('extras')) {
## significance interpretation!
alpha <- if(g) { nr/(nsim+1) } else { 2 * nr/(nsim+1) }
splat("Significance level of",
if(g) "simultaneous" else "pointwise",
"Monte Carlo test:",
paste0(if(g) nr else 2 * nr,
"/", nsim+1),
"=", signif(alpha, 3))
}
if(waxlyrical('gory') && !is.null(pwrong <- attr(x, "pwrong"))) {
splat("\t[Estimated significance level of pointwise excursions:",
paste0("pwrong=", signif(pwrong, 3), "]"))
}
NextMethod("print")
}
summary.envelope <- function(object, ...) {
e <- attr(object, "einfo")
g <- e$global
V <- e$VARIANCE
nr <- e$nrank
nsim <- e$nsim
csr <- e$csr
pois <- e$pois
if(is.null(pois)) pois <- csr
simtype <- e$simtype
constraints <- e$constraints
alternative <- e$alternative
use.theory <- e$use.theory
has.theo <- "theo" %in% fvnames(object, "*")
csr.theo <- csr && has.theo
use.theory <- if(is.null(use.theory)) csr.theo else (use.theory && has.theo)
fname <- short.deparse(attr(object, "ylab"))
type <- if(V) paste("Pointwise", e$nSD, "sigma") else
if(g) "Simultaneous" else "Pointwise"
splat(type, "critical envelopes for", fname,
"\nand observed value for", sQuote(e$Yname))
# determine *actual* type of simulation
descrip <-
if(csr) "simulations of CSR"
else if(!is.null(simtype)) {
switch(simtype,
csr="simulations of CSR",
rmh=paste("simulations of fitted",
if(pois) "Poisson" else "Gibbs",
"model"),
kppm="simulations of fitted cluster model",
slrm="simulations of fitted spatial logistic regression model",
expr="evaluations of user-supplied expression",
func="evaluations of user-supplied function",
list="point pattern datasets in user-supplied list",
funs="columns of user-supplied data",
"simulated point patterns")
} else "simulations of fitted model"
if(!is.null(constraints) && nzchar(constraints))
descrip <- paste(descrip, constraints)
#
splat("Obtained from", nsim, descrip)
#
if(waxlyrical('extras')) {
if(!is.null(e$dual) && e$dual)
splat("Theoretical (i.e. null) mean value of", fname,
"estimated from a separate set of",
e$nsim2, "simulations")
if(!is.null(attr(object, "simfuns")))
splat("(All", nsim, "simulated function values",
"are stored in attr(,", dQuote("simfuns"), ") )")
if(!is.null(attr(object, "simpatterns")))
splat("(All", nsim, "simulated point patterns",
"are stored in attr(,", dQuote("simpatterns"), ") )")
}
#
splat("Alternative:", alternative)
if(V) {
# nSD envelopes
splat(switch(alternative,
two.sided = "Envelopes",
"Critical boundary"),
"computed as sample mean",
switch(alternative,
two.sided="plus/minus",
less="minus",
greater="plus"),
e$nSD, "sample standard deviations")
} else {
# critical envelopes
lo.ord <- if(nr == 1L) "minimum" else paste(ordinal(nr), "smallest")
hi.ord <- if(nr == 1L) "maximum" else paste(ordinal(nr), "largest")
if(g)
splat(switch(alternative,
two.sided = "Envelopes",
"Critical boundary"),
"computed as",
if(use.theory) "theoretical curve" else "mean of simulations",
switch(alternative,
two.sided="plus/minus",
less="minus",
greater="plus"),
hi.ord,
"simulated value of maximum",
switch(alternative,
two.sided="absolute",
less="negative",
greater="positive"),
"deviation")
else {
if(alternative != "less")
splat("Upper envelope: pointwise", hi.ord, "of simulated curves")
if(alternative != "greater")
splat("Lower envelope: pointwise", lo.ord, "of simulated curves")
}
symmetric <- (alternative == "two.sided") && !g
alpha <- if(!symmetric) { nr/(nsim+1) } else { 2 * nr/(nsim+1) }
splat("Significance level of Monte Carlo test:",
paste0(if(!symmetric) nr else 2 * nr,
"/", nsim+1),
"=", alpha)
}
splat("Data:", e$Yname)
return(invisible(NULL))
}
# envelope.matrix
# core functionality to compute envelope values
# theory = funX[["theo"]]
# observed = fX
envelope.matrix <- function(Y, ...,
argvals=rvals, rvals=NULL, ## rvals is old name
observed=NULL, theory=NULL,
funX=NULL,
nsim=NULL, nsim2=NULL,
jsim=NULL, jsim.mean=NULL,
type=c("pointwise", "global", "variance"),
alternative=c("two.sided", "less", "greater"),
scale = NULL, clamp=FALSE,
csr=FALSE, use.theory = csr,
nrank=1, ginterval=NULL, nSD=2,
savefuns=FALSE,
check=TRUE,
Yname=NULL,
argname=NULL,
arg.desc=NULL,
do.pwrong=FALSE,
weights=NULL,
precomputed=NULL,
gaveup=FALSE) {
if(is.null(Yname))
Yname <- short.deparse(substitute(Y))
type <- match.arg(type)
alternative <- match.arg(alternative)
if(!is.null(funX))
stopifnot(is.fv(funX))
pwrong <- NULL
use.weights <- !is.null(weights)
cheat <- !is.null(precomputed)
if(is.null(argvals) && is.null(observed) && !is.null(funX)) {
## assume funX is summary function for observed data
argvals <- with(funX, .x)
observed <- with(funX, .y)
argname.orig <- fvnames(funX, ".x")
if(is.null(argname)) argname <- argname.orig
if(is.null(arg.desc))
arg.desc <- attr(funX, "desc")[match(argname.orig, colnames(funX))]
theory <- if(use.theory) (theory %orifnull% funX[["theo"]]) else NULL
if(check) stopifnot(nrow(funX) == nrow(Y))
} else {
## construct envelope from raw data
if(is.null(argname)) argname <- "r"
if(is.null(arg.desc)) arg.desc <- "distance argument r"
if(check) {
## validate vectors of data
if(is.null(argvals)) stop("argvals must be supplied")
if(is.null(observed)) stop("observed must be supplied")
stopifnot(length(argvals) == nrow(Y))
stopifnot(length(observed) == length(argvals))
}
}
use.theory <- use.theory && !is.null(theory)
if(use.theory && check) stopifnot(length(theory) == length(argvals))
simvals <- Y
fX <- observed
if(!is.null(funX)) {
atr <- attributes(funX)
fname <- atr$fname
yexp <- atr$yexp
} else {
fname <- "f"
yexp <- substitute(f(r), list(r=as.name(argname)))
atr <- list(alim=range(argvals),
ylab=yexp,
yexp=yexp,
fname="f")
}
NAvector <- rep(NA_real_, length(argvals))
if(!cheat) {
## ................ standard calculation .....................
## validate weights
if(use.weights && !gaveup)
check.nvector(weights, ncol(simvals),
things="simulated functions", naok=TRUE, vname="weights")
## determine numbers of columns used
Ncol <- if(!gaveup) ncol(simvals) else Inf
if(Ncol < 2)
stop("Need at least 2 columns of function values")
## all columns are used unless 'nsim' or 'jsim' given.
if(!(is.null(nsim) && is.null(jsim))) {
if(is.null(jsim)) {
jsim <- 1:nsim
} else if(is.null(nsim)) {
nsim <- length(jsim)
} else stopifnot(length(jsim) == nsim)
if(nsim > Ncol)
stop(paste(nsim, "simulations are not available; only",
Ncol, "columns provided"))
}
## nsim2 or jsim.mean may be given, and imply dual calculation
if(!(is.null(nsim2) && is.null(jsim.mean))) {
if(is.null(jsim.mean)) {
jsim.mean <- setdiff(seq_len(Ncol), jsim)[1:nsim2]
} else if(is.null(nsim2)) {
nsim2 <- length(jsim.mean)
} else stopifnot(length(jsim.mean) == nsim2)
if(nsim + nsim2 > Ncol)
stop(paste(nsim, "+", nsim2, "=", nsim+nsim2,
"simulations are not available; only",
Ncol, "columns provided"))
if(length(intersect(jsim, jsim.mean)))
warning("Internal warning: Indices in jsim and jsim.mean overlap")
}
restrict.columns <- !is.null(jsim)
dual <- !is.null(jsim.mean)
} else {
## ................ precomputed values ..................
## validate weights
if(use.weights)
check.nvector(weights, nsim,
things="simulations", naok=TRUE, vname="weights")
restrict.columns <- FALSE
dual <- FALSE
}
shadenames <- NULL
nsim.mean <- NULL
switch(type,
pointwise = {
## ....... POINTWISE ENVELOPES ...............................
if(gaveup) {
lo <- hi <- NAvector
} else if(cheat) {
stopifnot(checkfields(precomputed, c("lo", "hi")))
lo <- precomputed$lo
hi <- precomputed$hi
} else {
simvals[is.infinite(simvals)] <- NA
if(restrict.columns) {
simvals <- simvals[,jsim]
if(use.weights) weights <- weights[jsim]
}
nsim <- ncol(simvals)
if(nrank == 1L) {
lohi <- apply(simvals, 1L, range)
} else {
lohi <- apply(simvals, 1L,
# function(x, n) { sort(x)[n] },
orderstats,
k=c(nrank, nsim-nrank+1L))
}
lo <- lohi[1L,]
hi <- lohi[2L,]
}
lo.name <- "lower pointwise envelope of %s from simulations"
hi.name <- "upper pointwise envelope of %s from simulations"
##
if(!gaveup)
switch(alternative,
two.sided = { },
less = {
hi <- rep.int(Inf, length(hi))
hi.name <- "infinite upper limit"
},
greater = {
lo <- rep.int(-Inf, length(lo))
lo.name <- "infinite lower limit"
})
if(use.theory) {
results <- data.frame(r=argvals,
obs=fX,
theo=theory,
lo=lo,
hi=hi)
} else {
m <- if(gaveup) NAvector else
if(cheat) precomputed$mmean else
if(!use.weights) apply(simvals, 1L, mean, na.rm=TRUE) else
apply(simvals, 1L, weighted.mean, w=weights, na.rm=TRUE)
results <- data.frame(r=argvals,
obs=fX,
mmean=m,
lo=lo,
hi=hi)
}
colnames(results)[1] <- argname
shadenames <- c("lo", "hi")
if(do.pwrong) {
## estimate the p-value for the 'wrong test'
if(gaveup) {
pwrong <- NA_real_
} else if(cheat) {
pwrong <- precomputed$pwrong
do.pwrong <- !is.null(pwrong) && !badprobability(pwrong, FALSE)
} else {
dataranks <- t(apply(simvals, 1, rank, ties.method="random"))
upper.signif <- (dataranks <= nrank)
lower.signif <- (dataranks >= nsim-nrank+1L)
is.signif <- switch(alternative,
less = lower.signif,
greater = upper.signif,
two.sided = lower.signif | upper.signif)
is.signif.somewhere <- matcolany(is.signif)
pwrong <- sum(is.signif.somewhere)/nsim
}
}
},
global = {
## ..... SIMULTANEOUS ENVELOPES ..........................
if(gaveup) {
lo <- hi <- reference <- NAvector
} else if(cheat) {
## ... use precomputed values ..
stopifnot(checkfields(precomputed, c("lo", "hi")))
lo <- precomputed$lo
hi <- precomputed$hi
if(use.theory) {
reference <- theory
} else {
stopifnot(checkfields(precomputed, "mmean"))
reference <- precomputed$mmean
}
domain <- rep.int(TRUE, length(argvals))
} else {
## ... normal case: compute envelopes from simulations
if(!is.null(ginterval)) {
domain <- (argvals >= ginterval[1L]) & (argvals <= ginterval[2L])
funX <- funX[domain, ]
simvals <- simvals[domain, ]
} else domain <- rep.int(TRUE, length(argvals))
simvals[is.infinite(simvals)] <- NA
if(use.theory) {
reference <- theory[domain]
if(restrict.columns) {
simvals <- simvals[, jsim]
if(use.weights) weights <- weights[jsim]
}
} else if(dual) {
# Estimate the mean from one set of columns
# Form envelopes from another set of columns
simvals.mean <- simvals[, jsim.mean]
# mmean <-
reference <-
if(!use.weights) apply(simvals.mean, 1L, mean, na.rm=TRUE) else
apply(simvals.mean, 1L, weighted.mean, w=weights[jsim.mean],
na.rm=TRUE)
nsim.mean <- ncol(simvals.mean)
# retain only columns used for envelope
simvals <- simvals[, jsim]
} else {
# Compute the mean and envelopes using the same data
if(restrict.columns) {
simvals <- simvals[, jsim]
if(use.weights) weights <- weights[jsim]
}
# mmean <-
reference <-
if(!use.weights) apply(simvals, 1L, mean, na.rm=TRUE) else
apply(simvals, 1L, weighted.mean, w=weights, na.rm=TRUE)
}
nsim <- ncol(simvals)
# compute deviations
deviations <- sweep(simvals, 1L, reference)
deviations <-
switch(alternative,
two.sided = abs(deviations),
greater = if(clamp) pmax(0, deviations) else deviations,
less = if(clamp) pmax(0, -deviations) else (-deviations))
deviations <- matrix(deviations,
nrow=nrow(simvals), ncol=ncol(simvals))
## rescale ?
sc <- 1
if(!is.null(scale)) {
stopifnot(is.function(scale))
sc <- scale(argvals)
sname <- paste0("scale", paren(argname))
tname <- paste("values of", argname)
ans <- check.nvector(sc, length(argvals), things=tname,
fatal=FALSE, vname=sname)
if(!ans)
stop(attr(ans, "whinge"), call.=FALSE)
if(any(bad <- (sc <= 0))) {
## issue a warning unless this only happens at r=0
if(any(bad[argvals > 0]))
warning(paste("Some values of", sname,
"were negative or zero:",
"scale was reset to 1 for these values"),
call.=FALSE)
sc[bad] <- 1
}
deviations <- sweep(deviations, 1L, sc, "/")
}
## compute max (scaled) deviations
suprema <- apply(deviations, 2L, max, na.rm=TRUE)
# ranked deviations
dmax <- sort(suprema)[nsim-nrank+1L]
# simultaneous bands
lo <- reference - sc * dmax
hi <- reference + sc * dmax
}
lo.name <- "lower critical boundary for %s"
hi.name <- "upper critical boundary for %s"
if(!gaveup)
switch(alternative,
two.sided = { },
less = {
hi <- rep.int(Inf, length(hi))
hi.name <- "infinite upper boundary"
},
greater = {
lo <- rep.int(-Inf, length(lo))
lo.name <- "infinite lower boundary"
})
if(use.theory) {
results <- data.frame(r=argvals[domain],
obs=fX[domain],
theo=reference,
lo=lo,
hi=hi)
} else {
results <- data.frame(r=argvals[domain],
obs=fX[domain],
mmean=reference,
lo=lo,
hi=hi)
}
colnames(results)[1] <- argname
shadenames <- c("lo", "hi")
if(do.pwrong)
warning(paste("Argument", sQuote("do.pwrong=TRUE"), "ignored;",
"it is not relevant to global envelopes"))
},
variance={
## ....... POINTWISE MEAN, VARIANCE etc ......................
if(gaveup) {
Ef <- varf <- NAvector
} else if(cheat) {
# .... use precomputed values ....
stopifnot(checkfields(precomputed, c("Ef", "varf")))
Ef <- precomputed$Ef
varf <- precomputed$varf
} else {
## .... normal case: compute from simulations
simvals[is.infinite(simvals)] <- NA
if(restrict.columns) {
simvals <- simvals[, jsim]
if(use.weights) weights <- weights[jsim]
}
nsim <- ncol(simvals)
if(!use.weights) {
Ef <- apply(simvals, 1L, mean, na.rm=TRUE)
varf <- apply(simvals, 1L, var, na.rm=TRUE)
} else {
Ef <- apply(simvals, 1L, weighted.mean, w=weights, na.rm=TRUE)
varf <- apply(simvals, 1L, weighted.var, w=weights, na.rm=TRUE)
}
}
if(gaveup) {
sd <- stdres <- lo <- hi <- loCI <- hiCI <- NAvector
} else {
## derived quantities
sd <- sqrt(varf)
stdres <- (fX-Ef)/sd
stdres[!is.finite(stdres)] <- NA
## critical limits
lo <- Ef - nSD * sd
hi <- Ef + nSD * sd
## confidence interval
loCI <- Ef - nSD * sd/sqrt(nsim)
hiCI <- Ef + nSD * sd/sqrt(nsim)
}
lo.name <- paste("lower", nSD, "sigma critical limit for %s")
hi.name <- paste("upper", nSD, "sigma critical limit for %s")
loCI.name <- paste("lower", nSD, "sigma confidence bound",
"for mean of simulated %s")
hiCI.name <- paste("upper", nSD, "sigma confidence bound",
"for mean of simulated %s")
##
if(!gaveup)
switch(alternative,
two.sided = { },
less = {
hi <- hiCI <- rep.int(Inf, length(hi))
hi.name <- "infinite upper boundary"
hiCI.name <- "infinite upper confidence limit"
},
greater = {
lo <- loCI <- rep.int(-Inf, length(lo))
lo.name <- "infinite lower boundary"
loCI.name <- "infinite lower confidence limit"
})
## put together
if(use.theory) {
results <- data.frame(r=argvals,
obs=fX,
theo=theory,
lo=lo,
hi=hi)
colnames(results)[1] <- argname
shadenames <- c("lo", "hi")
morestuff <- data.frame(mmean=Ef,
var=varf,
res=fX-Ef,
stdres=stdres,
loCI=loCI,
hiCI=hiCI)
loCIlabel <-
if(alternative == "greater" && !gaveup) "-infinity" else
makefvlabel(NULL, NULL, fname, "loCI", argname=argname)
hiCIlabel <-
if(alternative == "less" && !gaveup) "infinity" else
makefvlabel(NULL, NULL, fname, "hiCI", argname=argname)
mslabl <- c(makefvlabel(NULL, "bar", fname, argname=argname),
makefvlabel("var", "hat", fname, argname=argname),
makefvlabel("res", "hat", fname, argname=argname),
makefvlabel("stdres", "hat", fname, argname=argname),
loCIlabel,
hiCIlabel)
wted <- if(use.weights) "weighted " else NULL
msdesc <- c(paste0(wted, "sample mean of %s from simulations"),
paste0(wted, "sample variance of %s from simulations"),
"raw residual",
"standardised residual",
loCI.name, hiCI.name)
} else {
results <- data.frame(r=argvals,
obs=fX,
mmean=Ef,
lo=lo,
hi=hi)
colnames(results)[1] <- argname
shadenames <- c("lo", "hi")
morestuff <- data.frame(var=varf,
res=fX-Ef,
stdres=stdres,
loCI=loCI,
hiCI=hiCI)
loCIlabel <-
if(alternative == "greater" && !gaveup) "-infinity" else
makefvlabel(NULL, NULL, fname, "loCI", argname=argname)
hiCIlabel <-
if(alternative == "less" && !gaveup) "infinity" else
makefvlabel(NULL, NULL, fname, "hiCI", argname=argname)
mslabl <- c(makefvlabel("var", "hat", fname, argname=argname),
makefvlabel("res", "hat", fname, argname=argname),
makefvlabel("stdres", "hat", fname, argname=argname),
loCIlabel,
hiCIlabel)
msdesc <- c(paste0(if(use.weights) "weighted " else NULL,
"sample variance of %s from simulations"),
"raw residual",
"standardised residual",
loCI.name, hiCI.name)
}
if(do.pwrong) {
## estimate the p-value for the 'wrong test'
if(gaveup) {
pwrong <- NA_real_
} else if(cheat) {
pwrong <- precomputed$pwrong
do.pwrong <- !is.null(pwrong) && !badprobability(pwrong, FALSE)
} else {
upper.signif <- (simvals > hi)
lower.signif <- (simvals < lo)
is.signif <- switch(alternative,
less = lower.signif,
greater = upper.signif,
two.sided = lower.signif | upper.signif)
# is.signif.somewhere <- apply(is.signif, 2, any)
is.signif.somewhere <- matcolany(is.signif)
pwrong <- sum(is.signif.somewhere)/nsim
}
}
}
)
############ WRAP UP #########################
if(use.theory) {
# reference is computed curve `theo'
reflabl <- makefvlabel(NULL, NULL, fname, "theo", argname=argname)
refdesc <- paste0("theoretical value of %s", if(csr) " for CSR" else NULL)
} else {
# reference is sample mean of simulations
reflabl <- makefvlabel(NULL, "bar", fname, argname=argname)
refdesc <- paste0(if(use.weights) "weighted " else NULL,
"sample mean of %s from simulations")
}
lolabl <- if(alternative == "greater" && !gaveup) "-infinity" else
makefvlabel(NULL, "hat", fname, "lo", argname=argname)
hilabl <- if(alternative == "less"&& !gaveup) "infinity" else
makefvlabel(NULL, "hat", fname, "hi", argname=argname)
result <- fv(results,
argu=argname,
ylab=atr$ylab,
valu="obs",
fmla= paste(". ~", argname),
alim=intersect.ranges(atr$alim, range(results[[argname]])),
labl=c(argname,
makefvlabel(NULL, "hat", fname, "obs", argname=argname),
reflabl,
lolabl,
hilabl),
desc=c(arg.desc,
"observed value of %s for data pattern",
refdesc, lo.name, hi.name),
fname=atr$fname,
yexp =atr$yexp)
# columns to be plotted by default
dotty <- c("obs", if(use.theory) "theo" else "mmean", "hi", "lo")
if(type == "variance") {
# add more stuff
result <- bind.fv(result, morestuff, mslabl, msdesc)
if(use.theory) dotty <- c(dotty, "mmean")
}
fvnames(result, ".") <- dotty
fvnames(result, ".s") <- shadenames
unitname(result) <- unitname(funX)
class(result) <- c("envelope", class(result))
# tack on envelope information
attr(result, "einfo") <- list(global = (type =="global"),
ginterval = ginterval,
alternative=alternative,
scale = scale,
clamp = clamp,
csr = csr,
use.theory = use.theory,
csr.theo = csr && use.theory,
simtype = "funs",
constraints = "",
nrank = nrank,
nsim = nsim,
VARIANCE = (type == "variance"),
nSD = nSD,
valname = NULL,
dual = dual,
nsim = nsim,
nsim2 = nsim.mean,
Yname = Yname,
do.pwrong=do.pwrong,
use.weights=use.weights,
gaveup = gaveup)
# tack on saved functions
if(savefuns && !gaveup) {
nSim <- ncol(Y)
alldata <- cbind(argvals, Y)
simnames <- paste("sim", 1:nSim, sep="")
colnames(alldata) <- c(argname, simnames)
alldata <- as.data.frame(alldata)
SimFuns <- fv(alldata,
argu=argname,
ylab=atr$ylab,
valu="sim1",
fmla= paste(". ~", argname),
alim=atr$alim,
labl=names(alldata),
desc=c(arg.desc,
paste("Simulation ", 1:nSim, sep="")),
unitname=unitname(funX))
fvnames(SimFuns, ".") <- simnames
attr(result, "simfuns") <- SimFuns
}
if(do.pwrong)
attr(result, "pwrong") <- pwrong
if(use.weights)
attr(result, "weights") <- weights
return(result)
}
envelope.envelope <- function(Y, fun=NULL, ...,
transform=NULL, global=FALSE, VARIANCE=FALSE) {
Yname <- short.deparse(substitute(Y))
stopifnot(inherits(Y, "envelope"))
Yorig <- Y
aargh <- list(...)
X <- attr(Y, "datapattern")
sf <- attr(Y, "simfuns")
sp <- attr(Y, "simpatterns")
wt <- attr(Y, "weights")
einfo <- attr(Y, "einfo")
csr <- aargh$internal$csr %orifnull% einfo$csr
if(is.null(fun) && is.null(sf)) {
# No simulated functions - must compute them from simulated patterns
if(is.null(sp))
stop(paste("Cannot compute envelope:",
"Y does not contain simulated functions",
"(was not generated with savefuns=TRUE)",
"and does not contain simulated patterns",
"(was not generated with savepatterns=TRUE)"))
# set default fun=Kest
fun <- Kest
}
if(!is.null(fun)) {
# apply new function
# point patterns are required
if(is.null(sp))
stop(paste("Object Y does not contain simulated point patterns",
"(attribute", dQuote("simpatterns"), ");",
"cannot apply a new", sQuote("fun")))
if(is.null(X))
stop(paste("Cannot apply a new", sQuote("fun"),
"; object Y generated by an older version of spatstat"))
## send signal if simulations were CSR
internal <- aargh$internal
if(csr) {
if(is.null(internal)) internal <- list()
internal$csr <- TRUE
}
## compute new envelope
result <- do.call(envelope,
resolve.defaults(list(Y=quote(X), fun=fun, simulate=sp),
aargh,
list(transform=transform,
global=global,
VARIANCE=VARIANCE,
internal=internal,
Yname=Yname,
nsim=einfo$nsim,
nsim2=einfo$nsim2,
weights=wt),
.StripNull=TRUE))
} else {
#' compute new envelope with existing simulated functions
if(is.null(sf))
stop(paste("Y does not contain a", dQuote("simfuns"), "attribute",
"(it was not generated with savefuns=TRUE)"))
if(!is.null(transform)) {
# Apply transformation to Y and sf
stopifnot(is.expression(transform))
## cc <- dotexpr.to.call(transform, "Y", "eval.fv")
cc <- inject.expr("with(Y, .)", transform)
Y <- eval(cc)
## cc <- dotexpr.to.call(transform, "sf", "eval.fv")
cc <- inject.expr("with(sf, .)", transform)
sf <- eval(cc)
}
#' catch discrepancy between domains of observed and simulated functions
if(nrow(sf) != nrow(Y)) {
rrsim <- sf[[fvnames(sf, ".x")]]
rrobs <- Y[[fvnames(Y, ".x")]]
ra <- intersect.ranges(range(rrsim), range(rrobs))
delta <- min(mean(diff(rrsim)), mean(diff(rrobs)))/2
oksim <- (rrsim >= ra[1] - delta) & (rrsim <= ra[2] + delta)
okobs <- (rrobs >= ra[1] - delta) & (rrobs <= ra[2] + delta)
if(sum(oksim) != sum(okobs))
stop("Internal error: Unable to reconcile the domains",
"of the observed and simulated functions", call.=FALSE)
if(mean(abs(rrsim[oksim] - rrobs[okobs])) > delta)
stop("Internal error: Unable to reconcile the r values",
"of the observed and simulated functions", call.=FALSE)
sf <- sf[oksim, ,drop=FALSE]
Y <- Y[okobs, ,drop=FALSE]
}
# extract simulated function values
df <- as.data.frame(sf)
rname <- fvnames(sf, ".x")
df <- df[, (names(df) != rname)]
# interface with 'envelope.matrix'
etype <- if(global) "global" else if(VARIANCE) "variance" else "pointwise"
dfm <- as.matrix(df)
dont.complain.about(dfm)
result <- do.call(envelope.matrix,
resolve.defaults(list(Y=quote(dfm)),
aargh,
list(type=etype,
csr=csr,
funX=Y,
Yname=Yname,
weights=wt),
.StripNull=TRUE))
}
if(!is.null(transform)) {
# post-process labels
labl <- attr(result, "labl")
dnames <- colnames(result)
dnames <- dnames[dnames %in% fvnames(result, ".")]
# expand "."
ud <- as.call(lapply(c("cbind", dnames), as.name))
dont.complain.about(ud)
expandtransform <- eval(substitute(substitute(tr, list(.=ud)),
list(tr=transform[[1L]])))
# compute new labels
attr(result, "fname") <- attr(Yorig, "fname")
mathlabl <- as.character(fvlegend(result, expandtransform))
# match labels to columns
evars <- all.vars(expandtransform)
used.dotnames <- evars[evars %in% dnames]
mathmap <- match(colnames(result), used.dotnames)
okmath <- !is.na(mathmap)
# update appropriate labels
labl[okmath] <- mathlabl[mathmap[okmath]]
attr(result, "labl") <- labl
}
# Tack on envelope info
copyacross <- c("Yname", "csr.theo", "use.theory", "simtype", "constraints")
attr(result, "einfo")[copyacross] <- attr(Yorig, "einfo")[copyacross]
attr(result, "einfo")$csr <- csr
# Save data
return(result)
}
pool.envelope <- local({
pool.envelope <- function(..., savefuns=FALSE, savepatterns=FALSE) {
Yname <- short.deparse(sys.call())
if(nchar(Yname) > 60) Yname <- paste(substr(Yname, 1L, 40L), "[..]")
Elist <- unname(list(...))
nE <- length(Elist)
if(nE == 0) return(NULL)
#' ........ validate envelopes .....................
#' All arguments must be envelopes
notenv <- !unlist(lapply(Elist, inherits, what="envelope"))
if(any(notenv)) {
n <- sum(notenv)
why <- paste(ngettext(n, "Argument", "Arguments"),
commasep(which(notenv)),
ngettext(n, "does not", "do not"),
"belong to the class",
dQuote("envelope"))
stop(why)
}
E1 <- Elist[[1L]]
## Only one envelope?
if(nE == 1)
return(E1)
## envelopes must be compatible
ok <- do.call(compatible, Elist)
if(!ok)
stop("Envelopes are not compatible")
## find name of function argument
argname <- fvnames(E1, ".x")
arg.desc <- attr(E1, "desc")[match(argname, colnames(E1))]
## ... reconcile parameters in different envelopes .......
eilist <- lapply(Elist, attr, which="einfo")
nrank <- resolveEinfo(eilist, "nrank", 1)
global <- resolveEinfo(eilist, "global", FALSE)
ginterval <- resolveEinfo(eilist, "ginterval", NULL, atomic=FALSE)
VARIANCE <- resolveEinfo(eilist, "VARIANCE", FALSE)
alternative <- resolveEinfo(eilist, "alternative", FALSE)
scale <- resolveEinfo(eilist, "scale", NULL, atomic=FALSE)
clamp <- resolveEinfo(eilist, "clamp", FALSE)
resolveEinfo(eilist, "simtype", "funs",
"Envelopes were generated using different types of simulation")
resolveEinfo(eilist, "constraints", "",
"Envelopes were generated using different types of conditioning")
resolveEinfo(eilist, "csr.theo", FALSE, NULL)
csr <- resolveEinfo(eilist, "csr", FALSE, NULL)
use.weights <- resolveEinfo(eilist, "use.weights" , FALSE,
"Weights were used in some, but not all, envelopes: they will be ignored")
use.theory <- resolveEinfo(eilist, "use.theory", csr, NULL)
##
weights <-
if(use.weights) unlist(lapply(Elist, attr, which="weights")) else NULL
type <- if(global) "global" else if(VARIANCE) "variance" else "pointwise"
## ........ validate saved functions .....................
if(savefuns || !VARIANCE) {
## Individual simulated functions are required
SFlist <- lapply(Elist, attr, which="simfuns")
isnul <- unlist(lapply(SFlist, is.null))
if(any(isnul)) {
n <- sum(isnul)
comply <- if(!VARIANCE) "compute the envelope:" else
"save the simulated functions:"
why <- paste("Cannot", comply,
ngettext(n, "argument", "arguments"),
commasep(which(isnul)),
ngettext(n, "does not", "do not"),
"contain a", dQuote("simfuns"), "attribute",
"(not generated with savefuns=TRUE)")
stop(why)
}
## Simulated functions must be the same function
fnames <- unique(lapply(SFlist, attr, which="fname"))
if(length(fnames) > 1L) {
fnames <- unlist(lapply(fnames, flatfname))
stop(paste("Envelope objects contain values",
"of different functions:",
commasep(sQuote(fnames))))
}
## vectors of r values must be identical
rlist <- lapply(SFlist, getargvals)
argvals <- rlist[[1L]]
samer <- unlist(lapply(rlist, identical, y=argvals))
if(!all(samer))
stop(paste("Simulated function values are not compatible",
"(different values of function argument)"))
## Extract function values and assemble into one matrix
matlist <- lapply(SFlist, getdotvals)
SFmatrix <- do.call(cbind, matlist)
}
## compute pooled envelope
switch(type,
pointwise = {
result <- envelope(SFmatrix, funX=E1,
type=type, alternative=alternative,
clamp=clamp,
nrank=nrank,
csr=csr, use.theory=use.theory,
Yname=Yname, weights=weights,
savefuns=savefuns)
},
global = {
simfunmatrix <- if(is.null(ginterval)) SFmatrix else {
## savefuns have not yet been clipped to ginterval
## while envelope data have been clipped.
domain <- (argvals >= ginterval[1L]) & (argvals <= ginterval[2L])
SFmatrix[domain, , drop=FALSE]
}
result <- envelope(simfunmatrix, funX=E1,
type=type, alternative=alternative,
scale=scale, clamp=clamp,
csr=csr, use.theory=use.theory,
nrank=nrank,
ginterval=ginterval,
Yname=Yname, weights=weights,
savefuns=savefuns)
},
variance = {
## Pool sample means and variances
nsims <- unlist(lapply(eilist, getElement, name="nsim"))
mmeans <- lapply(Elist, getElement, name="mmean")
vars <- lapply(Elist, getElement, name="var")
mmeans <- matrix(unlist(mmeans), ncol=nE)
vars <- matrix(unlist(vars), ncol=nE)
if(!use.weights) {
w.mean <- nsims
d.mean <- sum(nsims)
w.var <- nsims - 1
d.var <- sum(nsims) - 1
} else {
weightlist <- lapply(Elist, attr, which="weights")
w.mean <- unlist(lapply(weightlist, sum))
d.mean <- sum(w.mean)
ssw <- unlist(lapply(weightlist, meansqfrac))
## meansqfrac : function(x) {sum((x/sum(x))^2)}))
w.var <- w.mean * (1 - ssw)
d.var <- d.mean * (1 - sum(ssw))
}
poolmmean <- as.numeric(mmeans %*% matrix(w.mean/d.mean, ncol=1L))
within <- vars %*% matrix(w.var, ncol=1L)
between <- ((mmeans - poolmmean[])^2) %*% matrix(w.mean, ncol=1L)
poolvar <- as.numeric((within + between)/d.var)
## feed precomputed data to envelope.matrix
pc <- list(Ef=poolmmean[],
varf=poolvar[])
nsim <- sum(nsims)
result <- envelope.matrix(NULL, funX=E1,
type=type, alternative=alternative,
csr=csr, Yname=Yname,
weights=weights,
savefuns=savefuns,
nsim=nsim,
precomputed=pc)
})
## Copy envelope info that is not handled by envelope.matrix
copyacross <- c("Yname", "csr.theo", "use.theory", "simtype", "constraints")
attr(result, "einfo")[copyacross] <- attr(E1, "einfo")[copyacross]
## ..............saved patterns .....................
if(savepatterns) {
SPlist <- lapply(Elist, attr, which="simpatterns")
isnul <- unlist(lapply(SPlist, is.null))
if(any(isnul)) {
n <- sum(isnul)
why <- paste("Cannot save the simulated patterns:",
ngettext(n, "argument", "arguments"),
commasep(which(isnul)),
ngettext(n, "does not", "do not"),
"contain a", dQuote("simpatterns"), "attribute",
"(not generated with savepatterns=TRUE)")
warning(why)
} else {
attr(result, "simpatterns") <- Reduce(append, SPlist)
}
}
## ..............saved summary functions ................
if(savefuns) {
alldata <- cbind(argvals, SFmatrix)
Nsim <- ncol(SFmatrix)
simnames <- paste0("sim", 1:Nsim)
colnames(alldata) <- c(argname, simnames)
alldata <- as.data.frame(alldata)
SFtemplate <- SFlist[[1L]]
SimFuns <- fv(alldata,
argu=argname,
ylab=attr(SFtemplate, "ylab"),
valu="sim1",
fmla=paste(". ~", argname),
alim=attr(SFtemplate, "alim"),
labl=names(alldata),
desc=c(arg.desc,
paste("Simulation ", 1:Nsim, sep="")),
fname=attr(SFtemplate, "fname"),
yexp=attr(SFtemplate, "yexp"),
unitname=unitname(SFtemplate))
fvnames(SimFuns, ".") <- simnames
attr(result, "simfuns") <- SimFuns
}
dotnames <- lapply(Elist, fvnames, a=".")
dn <- dotnames[[1L]]
if(all(unlist(lapply(dotnames, identical, y=dn))))
fvnames(result, ".") <- dn
shadenames <- lapply(Elist, fvnames, a=".s")
sh <- shadenames[[1L]]
if(all(unlist(lapply(shadenames, identical, y=sh))))
fvnames(result, ".s") <- sh
return(result)
}
getargvals <- function(z) { as.matrix(z)[, fvnames(z, ".x")] }
getdotvals <- function(z) { as.matrix(z)[, fvnames(z, "."), drop=FALSE] }
meansqfrac <- function(x) {sum((x/sum(x))^2)}
pool.envelope
})
# resolve matching entries in different envelope objects
# x is a list of envelope info objects
resolveEinfo <- function(x, what, fallback, warn, atomic=TRUE) {
if(atomic) {
y <- unique(unlist(lapply(x, getElement, name=what)))
if(length(y) == 1L)
return(y)
} else {
y <- unique(lapply(x, getElement, name=what))
if(length(y) == 1L)
return(y[[1L]])
}
if(missing(warn))
warn <- paste("Envelopes were generated using different values",
"of argument", paste(sQuote(what), ";", sep=""),
"reverting to default value")
if(!is.null(warn))
warning(warn, call.=FALSE)
return(fallback)
}
as.data.frame.envelope <- function(x, ..., simfuns=FALSE) {
if(simfuns && !is.null(sf <- attr(x, "simfuns"))) {
# tack on the simulated functions as well
y <- as.data.frame(bind.fv(x, sf, clip=TRUE))
return(y)
}
NextMethod("as.data.frame")
}
Try the spatstat.explore package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.