https://github.com/cran/spatstat
Tip revision: ace26c246ee6feb8779515fa668bec59b24a1fcc authored by Adrian Baddeley on 12 March 2007, 13:35:27 UTC
version 1.11-2
version 1.11-2
Tip revision: ace26c2
lurking.R
# Lurking variable plot for arbitrary covariate.
#
#
# $Revision: 1.12 $ $Date: 2007/02/16 07:08:38 $
#
lurking <- function(object, covariate, type="eem",
cumulative=TRUE,
clipwindow=default.clipwindow(object),
rv = NULL,
plot.sd, plot.it=TRUE,
typename,
covname, oldstyle=FALSE, ...) {
# validate
verifyclass(object, "ppm")
# extract spatial locations
Q <- quad.ppm(object)
datapoints <- Q$data
quadpoints <- union.quad(Q)
Z <- is.data(Q)
wts <- w.quad(Q)
#################################################################
# compute the covariate
if(is.vector(covariate) && is.numeric(covariate))
covvalues <- covariate
else if(is.im(covariate))
covvalues <- covariate[quadpoints]
else if(is.expression(covariate)) {
# Evaluate the desired covariate at all quadrature points
# Set up environment for evaluating expression
glmdata <- object$internal$glmdata
# Fix special cases
if(is.null(glmdata)) {
# default
glmdata <- data.frame(x=quadpoints$x, y=quadpoints$y)
if(is.marked(quadpoints))
glmdata$marks <- marks(quadpoints)
}
# ensure x and y are in data frame
if(!all(c("x","y") %in% names(glmdata))) {
glmdata$x <- quadpoints$x
glmdata$y <- quadpoints$y
}
# Evaluate expression
evaluate <- function(a, b) {
if (exists("is.R") && is.R())
eval(a, envir = b)
else eval(a, local = b)
}
covvalues <- evaluate(covariate, glmdata)
if(!is.numeric(covvalues))
stop("The evaluated covariate is not numeric")
} else
stop(paste("The", sQuote("covariate"),
"should be either a numeric vector or an expression"))
# Validate covariate values
if(any(is.na(covvalues)))
stop("covariate contains NA's")
if(any(is.infinite(covvalues) | is.nan(covvalues)))
stop("covariate contains Inf or NaN values")
if(length(covvalues) != quadpoints$n)
stop("Length of covariate vector,", length(covvalues), "!=",
quadpoints$n, ", number of quadrature points")
# Quadrature points marked by covariate value
covq <- quadpoints %mark% covvalues
################################################################
# Residuals/marks attached to appropriate locations.
# Stoyan-Grabarnik weights are attached to the data points only.
# Others (residuals) are attached to all quadrature points.
resvalues <-
if(!is.null(rv)) rv
else if(type=="eem") eem(object)
else residuals.ppm(object,type=type)
res <- (if(type == "eem") datapoints else quadpoints) %mark% resvalues
# ... and the same locations marked by the covariate
covres <- if(type == "eem") covq[Z] else covq
# NAMES OF THINGS
# name of the covariate
if(missing(covname))
covname <- if(is.expression(covariate)) paste(covariate) else "covariate"
# type of residual/mark
if(missing(typename))
typename <- if(!is.null(rv)) "rv" else attr(resvalues, "typename")
#######################################################################
# START ANALYSIS
# Clip to subwindow if needed
clip <- !is.poisson.ppm(object) ||
(!missing(clipwindow) && !is.null(clipwindow))
if(clip) {
covq <- covq[clipwindow]
res <- res[clipwindow]
covres <- covres[clipwindow]
clipquad <- inside.owin(quadpoints$x, quadpoints$y, clipwindow)
wts <- wts[ clipquad ]
}
# -----------------------------------------------------------------------
# (A) EMPIRICAL CUMULATIVE FUNCTION
# based on data points if type="eem", otherwise on quadrature points
# cumulative sums which ignore NA's
cumsumna <- function(x) {
x[is.na(x)] <- 0
return(cumsum(x))
}
# Reorder the data/quad points in order of increasing covariate value
# and then compute the cumulative sum of their residuals/marks
markscovres <- marks(covres)
o <- order(markscovres)
covsort <- markscovres[o]
cummark <- cumsumna(marks(res)[o])
# we'll plot(covsort, cummark) in the cumulative case
# (B) THEORETICAL MEAN CUMULATIVE FUNCTION
# based on all quadrature points
# Range of covariate values
covqmarks <- marks(covq)
covrange <- range(covqmarks, na.rm=TRUE)
# Suitable breakpoints
cvalues <- seq(covrange[1], covrange[2], length=100)
csmall <- cvalues[1] - diff(cvalues[1:2])
cbreaks <- c(csmall, cvalues)
# cumulative area as function of covariate values
covclass <- cut(covqmarks, breaks=cbreaks)
increm <- tapply(wts, covclass, sum)
cumarea <- cumsumna(increm)
# compute theoretical mean (when model is true)
mean0 <- if(type == "eem") cumarea else rep(0, length(cumarea))
# we'll plot(cvalues, mean0) in the cumulative case
# (A'),(B') DERIVATIVES OF (A) AND (B)
# Required if cumulative=FALSE
# Estimated by spline smoothing
if(!cumulative) {
# fit smoothing spline to (A)
ss <- smooth.spline(covsort, cummark, ...)
# estimate derivative of (A)
derivmark <- predict(ss, covsort, deriv=1)$y
# similarly for (B)
ss <- smooth.spline(cvalues, mean0, ...)
derivmean <- predict(ss, cvalues, deriv=1)$y
}
# -----------------------------------------------------------------------
# Store what will be plotted
if(cumulative) {
empirical <- data.frame(covariate=covsort, value=cummark)
theoretical <- data.frame(covariate=cvalues, mean=mean0)
} else {
empirical <- data.frame(covariate=covsort, value=derivmark)
theoretical <- data.frame(covariate=cvalues, mean=derivmean)
}
# ------------------------------------------------------------------------
# (C) STANDARD DEVIATION if desired
# (currently implemented only for Poisson)
# (currently implemented only for cumulative case)
if(missing(plot.sd))
plot.sd <- is.poisson.ppm(object)
else if(plot.sd && !is.poisson.ppm(object))
warning("standard deviation is calculated for Poisson model; not valid for this model")
if(plot.sd && cumulative) {
# Fitted intensity at quadrature points
lambda <- fitted.ppm(object, type="trend")
# Fisher information for coefficients
asymp <- vcov(object,what="internals")
Fisher <- asymp$fisher
# Local sufficient statistic at quadrature points
suff <- asymp$suff
# Clip if required
if(clip) {
lambda <- lambda[clipquad]
suff <- suff[clipquad, , drop=FALSE] # suff is a matrix
}
# First term: integral of lambda^(2p+1)
switch(type,
pearson={
varI <- cumarea
},
raw={
# Compute sum of w*lambda for quadrature points in each interval
dvar <- tapply(wts * lambda, covclass, sum)
# tapply() returns NA when the table is empty
dvar[is.na(dvar)] <- 0
# Cumulate
varI <- cumsum(dvar)
},
inverse=, # same as eem
eem={
# Compute sum of w/lambda for quadrature points in each interval
dvar <- tapply(wts / lambda, covclass, sum)
# tapply() returns NA when the table is empty
dvar[is.na(dvar)] <- 0
# Cumulate
varI <- cumsum(dvar)
})
# variance-covariance matrix of coefficients
V <- try(solve(Fisher), silent=TRUE)
if(inherits(V, "try-error")) {
warning("Fisher information is singular; reverting to oldstyle=TRUE")
oldstyle <- TRUE
}
# Second term: B' V B
if(oldstyle) {
varII <- 0
} else {
# lamp = lambda^(p + 1)
lamp <- switch(type,
raw = lambda,
pearson = sqrt(lambda),
inverse =,
eem = ifelse(lambda > 0, 1, 0))
# Compute sum of w * lamp * suff for quad points in intervals
Bcontrib <- as.vector(wts * lamp) * suff
dB <- matrix(, nrow=length(cumarea), ncol=ncol(Bcontrib))
for(j in seq(ncol(dB)))
dB[,j] <- tapply(Bcontrib[,j], covclass, sum)
# tapply() returns NA when the table is empty
dB[is.na(dB)] <- 0
# Cumulate columns
B <- apply(dB, 2, cumsum)
# compute B' V B for each i
varII <- diag(B %*% V %*% t(B))
}
#
# variance of residuals
varR <- varI - varII
# trap numerical errors
nbg <- (varR < 0)
if(any(nbg)) {
ran <- range(varR)
varR[nbg] <- 0
relerr <- abs(ran[1]/ran[2])
nerr <- sum(nbg)
if(relerr > 1e-6)
warning(paste(nerr, "negative",
ngettext(nerr, "value (", "values (min="),
signif(ran[1], 4), ")",
"of residual variance reset to zero",
"(out of", length(varR), "values)"))
}
theoretical$sd <- sqrt(varR)
}
# --------------- PLOT THEM ----------------------------------
if(plot.it) {
# work out plot range
mr <- range(c(0, empirical$value, theoretical$mean), na.rm=TRUE)
if(!is.null(theoretical$sd))
mr <- range(c(mr, theoretical$mean + 2 * theoretical$sd,
theoretical$mean - 2 * theoretical$sd),
na.rm=TRUE)
# start plot
plot(covrange, mr, type="n", xlab=covname,
ylab=paste(if(cumulative)"cumulative" else "marginal", typename))
# (A)/(A') Empirical
lines(value ~ covariate, empirical)
# (B)/(B') Theoretical mean
lines(mean ~ covariate, theoretical, lty=2)
# (C) Standard deviation
if(!is.null(theoretical$sd)) {
lines(mean + 2 * sd ~ covariate, theoretical, lty=3)
lines(mean - 2 * sd ~ covariate, theoretical, lty=3)
}
}
# ---------------- RETURN COORDINATES ----------------------------
stuff <- list(empirical=empirical, theoretical=theoretical)
return(invisible(stuff))
}