https://github.com/cran/spatstat
Tip revision: 3cabcc19f5c2170689ebcc8ae0a5f5dea2978dde authored by Adrian Baddeley on 26 April 2013, 07:28:25 UTC
version 1.31-2
version 1.31-2
Tip revision: 3cabcc1
scanstat.R
#
# scanstat.R
#
# Spatial scan statistics
#
# $Revision: 1.8 $ $Date: 2013/04/25 06:37:43 $
#
scanmeasure <- function(X, ...){
UseMethod("scanmeasure")
}
scanmeasure.ppp <- function(X, r, ..., method=c("counts", "fft")) {
method <- match.arg(method)
# enclosing window
W <- as.rectangle(as.owin(X))
# expand domain to include all circles
W <- grow.rectangle(W, r)
# determine pixel resolution
W <- as.mask(W, ...)
#
switch(method,
counts = {
# direct calculation using C code
# get new dimensions
dimyx <- W$dim
xr <- W$xrange
yr <- W$yrange
nr <- dimyx[1]
nc <- dimyx[2]
#
n <- npoints(X)
zz <- .C("scantrans",
x=as.double(X$x),
y=as.double(X$y),
n=as.integer(n),
xmin=as.double(xr[1]),
ymin=as.double(yr[1]),
xmax=as.double(xr[2]),
ymax=as.double(yr[2]),
nr=as.integer(nr),
nc=as.integer(nc),
R=as.double(r),
counts=as.integer(numeric(prod(dimyx))),
PACKAGE="spatstat")
zzz <- matrix(zz$counts, nrow=dimyx[1], ncol=dimyx[2], byrow=TRUE)
Z <- im(zzz, xrange=xr, yrange=yr, unitname=unitname(X))
},
fft = {
# Previous version of scanmeasure.ppp had
# Y <- pixellate(X, ..., padzero=TRUE)
# but this is liable to Gibbs phenomena.
# Instead, convolve with small Gaussian (sd = 1 pixel width)
sigma <- with(W, unique(c(xstep, ystep)))
Y <- density(X, ..., sigma=sigma)
# invoke scanmeasure.im
Z <- scanmeasure(Y, r)
Z <- eval.im(as.integer(round(Z)))
})
return(Z)
}
scanmeasure.im <- function(X, r, ...) {
D <- disc(radius=r)
eps <- with(X, c(xstep,ystep))
D <- as.im(as.mask(D, eps=eps))
Z <- imcov(X, D)
return(Z)
}
scanPoisLRTS <- function(nZ, nG, muZ, muG, alternative) {
nZco <- nG - nZ
muZco <- muG - muZ
nlogn <- function(n, a) ifelse(n == 0, 0, n * log(n/a))
ll <- nlogn(nZ, muZ) + nlogn(nZco, muZco) - nlogn(nG, muG)
criterion <- (nZ * muZco - muZ * nZco)
switch(alternative,
less={
ll[criterion > 0] <- 0
},
greater={
ll[criterion < 0] <- 0
},
two.sided={})
return(2 * ll)
}
scanBinomLRTS <- function(nZ, nG, muZ, muG, alternative) {
nZco <- nG - nZ
muZco <- muG - muZ
nlogn <- function(n, a) ifelse(n == 0, 0, n * log(n/a))
logbin <- function(k, n) { nlogn(k, n) + nlogn(n-k, n) }
ll <- logbin(nZ, muZ) + logbin(nZco, muZco) - logbin(nG, muG)
criterion <- (nZ * muZco - muZ * nZco)
switch(alternative,
less={
ll[criterion > 0] <- 0
},
greater={
ll[criterion < 0] <- 0
},
two.sided={})
return(2 * ll)
}
scanLRTS <- function(X, r, ...,
method=c("poisson", "binomial"),
baseline=NULL,
case=2,
alternative=c("greater", "less", "two.sided")) {
stopifnot(is.ppp(X))
method <- match.arg(method)
alternative <- match.arg(alternative)
switch(method,
poisson={
Y <- X
Xmask <- as.mask(as.owin(X), ...)
if(is.null(baseline)) {
mu <- as.im(Xmask, value=1)
} else if(is.ppm(baseline)) {
if(is.marked(baseline))
stop("baseline is a marked point process: not supported")
mu <- predict(baseline, locations=Xmask)
} else if(is.im(baseline) || is.function(baseline)) {
mu <- as.im(baseline, W=Xmask)
} else stop(paste("baseline should be",
"a pixel image, a function, or a fitted model"))
},
binomial={
stopifnot(is.multitype(X))
lev <- levels(marks(X))
if(length(lev) != 2)
warning("X should usually be a bivariate (2-type) point pattern")
if(!is.null(baseline))
stop("baseline is not supported in the binomial case")
if(is.character(case) && !(case %in% lev))
stop(paste("Unrecognised label for cases:", sQuote(case)))
if(is.numeric(case) && !(case %in% seq_along(lev)))
stop(paste("Undefined level:", case))
Y <- split(X)[[case]]
mu <- unmark(X)
})
nZ <- scanmeasure(Y, r, ...)
muZ <- scanmeasure(mu, r)
if(!compatible.im(nZ, muZ)) {
ha <- harmonise.im(nZ, muZ)
nZ <- ha[[1]]
muZ <- ha[[2]]
}
nG <- npoints(Y)
switch(method,
poisson = {
muG <- integral.im(mu)
result <- eval.im(scanPoisLRTS(nZ, nG, muZ, muG, alternative))
},
binomial = {
muG <- npoints(mu)
result <- eval.im(scanBinomLRTS(nZ, nG, muZ, muG, alternative))
},
{ result <- NULL })
return(result)
}
scan.test <- function(X, r, ...,
method=c("poisson", "binomial"),
nsim = 19,
baseline=NULL,
case = 2,
alternative=c("greater", "less", "two.sided"),
verbose=TRUE) {
dataname <- short.deparse(substitute(X))
stopifnot(is.ppp(X))
method <- match.arg(method)
alternative <- match.arg(alternative)
check.1.real(r)
check.1.real(nsim)
if(!(round(nsim) == nsim && nsim > 1))
stop("nsim should be an integer > 1")
regionname <- paste("circles of radius", r)
#
# compute observed loglikelihood function
# This also validates the arguments.
obsLRTS <- scanLRTS(X=X, r=r,
method=method,
alternative=alternative, baseline=baseline,
case=case, ...)
obs <- max(obsLRTS)
sim <- numeric(nsim)
# determine how to simulate
switch(method,
binomial={
methodname <- c("Spatial scan test",
"Null hypothesis: constant relative risk",
paste("Candidate cluster regions:", regionname),
"Likelihood: binomial",
paste("Monte Carlo p-value based on",
nsim, "simulations"))
lev <- levels(marks(X))
names(lev) <- lev
casename <- lev[case]
counted <- paste("points with mark", sQuote(casename),
"inside cluster region")
simexpr <- expression(rlabel(X))
},
poisson={
counted <- paste("points inside cluster region")
X <- unmark(X)
Xwin <- as.owin(X)
Xmask <- as.mask(Xwin, ...)
if(is.null(baseline)) {
nullname <- "Complete Spatial Randomness (CSR)"
lambda <- summary(X)$intensity
simexpr <- expression(runifpoispp(lambda, Xwin))
} else if(is.ppm(baseline)) {
nullname <- baseline$callstring
rmhstuff <- rmh(baseline, preponly=TRUE, verbose=FALSE)
simexpr <- expression(rmhEngine(rmhstuff))
} else if(is.im(baseline) || is.function(baseline)) {
nullname <- "Poisson process with intensity proportional to baseline"
base <- as.im(baseline, W=Xmask)
alpha <- npoints(X)/integral.im(base)
lambda <- eval.im(alpha * base)
simexpr <- expression(rpoispp(lambda))
} else stop(paste("baseline should be",
"a pixel image, a function, or a fitted model"))
methodname <- c("Spatial scan test",
paste("Null hypothesis:", nullname),
paste("Candidate cluster regions:", regionname),
"Likelihood: Poisson",
paste("Monte Carlo p-value based on",
nsim, "simulations"))
})
if(verbose) cat("Simulating...")
for(i in 1:nsim) {
if(verbose) progressreport(i, nsim)
Xsim <- eval(simexpr)
simLRTS <- scanLRTS(X=Xsim, r=r,
method=method, alternative=alternative,
baseline=baseline,
case=case,
...)
sim[i] <- max(simLRTS)
}
pval <- mean(c(sim,obs) >= obs, na.rm=TRUE)
names(obs) <- "maxLRTS"
nm.alternative <- switch(alternative,
greater="Excess of",
less="Deficit of",
two.sided="Two-sided: excess or deficit of",
stop("Unknown alternative"))
nm.alternative <- paste(nm.alternative, counted)
result <- list(statistic = obs,
p.value = pval,
alternative = nm.alternative,
method = methodname,
data.name = dataname)
class(result) <- c("scan.test", "htest")
attr(result, "obsLRTS") <- obsLRTS
attr(result, "X") <- X
return(result)
}
plot.scan.test <- function(x, ..., do.window=TRUE) {
xname <- short.deparse(substitute(x))
Z <- as.im(x)
do.call("plot",
resolve.defaults(list(x=Z), list(...), list(main=xname)))
if(do.window) {
X <- attr(x, "X")
plot(as.owin(X), add=TRUE)
}
invisible(NULL)
}
as.im.scan.test <- function(X, ...) {
X <- attr(X, "obsLRTS")
return(as.im(X, ...))
}