https://github.com/cran/spatstat
Revision 967fc2de73b44e75d7d8497a8d21770604051e17 authored by Adrian Baddeley on 21 December 2011, 08:53:46 UTC, committed by cran-robot on 21 December 2011, 08:53:46 UTC
1 parent 266eff5
Tip revision: 967fc2de73b44e75d7d8497a8d21770604051e17 authored by Adrian Baddeley on 21 December 2011, 08:53:46 UTC
version 1.25-1
version 1.25-1
Tip revision: 967fc2d
envelope.R
#
# envelope.R
#
# computes simulation envelopes
#
# $Revision: 2.22 $ $Date: 2011/10/07 08:52:09 $
#
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
# "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
#
# ...................................................................
simulrecipe <- function(type, expr, envir, csr, pois=csr) {
if(csr && !pois) warning("Internal error: csr=TRUE but pois=FALSE")
out <- list(type=type,
expr=expr,
envir=envir,
csr=csr,
pois=pois)
class(out) <- "simulrecipe"
out
}
envelope.ppp <-
function(Y, fun=Kest, nsim=99, nrank=1, ...,
simulate=NULL, verbose=TRUE, clipdata=TRUE,
transform=NULL, global=FALSE, ginterval=NULL,
savefuns=FALSE, savepatterns=FALSE, nsim2=nsim,
VARIANCE=FALSE, nSD=2,
Yname=NULL, maxnerr=nsim) {
cl <- match.call()
if(is.null(Yname)) Yname <- deparse(substitute(Y))
envir.user <- parent.frame()
envir.here <- sys.frame(sys.nframe())
if(is.null(simulate)) {
# ...................................................
# Realisations of complete spatial randomness
# will be generated by rpoispp
# Data pattern X is argument Y
# Data pattern determines intensity of Poisson process
X <- Y
sY <- summary(Y, checkdup=FALSE)
Yintens <- sY$intensity
Ywin <- Y$window
# expression that will be evaluated
simexpr <-
if(!is.marked(Y)) {
# unmarked point pattern
expression(rpoispp(Yintens, win=Ywin))
} else {
# marked point pattern
Ymarx <- marks(Y, dfok=FALSE)
expression({A <- rpoispp(Yintens, win=Ywin);
A %mark% sample(Ymarx, A$n, replace=TRUE)})
}
# evaluate in THIS environment
simrecipe <- simulrecipe(type = "csr",
expr = simexpr,
envir = envir.here,
csr = TRUE,
pois = TRUE)
} else {
# ...................................................
# Simulations are determined by 'simulate' argument
# Processing is deferred to envelopeEngine
simrecipe <- simulate
# Data pattern is argument Y
X <- Y
}
envelopeEngine(X=X, fun=fun, simul=simrecipe,
nsim=nsim, nrank=nrank, ...,
verbose=verbose, clipdata=clipdata,
transform=transform, global=global, ginterval=ginterval,
savefuns=savefuns, savepatterns=savepatterns, nsim2=nsim2,
VARIANCE=VARIANCE, nSD=nSD,
Yname=Yname, maxnerr=maxnerr, cl=cl,
envir.user=envir.user)
}
envelope.ppm <-
function(Y, fun=Kest, nsim=99, nrank=1, ...,
simulate=NULL, verbose=TRUE, clipdata=TRUE,
start=NULL, control=default.rmhcontrol(Y, nrep=nrep), nrep=1e5,
transform=NULL, global=FALSE, ginterval=NULL,
savefuns=FALSE, savepatterns=FALSE, nsim2=nsim,
VARIANCE=FALSE, nSD=2,
Yname=NULL, maxnerr=nsim) {
cl <- match.call()
if(is.null(Yname)) Yname <- deparse(substitute(Y))
envir.user <- parent.frame()
envir.here <- sys.frame(sys.nframe())
# Extract data pattern X from fitted model Y
X <- data.ppm(Y)
if(is.null(simulate)) {
# ...................................................
# Simulated realisations of the fitted model Y
# will be generated
pois <- is.poisson(Y)
csr <- is.stationary(Y) && pois
type <- if(csr) "csr" else "rmh"
# Set up parameters for rmh
rmodel <- rmhmodel(Y, verbose=FALSE)
if(is.null(start))
start <- list(n.start=X$n)
rstart <- rmhstart(start)
rcontr <- rmhcontrol(control)
# pre-digest arguments
rmhinfolist <- rmh(rmodel, rstart, rcontr, preponly=TRUE, verbose=FALSE)
# expression that will be evaluated
simexpr <- expression(rmhEngine(rmhinfolist, verbose=FALSE))
envir <- envir.here
# evaluate in THIS environment
simrecipe <- simulrecipe(type = type,
expr = simexpr,
envir = envir.here,
csr = csr,
pois = pois)
} else {
# ...................................................
# Simulations are determined by 'simulate' argument
# Processing is deferred to envelopeEngine
simrecipe <- simulate
}
envelopeEngine(X=X, fun=fun, simul=simrecipe,
nsim=nsim, nrank=nrank, ...,
verbose=verbose, clipdata=clipdata,
transform=transform, global=global, ginterval=ginterval,
savefuns=savefuns, savepatterns=savepatterns, nsim2=nsim2,
VARIANCE=VARIANCE, nSD=nSD,
Yname=Yname, maxnerr=maxnerr, cl=cl,
envir.user=envir.user)
}
envelope.kppm <-
function(Y, fun=Kest, nsim=99, nrank=1, ...,
simulate=NULL, verbose=TRUE, clipdata=TRUE,
transform=NULL, global=FALSE, ginterval=NULL,
savefuns=FALSE, savepatterns=FALSE, nsim2=nsim,
VARIANCE=FALSE, nSD=2, Yname=NULL, maxnerr=nsim)
{
cl <- match.call()
if(is.null(Yname)) Yname <- deparse(substitute(Y))
envir.user <- parent.frame()
envir.here <- sys.frame(sys.nframe())
# Extract data pattern X from fitted model Y
X <- Y$X
if(is.null(simulate)) {
# Simulated realisations of the fitted model Y
# will be generated using simulate.kppm
kmodel <- Y
# expression that will be evaluated
simexpr <- expression(simulate(kmodel)[[1]])
# evaluate in THIS environment
simrecipe <- simulrecipe(type = "kppm",
expr = simexpr,
envir = envir.here,
csr = FALSE,
pois = FALSE)
} else {
# ...................................................
# Simulations are determined by 'simulate' argument
# Processing is deferred to envelopeEngine
simrecipe <- simulate
}
envelopeEngine(X=X, fun=fun, simul=simrecipe,
nsim=nsim, nrank=nrank, ...,
verbose=verbose, clipdata=clipdata,
transform=transform, global=global, ginterval=ginterval,
savefuns=savefuns, savepatterns=savepatterns, nsim2=nsim2,
VARIANCE=VARIANCE, nSD=nSD,
Yname=Yname, maxnerr=maxnerr, cl=cl,
envir.user=envir.user)
}
## .................................................................
## Engine for simulating and computing envelopes
## .................................................................
#
# X is the data point pattern, which could be ppp, pp3, ppx etc
# X determines the class of pattern expected from the simulations
#
envelopeEngine <-
function(X, fun, simul,
nsim=99, nrank=1, ...,
verbose=TRUE, clipdata=TRUE,
transform=NULL, global=FALSE, ginterval=NULL,
savefuns=FALSE, savepatterns=FALSE, nsim2=nsim,
VARIANCE=FALSE, nSD=2,
Yname=NULL, maxnerr=nsim, internal=NULL, cl=NULL,
envir.user=envir.user) {
#
envir.here <- sys.frame(sys.nframe())
# ----------------------------------------------------------
# Determine Simulation
# ----------------------------------------------------------
# Identify class of patterns to be simulated, from data pattern X
Xclass <- if(is.ppp(X)) "ppp" else
if(is.pp3(X)) "pp3" else
if(is.ppx(X)) "ppx" else
stop("Unrecognised class of point pattern")
Xobjectname <- paste("point pattern of class", sQuote(Xclass))
# Identify type of simulation from argument 'simul'
if(inherits(simul, "simulrecipe")) {
# ..................................................
# simulation recipe is given
simtype <- simul$type
simexpr <- simul$expr
envir <- simul$envir
csr <- simul$csr
pois <- simul$pois
} else {
# ...................................................
# simulation is specified by argument `simulate' to envelope()
simulate <- simul
# which should be an expression, or a list of point patterns,
# or an envelope object.
csr <- FALSE
# override
if(!is.null(icsr <- internal$csr)) csr <- icsr
pois <- csr
model <- NULL
if(inherits(simulate, "envelope")) {
# envelope object: see if it contains stored point patterns
simpat <- attr(simulate, "simpatterns")
if(!is.null(simpat))
simulate <- simpat
else
stop(paste("The argument", sQuote("simulate"),
"is an envelope object but does not contain",
"any saved point patterns."))
}
if(is.expression(simulate)) {
# The user-supplied expression 'simulate' will be evaluated repeatedly
simtype <- "expr"
simexpr <- simulate
envir <- envir.user
} else if(is.list(simulate) &&
( (is.ppp(X) && all(unlist(lapply(simulate, is.ppp))))
|| (is.pp3(X) && all(unlist(lapply(simulate, is.pp3))))
|| (is.ppx(X) && all(unlist(lapply(simulate, is.ppx)))))) {
# The user-supplied list of point patterns will be used
simtype <- "list"
SimDataList <- simulate
# expression that will be evaluated
simexpr <- expression(SimDataList[[i]])
envir <- envir.here
# ensure that `i' is defined
i <- 1
# any messages?
if(!is.null(mess <- attr(simulate, "internal"))) {
# determine whether these point patterns are realisations of CSR
csr <- !is.null(mc <- mess$csr) && mc
}
} else stop(paste(sQuote("simulate"),
"should be an expression, or a list of point patterns"))
}
# -------------------------------------------------------------------
# Determine clipping window
# ------------------------------------------------------------------
if(clipdata) {
# Generate one realisation
Xsim <- eval(simexpr, envir=envir)
if(!inherits(Xsim, Xclass))
switch(simtype,
csr=stop(paste("Internal error:", Xobjectname, "not generated")),
rmh=stop(paste("Internal error: rmh did not return an",
Xobjectname)),
kppm=stop(paste("Internal error: simulate.kppm did not return an",
Xobjectname)),
expr=stop(paste("Evaluating the expression", sQuote("simulate"),
"did not yield an", Xobjectname)),
list=stop(paste("Internal error: list entry was not an",
Xobjectname)),
stop(paste("Internal error:", Xobjectname, "not generated"))
)
# Extract window
clipwin <- Xsim$window
if(!is.subset.owin(clipwin, X$window))
warning("Window containing simulated patterns is not a subset of data window")
}
# ------------------------------------------------------------------
# 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))) deparse(substitute(fun))
else if(is.character(fun)) fun else "fun"
fname <- sQuote(fname)
# R function to apply
if(is.character(fun))
fun <- get(fun)
if(!is.function(fun))
stop(paste("unrecognised format for function", fname))
fargs <- names(formals(fun))
if(!any(c("r", "...") %in% fargs))
stop(paste(fname, "should have an argument",
sQuote("r"), "or", sQuote("...")))
usecorrection <- any(c("correction", "...") %in% fargs)
# ---------------------------------------------------------------------
# validate other arguments
if((nrank %% 1) != 0)
stop("nrank must be an integer")
stopifnot(nrank > 0 && nrank < nsim/2)
rgiven <- "r" %in% names(list(...))
if(tran <- !is.null(transform)) {
stopifnot(is.expression(transform))
# 'transform' is an expression
aa <- substitute(substitute(tt, list(.=as.name("funX"))),
list(tt=transform))
# 'aa' is a language expression invoking 'substitute'
bb <- eval(parse(text=deparse(aa)))
# 'bb' is an expression obtained by replacing "." by "funX"
transform.funX <- as.call(bb)
transform.funX[[1]] <- as.name("eval.fv")
# 'transform.funX' is an unevaluated call to eval.fv
aa <- substitute(substitute(tt, list(.=as.name("funXsim"))),
list(tt=transform))
bb <- eval(parse(text=deparse(aa)))
transform.funXsim <- as.call(bb)
transform.funXsim[[1]] <- as.name("eval.fv")
}
if(!is.null(ginterval))
stopifnot(is.numeric(ginterval) && length(ginterval) == 2)
# ---------------------------------------------------------------------
# Evaluate function for data pattern X
# ------------------------------------------------------------------
Xarg <- if(!clipdata) X else X[clipwin]
funX <- do.call(fun,
resolve.defaults(list(Xarg),
list(...),
if(usecorrection) list(correction="best") else NULL))
if(!inherits(funX, "fv"))
stop(paste("The function", fname,
"must return an object of class", sQuote("fv")))
argname <- fvnames(funX, ".x")
valname <- fvnames(funX, ".y")
has.theo <- "theo" %in% fvnames(funX, "*")
csr.theo <- csr && has.theo
if(tran) {
# extract only the recommended value
if(csr.theo)
funX <- funX[, c(argname, valname, "theo")]
else
funX <- funX[, c(argname, valname)]
# apply the transformation to it
funX <- eval(transform.funX)
}
rvals <- funX[[argname]]
fX <- funX[[valname]]
# default domain over which to maximise
alim <- attr(funX, "alim")
if(global && is.null(ginterval))
ginterval <- if(rgiven) range(rvals) else alim
#--------------------------------------------------------------------
# Determine number of simulations
# ------------------------------------------------------------------
#
## determine whether dual simulations are required
## (one set of simulations to calculate the theoretical mean,
## another independent set of simulations to obtain the critical point.)
dual <- (global && !csr.theo && !VARIANCE)
Nsim <- if(!dual) nsim else (nsim + nsim2)
# if taking data from a list of point patterns,
# check there are enough of them
if(simtype == "list" && Nsim > length(SimDataList))
stop(paste("Number of simulations",
paren(if(!dual)
paste(nsim) else
paste(nsim, "+", nsim2, "=", Nsim)
),
"exceeds number of point pattern datasets supplied"))
# Undocumented secret exit
# ------------------------------------------------------------------
if(!is.null(eject <- internal$eject) && eject == "patterns") {
# generate simulated realisations and return only these patterns
if(verbose) {
action <- if(simtype == "list") "Extracting" else "Generating"
descrip <- switch(simtype,
csr = "simulations of CSR",
rmh = paste("simulated realisations of fitted",
if(pois) "Poisson" else "Gibbs",
"model"),
kppm = "simulated realisations of fitted cluster model",
expr = "simulations by evaluating expression",
list = "point patterns from list",
"simulated realisations")
explan <- if(dual) paren(paste(nsim2, "to estimate the mean and",
nsim, "to calculate envelopes")) else ""
cat(paste(action, Nsim, descrip, explan, "...\n"))
}
XsimList <- list()
# start simulation loop
for(i in 1:Nsim) {
if(verbose) progressreport(i, Nsim)
Xsim <- eval(simexpr, envir=envir)
if(!inherits(Xsim, Xclass))
switch(simtype,
csr={
stop(paste("Internal error:", Xobjectname, "not generated"))
},
rmh={
stop(paste("Internal error: rmh did not return an",
Xobjectname))
},
kppm={
stop(paste("Internal error: simulate.kppm did not return an",
Xobjectname))
},
expr={
stop(paste("Evaluating the expression", sQuote("simulate"),
"did not yield an", Xobjectname))
},
list={
stop(paste("Internal error: list entry was not an",
Xobjectname))
},
stop(paste("Internal error:", Xobjectname, "not generated"))
)
XsimList[[i]] <- Xsim
}
if(verbose) {
cat(paste("Done.\n"))
flush.console()
}
attr(XsimList, "internal") <- list(csr=csr)
return(XsimList)
}
# capture main decision parameters
EnvelopeInfo <- list(call=cl,
Yname=Yname,
valname=valname,
csr=csr,
csr.theo=csr.theo,
pois=pois,
simtype=simtype,
nrank=nrank,
nsim=nsim,
Nsim=Nsim,
global=global,
dual=dual,
nsim2=nsim2,
VARIANCE=VARIANCE,
nSD=nSD)
# ----------------------------------------
######### SIMULATE #######################
# ----------------------------------------
if(verbose) {
action <- if(simtype == "list") "Extracting" else "Generating"
descrip <- switch(simtype,
csr = "simulations of CSR",
rmh = paste("simulated realisations of fitted",
if(pois) "Poisson" else "Gibbs",
"model"),
kppm = "simulated realisations of fitted cluster model",
expr = "simulations by evaluating expression",
list = "point patterns from list",
"simulated patterns")
explan <- if(dual) paren(paste(nsim2, "to estimate the mean and",
nsim, "to calculate envelopes")) else ""
cat(paste(action, Nsim, descrip, explan, "...\n"))
}
# determine whether simulated point patterns should be saved
catchpatterns <- savepatterns && simtype != "list"
Caughtpatterns <- list()
# allocate space for computed function values
nrvals <- length(rvals)
simvals <- matrix(, nrow=nrvals, ncol=Nsim)
# arguments for function
funargs <-
resolve.defaults(if(rgiven) NULL else list(r=rvals),
list(...),
if(usecorrection) list(correction="best") else NULL)
# start simulation loop
nerr <- 0
for(i in 1:Nsim) {
ok <- FALSE
# safely generate a random pattern and apply function
while(!ok) {
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)),
expr=stop(paste("Evaluating the expression", sQuote("simulate"),
"did not yield an", Xobjectname)),
list=stop(paste("Internal error: list entry was not an",
Xobjectname)),
stop(paste("Internal error:", Xobjectname, "not generated"))
)
if(catchpatterns)
Caughtpatterns[[i]] <- Xsim
# apply function safely
funXsim <- try(do.call(fun, append(list(Xsim), funargs)))
ok <- !inherits(funXsim, "try-error")
if(!ok) {
nerr <- nerr + 1
if(nerr > maxnerr)
stop("Exceeded maximum number of errors")
cat("[retrying]\n")
}
}
# sanity checks
if(!inherits(funXsim, "fv"))
stop(paste("When applied to a simulated pattern, the function",
fname, "did not return an object of class",
sQuote("fv")))
argname.sim <- fvnames(funXsim, ".x")
valname.sim <- fvnames(funXsim, ".y")
if(argname.sim != argname)
stop(paste("The objects returned by", fname,
"when applied to a simulated pattern",
"and to the data pattern",
"are incompatible. They have different argument names",
sQuote(argname.sim), "and", sQuote(argname),
"respectively"))
if(valname.sim != valname)
stop(paste("When", fname, "is applied to a simulated pattern",
"it provides an estimate named", sQuote(valname.sim),
"whereas the estimate for the data pattern is named",
sQuote(valname),
". Try using the argument", sQuote("correction"),
"to make them compatible"))
if(tran) {
# extract only the recommended value
if(csr.theo)
funXsim <- funXsim[, c(argname, valname, "theo")]
else
funXsim <- funXsim[, c(argname, valname)]
# apply the transformation to it
funXsim <- eval(transform.funXsim)
}
# extract the values for simulation i
simvals.i <- funXsim[[valname]]
if(length(simvals.i) != nrvals)
stop("Vectors of function values have incompatible lengths")
simvals[ , i] <- funXsim[[valname]]
if(verbose)
progressreport(i, Nsim)
}
## end simulation loop
if(verbose) {
cat("\nDone.\n")
flush.console()
}
# ...........................................................
# save functions and/or patterns if so commanded
if(savefuns) {
alldata <- cbind(rvals, simvals)
simnames <- paste("sim", 1:nsim, sep="")
colnames(alldata) <- c("r", simnames)
alldata <- as.data.frame(alldata)
SimFuns <- fv(alldata,
argu="r",
ylab=attr(funX, "ylab"),
valu="sim1",
fmla= deparse(. ~ r),
alim=attr(funX, "alim"),
labl=names(alldata),
desc=c("distance argument r",
paste("Simulation ", 1:nsim, sep="")))
fvnames(SimFuns, ".") <- simnames
}
if(savepatterns)
SimPats <- if(simtype == "list") SimDataList else Caughtpatterns
######### COMPUTE ENVELOPES #######################
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, csr=csr, use.theory=csr.theo,
nrank=nrank, ginterval=ginterval, nSD=nSD,
Yname=Yname)
# tack on envelope information
attr(result, "einfo") <- EnvelopeInfo
# tack on functions and/or patterns if so commanded
if(savefuns)
attr(result, "simfuns") <- SimFuns
if(savepatterns) {
attr(result, "simpatterns") <- SimPats
attr(result, "datapattern") <- X
}
return(result)
}
plot.envelope <- function(x, ..., shade=c("hi", "lo")) {
xname <- deparse(substitute(x))
argh <- list(...)
if(missing(shade) &&
length(argh) > 0 &&
any(isfo <- unlist(lapply(argh, inherits, what="formula")))) {
# the plot is specified by a formula;
# check whether columns 'hi' and 'lo' are used
fmla <- argh[[min(which(isfo))]]
vars <- variablesinformula(fmla)
hilo <- c("hi", "lo")
dotnames <- fvnames(x, ".")
present <- (all(hilo %in% vars) ||
(all(hilo %in% dotnames) && ("." %in% vars)))
if(!present)
shade <- NULL
}
do.call("plot.fv", resolve.defaults(list(x), list(...),
list(main=xname, shade=shade)))
}
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
nr <- e$nrank
nsim <- e$nsim
V <- e$VARIANCE
fname <- deparse(attr(x, "ylab"))
type <- if(V) paste("Pointwise", e$nSD, "sigma") else
if(g) "Simultaneous" else "Pointwise"
cat(paste(type, "critical envelopes for", fname, "\n"))
if(!is.null(valname <- e$valname))
cat(paste("Edge correction:",
dQuote(valname), "\n"))
# determine *actual* type of simulation
descrip <-
if(csr) "simulations of CSR"
else if(!is.null(simtype)) {
switch(simtype,
csr="simulations of CSR",
rmh=paste("simulations of fitted",
if(pois) "Poisson" else "Gibbs",
"model"),
kppm="simulations of fitted cluster model",
expr="evaluations of user-supplied expression",
list="point pattern datasets in user-supplied list",
funs="columns of user-supplied data")
} else "simulations of fitted model"
#
cat(paste("Obtained from", nsim, descrip, "\n"))
#
if(!is.null(e$dual) && e$dual)
cat(paste("Theoretical (i.e. null) mean value of", fname,
"estimated from a separate set of",
e$nsim2, "simulations\n"))
if(!is.null(attr(x, "simfuns")))
cat("(All simulated function values are stored)\n")
if(!is.null(attr(x, "simpatterns")))
cat("(All simulated point patterns are stored)\n")
alpha <- if(g) { nr/(nsim+1) } else { 2 * nr/(nsim+1) }
if(!V) {
# significance interpretation!
cat(paste("Significance level of",
if(!g) "pointwise",
"Monte Carlo test:",
paste(if(g) nr else 2 * nr,
"/", nsim+1, sep=""),
"=", alpha, "\n"))
}
cat(paste("Data:", e$Yname, "\n"))
NextMethod("print")
}
summary.envelope <- function(object, ...) {
e <- attr(object, "einfo")
g <- e$global
nr <- e$nrank
nsim <- e$nsim
csr <- e$csr
pois <- e$pois
if(is.null(pois)) pois <- csr
has.theo <- "theo" %in% fvnames(object, "*")
csr.theo <- csr && has.theo
simtype <- e$simtype
fname <- deparse(attr(object, "ylab"))
V <- e$VARIANCE
type <- if(V) paste("Pointwise", e$nSD, "sigma") else
if(g) "Simultaneous" else "Pointwise"
cat(paste(type, "critical envelopes for", fname, "\n"))
# determine *actual* type of simulation
descrip <-
if(csr) "simulations of CSR"
else if(!is.null(simtype)) {
switch(simtype,
csr="simulations of CSR",
rmh=paste("simulations of fitted",
if(pois) "Poisson" else "Gibbs",
"model"),
kppm="simulations of fitted cluster model",
expr="evaluations of user-supplied expression",
list="point pattern datasets in user-supplied list",
funs="columns of user-supplied data",
"simulated point patterns")
} else "simulations of fitted model"
#
cat(paste("Obtained from", nsim, descrip, "\n"))
#
if(!is.null(attr(object, "simfuns")))
cat(paste("(All", nsim, "simulated function values",
"are stored in attr(,", dQuote("simfuns"), ") )\n"))
if(!is.null(attr(object, "simpatterns")))
cat(paste("(All", nsim, "simulated point patterns",
"are stored in attr(,", dQuote("simpatterns"), ") )\n"))
#
if(V) {
# nSD envelopes
cat(paste("Envelopes computed as sample mean plus/minus",
e$nSD, "sample standard deviations\n"))
} else {
# critical envelopes
lo.ord <- if(nr == 1) "minimum" else paste(ordinal(nr), "smallest")
hi.ord <- if(nr == 1) "maximum" else paste(ordinal(nsim-nr+1), "largest")
if(g)
cat(paste("Envelopes computed as",
if(csr.theo) "theoretical curve" else "mean of simulations",
"plus/minus", hi.ord,
"simulated value of maximum absolute deviation\n"))
else {
cat(paste("Upper envelope: pointwise", hi.ord,"of simulated curves\n"))
cat(paste("Lower envelope: pointwise", lo.ord,"of simulated curves\n"))
}
alpha <- if(g) { nr/(nsim+1) } else { 2 * nr/(nsim+1) }
cat(paste("Significance level of Monte Carlo test:",
paste(if(g) nr else 2 * nr,
"/", nsim+1, sep=""),
"=", alpha, "\n"))
}
cat(paste("Data:", e$Yname, "\n"))
return(invisible(NULL))
}
# envelope.matrix
# core functionality to compute envelope values
# theory = funX[["theo"]]
# observed = fX
envelope.matrix <- function(Y, ...,
rvals=NULL, observed=NULL, theory=NULL,
funX=NULL,
nsim=NULL, nsim2=NULL,
jsim=NULL, jsim.mean=NULL,
type=c("pointwise", "global", "variance"),
csr=FALSE, use.theory = csr,
nrank=1, ginterval=NULL, nSD=2,
check=TRUE,
Yname=NULL) {
if(is.null(Yname))
Yname <- deparse(substitute(Y), 60, nlines=1)
type <- match.arg(type)
if(!is.null(funX)) stopifnot(is.fv(funX))
if(is.null(rvals) && is.null(observed) && !is.null(funX)) {
# assume funX is summary function for observed data
rvals <- with(funX, .x)
observed <- with(funX, .y)
theory <- if(use.theory && "theo" %in% names(funX)) with(funX, theo) else NULL
} else if(check) {
# validate vectors of data
if(is.null(rvals)) stop("rvals must be supplied")
if(is.null(observed)) stop("observed must be supplied")
stopifnot(length(rvals) == nrow(Y))
stopifnot(length(observed) == length(rvals))
}
if(use.theory) {
use.theory <- !is.null(theory)
if(use.theory && check) stopifnot(length(theory) == length(rvals))
}
atr <- if(!is.null(funX)) attributes(funX) else
list(alim=range(rvals),
ylab=quote(f(r)),
yexp=quote(f(r)),
fname="f")
simvals <- Y
fX <- observed
# determine numbers of columns used
Ncol <- ncol(simvals)
if(Ncol < 2)
stop("Need at least 2 columns of function values")
if(is.null(jsim) && !is.null(nsim)) {
# usual case - 'nsim' determines 'jsim'
if(nsim > Ncol)
stop(paste(nsim, "simulations are not available; only",
Ncol, "columns provided"))
jsim <- 1:nsim
if(!is.null(nsim2)) {
# 'nsim2' determines 'jsim.mean'
if(nsim + nsim2 > Ncol)
stop(paste(nsim, "+", nsim2, "=", nsim+nsim2,
"simulations are not available; only",
Ncol, "columns provided"))
jsim.mean <- nsim + 1:nsim2
}
}
restrict.columns <- !is.null(jsim)
dual <- !is.null(jsim.mean)
switch(type,
pointwise = {
# ....... POINTWISE ENVELOPES ...............................
simvals[is.infinite(simvals)] <- NA
if(restrict.columns)
simvals <- simvals[,jsim]
nsim <- ncol(simvals)
nsim.mean <- NULL
if(nrank == 1) {
lohi <- apply(simvals, 1, range)
} else {
lohi <- apply(simvals, 1,
function(x, n) { sort(x)[n] },
n=c(nrank, nsim-nrank+1))
}
lo <- lohi[1,]
hi <- lohi[2,]
lo.name <- paste("lower pointwise envelope of %s from simulations")
hi.name <- paste("upper pointwise envelope of %s from simulations")
#
if(use.theory) {
results <- data.frame(r=rvals,
obs=fX,
theo=theory,
lo=lo,
hi=hi)
} else {
m <- apply(simvals, 1, mean, na.rm=TRUE)
results <- data.frame(r=rvals,
obs=fX,
mmean=m,
lo=lo,
hi=hi)
}
},
global = {
# ..... SIMULTANEOUS ENVELOPES ..........................
if(!is.null(ginterval)) {
domain <- (rvals >= ginterval[1]) & (rvals <= ginterval[2])
funX <- funX[domain, ]
simvals <- simvals[domain, ]
} else domain <- rep(TRUE, length(rvals))
simvals[is.infinite(simvals)] <- NA
if(use.theory) {
reference <- theory[domain]
if(restrict.columns)
simvals <- simvals[, jsim]
nsim.mean <- NULL
} else if(dual) {
# Estimate the mean from one set of columns
# Form envelopes from another set of columns
simvals.mean <- simvals[, jsim.mean]
reference <- mmean <- apply(simvals.mean, 1, mean, na.rm=TRUE)
nsim.mean <- ncol(simvals.mean)
simvals <- simvals[, jsim]
} else {
# Estimate the mean and form the envelopes using the same data
if(restrict.columns)
simvals <- simvals[, jsim]
reference <- mmean <- apply(simvals, 1, mean, na.rm=TRUE)
nsim.mean <- NULL
}
nsim <- ncol(simvals)
# compute max absolute deviations
deviations <- sweep(simvals, 1, reference)
suprema <- apply(abs(deviations), 2, max, na.rm=TRUE)
# ranked deviations
dmax <- sort(suprema)[nsim-nrank+1]
# simultaneous bands
lo <- reference - dmax
hi <- reference + dmax
lo.name <- "lower critical boundary for %s"
hi.name <- "upper critical boundary for %s"
if(use.theory)
results <- data.frame(r=rvals[domain],
obs=fX[domain],
theo=reference,
lo=lo,
hi=hi)
else
results <- data.frame(r=rvals[domain],
obs=fX[domain],
mmean=reference,
lo=lo,
hi=hi)
},
variance={
# ....... POINTWISE MEAN, VARIANCE etc ......................
simvals[is.infinite(simvals)] <- NA
if(restrict.columns)
simvals <- simvals[, jsim]
nsim <- ncol(simvals)
nsim.mean <- NULL
Ef <- apply(simvals, 1, mean, na.rm=TRUE)
varf <- apply(simvals, 1, var, na.rm=TRUE)
# derived quantities
sd <- sqrt(varf)
stdres <- (fX-Ef)/sd
stdres[!is.finite(stdres)] <- NA
# critical limits
lo <- Ef - nSD * sd
hi <- Ef + nSD * sd
lo.name <- paste("lower", nSD, "sigma critical limit for %s")
hi.name <- paste("upper", nSD, "sigma critical limit for %s")
# confidence interval
loCI <- Ef - nSD * sd/sqrt(nsim)
hiCI <- Ef + nSD * sd/sqrt(nsim)
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")
# put together
if(use.theory) {
results <- data.frame(r=rvals,
obs=fX,
theo=theory,
lo=lo,
hi=hi)
morestuff <- data.frame(mmean=Ef,
var=varf,
res=fX-Ef,
stdres=stdres,
loCI=loCI,
hiCI=hiCI)
mslabl <- c("bar(%s)(r)",
"paste(var,%s)(r)",
"paste(res,%s)(r)",
"paste(stdres,%s)(r)",
"%s[loCI](r)", "%s[hiCI](r)")
msdesc <- c("sample mean of %s from simulations",
"sample variance of %s from simulations",
"raw residual",
"standardised residual",
loCI.name, hiCI.name)
} else {
results <- data.frame(r=rvals,
obs=fX,
mmean=Ef,
lo=lo,
hi=hi)
morestuff <- data.frame(var=varf,
res=fX-Ef,
stdres=stdres,
loCI=loCI,
hiCI=hiCI)
mslabl <- c("paste(var,%s)(r)",
"paste(res,%s)(r)",
"paste(stdres,%s)(r)",
"%s[loCI](r)", "%s[hiCI](r)")
msdesc <- c("sample variance of %s from simulations",
"raw residual",
"standardised residual",
loCI.name, hiCI.name)
}
}
)
############ WRAP UP #########################
if(use.theory) {
# reference is computed curve `theo'
reflabl <- "%s[theo](r)"
refdesc <- "theoretical value of %s"
if(csr) refdesc <- paste(refdesc, "for CSR")
} else {
# reference is sample mean of simulations
reflabl <- "bar(%s)(r)"
refdesc <- "sample mean of %s from simulations"
}
result <- fv(results,
argu="r",
ylab=atr$ylab,
valu="obs",
fmla= deparse(. ~ r),
alim=atr$alim,
labl=c("r", "%s[obs](r)",
reflabl,
"%s[lo](r)", "%s[hi](r)"),
desc=c("distance argument r",
"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
unitname(result) <- unitname(funX)
class(result) <- c("envelope", class(result))
# tack on envelope information
attr(result, "einfo") <- list(global = (type =="global"),
csr = csr,
simtype = "funs",
nrank = nrank,
nsim = nsim,
VARIANCE = (type == "variance"),
nSD = nSD,
valname = NULL,
dual = dual,
nsim = nsim,
nsim2 = nsim.mean,
Yname = Yname)
return(result)
}
envelope.envelope <- function(Y, fun=NULL, ...) {
Yname <- deparse(substitute(Y), 60, nlines=1)
stopifnot(inherits(Y, "envelope"))
csr <- attr(Y, "einfo")$csr
if(!is.null(fun)) {
# apply new function
# point patterns are required
sp <- attr(Y, "simpatterns")
X <- attr(Y, "datapattern")
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"))
result <- do.call(envelope,
resolve.defaults(list(Y=X, fun=fun, simulate=sp),
list(...),
list(Yname=Yname,
nsim=length(sp)),
.StripNull=TRUE))
copyacross <- c("Yname", "csr", "csr.theo", "simtype")
attr(result, "einfo")[copyacross] <- attr(Y, "einfo")[copyacross]
return(result)
}
# compute new envelope with existing simulated functions
sf <- attr(Y, "simfuns")
if(is.null(sf))
stop(paste("Y does not contain a", dQuote("simfuns"), "attribute",
"(it was not generated with savefuns=TRUE)"))
# extract simulated function values
df <- as.data.frame(sf)
rname <- fvnames(sf, ".x")
df <- df[, (names(df) != rname)]
result <- do.call(envelope.matrix,
resolve.defaults(list(Y=as.matrix(df)),
list(...),
list(funX=Y,
Yname=Yname),
.StripNull=TRUE))
copyacross <- c("Yname", "csr", "csr.theo", "simtype")
attr(result, "einfo")[copyacross] <- attr(Y, "einfo")[copyacross]
return(result)
}
pool <- function(...) {
UseMethod("pool")
}
pool.envelope <- function(...) {
Yname <- paste(deparse(sys.call()), collapse="")
if(nchar(Yname) > 60) Yname <- paste(substr(Yname, 1, 40), "[..]")
Elist <- list(...)
nE <- length(Elist)
if(nE == 0) return(NULL)
# validate....
# 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)
}
if(nE == 1) return(Elist[[1]])
# All arguments must have 'simfuns'
SFlist <- lapply(Elist, attr, which="simfuns")
isnul <- unlist(lapply(SFlist, is.null))
if(any(isnul)) {
n <- sum(isnul)
why <- paste(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)
}
# vectors of r values must be identical
rlist <- lapply(SFlist, function(z) { with(z, .x) })
rvals <- rlist[[1]]
samer <- unlist(lapply(rlist, identical, y=rvals))
if(!all(samer))
stop(paste("Simulated function values are not compatible",
"(different values of function argument)"))
# functions must be the same
fnames <- unique(unlist(lapply(SFlist, attr, which="fname")))
if(length(fnames) > 1)
stop(paste("Envelopes were computed for different functions:",
commasep(sQuote(fnames))))
# ... reconcile parameters .......
einfolist <- lapply(Elist, attr, which="einfo")
global <- unique(unlist(lapply(einfolist, function(z) { z$global } )))
VARIANCE <- unique(unlist(lapply(einfolist, function(z) { z$VARIANCE } )))
simtype <- unique(unlist(lapply(einfolist, function(z) { z$simtype } )))
csr <- unique(unlist(lapply(einfolist, function(z) { z$csr } )))
csr.theo <- unique(unlist(lapply(einfolist, function(z) { z$csr.theo } )))
resolve <- function(x, x0, warn) {
xname <- deparse(substitute(x))
if(length(x) == 1)
return(x)
if(missing(warn))
warn <- paste("Envelopes were generated using different values",
"of argument", paste(sQuote(xname), ";", sep=""),
"reverting to default value")
if(!is.null(warn))
warning(warn, call.=FALSE)
return(x0)
}
global <- resolve(global, FALSE)
VARIANCE <- resolve(VARIANCE, FALSE)
simtype <- resolve(simtype, "funs",
"Envelopes were generated using different types of simulation")
csr <- resolve(csr, FALSE, NULL)
csr.theo <- resolve(csr.theo, FALSE, NULL)
type <- if(global) "global" else if(VARIANCE) "variance" else "pointwise"
# ..... ready to compute
getsimvals <- function(z) {
rname <- fvnames(z, ".x")
as.matrix(as.data.frame(z)[, names(z) != rname])
}
matlist <- lapply(SFlist, getsimvals)
bigmat <- do.call(cbind, matlist)
result <- envelope(bigmat, funX=Elist[[1]], type=type, csr=csr, Yname=Yname)
return(result)
}
Computing file changes ...
