slrm.R
#
# slrm.R
#
# Spatial Logistic Regression
#
# $Revision: 1.12 $ $Date: 2012/02/04 11:46:38 $
#
slrm <- function(formula, ..., data=NULL, offset=TRUE, link="logit",
dataAtPoints=NULL, splitby=NULL) {
# remember call
CallInfo <- list(callstring = paste(deparse(sys.call()), collapse=""),
cl = match.call(),
formula = formula,
offset=offset,
link=link,
splitby=splitby,
dotargs=list(...))
if(!(link %in% c("logit", "cloglog")))
stop(paste("Unrecognised link", dQuote(link)))
########### INTERPRET FORMULA ##############################
if(!inherits(formula, "formula"))
stop(paste("Argument", dQuote("formula"), "should be a formula"))
# check formula has LHS and RHS. Extract them
if(length(formula) < 3)
stop(paste("Argument", sQuote("formula"),
"must have a left hand side"))
Yname <- lhs <- formula[[2]]
trend <- rhs <- formula[c(1,3)]
if(!is.name(Yname))
stop("Left hand side of formula should be a single name")
Yname <- paste(Yname)
if(!inherits(trend, "formula"))
stop("Internal error: failed to extract RHS of formula")
varnames <- unique(variablesinformula(trend))
specials <- c("x", "y", "logpixelarea")
covnames <- varnames[!(varnames %in% specials)]
# add 'splitby' to covariate names
if(!is.null(splitby)) {
if(!is.character(splitby) || length(splitby) != 1)
stop("splitby should be a single character string")
covnames <- unique(c(covnames, splitby))
}
CallInfo$responsename <- Yname
CallInfo$varnames <- varnames
CallInfo$covnames <- covnames
# Parent environment
parenv <- environment(formula)
######## FIND DATA AND RESHAPE #######################
Data <- slr.prepare(CallInfo, parenv, data, dataAtPoints, splitby)
W <- Data$W
df <- Data$df
######## FIT MODEL ###############################
dformula <- formula
if(offset) {
# insert offset term in formula
rhs <- paste(as.character(rhs), collapse=" ")
rhs <- paste(c(rhs, "offset(logpixelarea)"), collapse="+")
dformula <- as.formula(paste(Yname, rhs))
}
linkname <- link
FIT <- glm(dformula, family=binomial(link=linkname),
data=df)
result <- list(call = CallInfo$cl,
CallInfo = CallInfo,
Data = Data,
Fit = list(FIT=FIT, dformula=dformula),
terms = terms(formula))
class(result) <- c("slrm", class(result))
return(result)
}
################ UTILITY TO FIND AND RESHAPE DATA #################
slr.prepare <- function(CallInfo, envir, data,
dataAtPoints=NULL, splitby=NULL,
clip=TRUE) {
# CallInfo is produced by slrm()
# envir is parent environment of model formula
# data is 'data' argument that takes precedence over 'envir'
# 'clip' is TRUE if the data should be clipped to the domain of Y
Yname <- CallInfo$responsename
varnames <- CallInfo$varnames
covnames <- CallInfo$covnames
dotargs <- CallInfo$dotargs
#
getobj <- function(nama, env, dat) {
if(!is.null(dat) && !is.null(x <- dat[[nama]]))
return(x)
else return(get(nama, envir=env))
}
# Get the response point pattern Y
Y <- getobj(Yname, envir, data)
if(!is.ppp(Y))
stop(paste("The response", sQuote(Yname), "must be a point pattern"))
#
if(!is.null(dataAtPoints)) {
dataAtPoints <- as.data.frame(dataAtPoints)
if(nrow(dataAtPoints) != npoints(Y))
stop(paste("dataAtPoints should have one row for each point in",
dQuote(Yname)))
}
# Find the covariates
ncov <- length(covnames)
covlist <- lapply(as.list(covnames), getobj, env = envir, dat=data)
names(covlist) <- covnames
# Each covariate should be an image, a window, a function, or a single number
if(ncov == 0) {
isim <- isowin <- ismask <- isfun <- isnum <- isspatial <- israster <- logical(0)
} else {
isim <- unlist(lapply(covlist, is.im))
isowin <- unlist(lapply(covlist, is.owin))
ismask <- unlist(lapply(covlist, is.mask))
isfun <- unlist(lapply(covlist, is.function))
isspatial <- isim | isowin | isfun
israster <- isim | ismask
isnum <- unlist(lapply(covlist,
function(x) { is.numeric(x) && length(x) == 1} ))
}
if(!all(ok <- (isspatial | isnum))) {
n <- sum(!ok)
stop(paste(ngettext(n, "The argument", "Each of the arguments"),
commasep(sQuote(covnames[!ok])),
"should be either an image, a window, or a single number"))
}
# 'splitby'
if(!is.null(splitby)) {
splitwin <- covlist[[splitby]]
if(!is.owin(splitwin))
stop("The splitting covariate must be a window")
# ensure it is a polygonal window
covlist[[splitby]] <- splitwin <- as.polygonal(splitwin)
# delete splitting covariate from lists to be processed
issplit <- (covnames == splitby)
isspatial[issplit] <- FALSE
israster[issplit] <- FALSE
}
#
nnum <- sum(isnum)
nspatial <- sum(isspatial)
nraster <- sum(israster)
#
numlist <- covlist[isnum]
spatiallist <- covlist[isspatial]
rasterlist <- covlist[israster]
#
numnames <- names(numlist)
spatialnames <- names(spatiallist)
rasternames <- names(rasterlist)
#
######## CONVERT TO RASTER DATA ###############################
convert <- function(x,W) {
if(is.im(x) || is.function(x)) return(as.im(x,W))
if(is.owin(x)) return(as.im(x, W, value=TRUE, na.replace=FALSE))
return(NULL)
}
# determine spatial domain & common resolution: convert all data to it
if(length(dotargs) > 0 || nraster == 0) {
# Pixel resolution is determined by explicit arguments
if(clip) {
# Window extent is determined by response point pattern
D <- as.owin(Y)
} else {
# Window extent is union of domains of data
domains <- lapply(append(spatiallist, list(Y)), as.owin)
D <- do.call("union.owin", domains)
}
# Create template mask
W <- do.call("as.mask", append(list(D), dotargs))
# Convert all spatial objects to this resolution
spatiallist <- lapply(spatiallist, convert, W=W)
} else {
# Pixel resolution is determined implicitly by covariate data
W <- do.call("commonGrid", rasterlist)
if(clip) {
# Restrict data to spatial extent of response point pattern
W <- intersect.owin(W, as.owin(Y))
}
# Adjust spatial objects to this resolution
spatiallist <- lapply(spatiallist, convert, W=W)
}
# images containing coordinate values
xcoordim <- as.im(function(x,y){x}, W=W)
ycoordim <- as.im(function(x,y){y}, W=W)
#
# create a list of covariate images, with names as in formula
covimages <- append(list(x=xcoordim, y=ycoordim), spatiallist)
basepixelarea <- W$xstep * W$ystep
######## ASSEMBLE DATA FRAME ###############################
if(is.null(splitby)) {
df <- slrAssemblePixelData(Y, Yname, W,
covimages, dataAtPoints, basepixelarea)
sumYloga <- Y$n * log(basepixelarea)
} else {
# fractional pixel areas
pixsplit <- pixellate(splitwin, W)
splitpixelarea <- as.vector(as.matrix(pixsplit))
# determine which points of Y are inside/outside window
ins <- inside.owin(Y$x, Y$y, splitwin)
# split processing
dfIN <- slrAssemblePixelData(Y[ins], Yname, W, covimages,
dataAtPoints[ins, ], splitpixelarea)
dfIN[[splitby]] <- TRUE
dfOUT <- slrAssemblePixelData(Y[!ins], Yname, W, covimages,
dataAtPoints[!ins, ],
basepixelarea - splitpixelarea)
dfOUT[[splitby]] <- FALSE
df <- rbind(dfIN, dfOUT)
# sum of log pixel areas associated with points
Ysplit <- pixsplit[Y]
sumYloga <- sum(log(ifelse(ins, Ysplit, basepixelarea - Ysplit)))
}
# tack on any numeric values
df <- do.call("cbind", append(list(df), numlist))
### RETURN ALL
Data <- list(response=Y,
covariates=covlist,
spatialnames=spatialnames,
numnames=numnames,
W=W,
df=df,
sumYloga=sumYloga,
dataAtPoints=dataAtPoints)
return(Data)
}
#
slrAssemblePixelData <- function(Y, Yname, W,
covimages, dataAtPoints, pixelarea) {
# pixellate point pattern
Z <- pixellate(Y, W=W)
Z <- eval.im(as.integer(Z>0))
# overwrite pixel entries for data points using exact values
# coordinates
xcoordim <- covimages[["x"]]
ycoordim <- covimages[["y"]]
xcoordim[Y] <- Y$x
ycoordim[Y] <- Y$y
covimages[["x"]] <- xcoordim
covimages[["y"]] <- ycoordim
# overwrite pixel entries
if(!is.null(dataAtPoints)) {
enames <- colnames(dataAtPoints)
relevant <- enames %in% names(covimages)
for(v in enames[relevant]) {
cova <- covimages[[v]]
cova[Y] <- dataAtPoints[, v, drop=TRUE]
covimages[[v]] <- cova
}
}
# assemble list of all images
Ylist <- list(Z)
names(Ylist) <- Yname
allimages <- append(Ylist, covimages)
# extract pixel values of each image
pixelvalues <-
function(z) {
v <- as.vector(as.matrix(z))
if(z$type != "factor") return(v)
lev <- levels(z)
return(factor(v, levels=seq_along(lev), labels=lev))
}
pixdata <- lapply(allimages, pixelvalues)
df <- as.data.frame(pixdata)
# add log(pixel area) column
if(length(pixelarea) == 1) {
df <- cbind(df, logpixelarea=log(pixelarea))
} else {
ok <- (pixelarea > 0)
df <- cbind(df[ok, ], logpixelarea=log(pixelarea[ok]))
}
return(df)
}
is.slrm <- function(x) {
inherits(x, "slrm")
}
coef.slrm <- function(object, ...) {
coef(object$Fit$FIT)
}
print.slrm <- function(x, ...) {
lk <- x$CallInfo$link
switch(lk,
logit= {
cat("Fitted spatial logistic regression model\n")
},
cloglog= {
cat("Fitted spatial regression model (complementary log-log) \n")
},
{
cat("Fitted spatial regression model\n")
cat(paste("Link =", dQuote(lk), "\n"))
})
cat("Formula:\t")
print(x$CallInfo$formula)
cat("Fitted coefficients:\n")
print(coef(x))
return(invisible(NULL))
}
logLik.slrm <- function(object, ..., adjust=TRUE) {
FIT <- object$Fit$FIT
ll <- -deviance(FIT)/2
if(adjust) {
sumYloga <- object$Data$sumYloga
ll <- ll - sumYloga
}
attr(ll, "df") <- length(coef(object))
class(ll) <- "logLik"
return(ll)
}
fitted.slrm <- function(object, ...) {
if(length(list(...)) > 0)
warning("second argument (and any subsequent arguments) ignored")
predict(object, type="probabilities")
}
predict.slrm <- function(object, ..., type="intensity",
newdata=NULL, window=NULL) {
type <- pickoption("type", type,
c(probabilities="probabilities",
link="link",
intensity="intensity",
lambda="intensity"))
FIT <- object$Fit$FIT
link <- object$CallInfo$link
W <- object$Data$W
df <- object$Data$df
loga <- df$logpixelarea
if(is.null(newdata) && is.null(window)) {
# fitted values from existing fit
switch(type,
probabilities={
values <- fitted(FIT)
},
link={
values <- predict(FIT, type="link")
},
intensity={
# this calculation applies whether an offset was included or not
if(link == "cloglog") {
linkvalues <- predict(FIT, type="link")
values <- exp(linkvalues - loga)
} else {
probs <- fitted(FIT)
values <- -log(1-probs)/exp(loga)
}
}
)
v <- rep(NA, nrow(df))
which <- complete.cases(df)
v[which] <- values
out <- im(v, xcol=W$xcol, yrow=W$yrow, unitname=unitname(W))
return(out)
} else {
# prediction using new values
# update arguments that may affect pixel resolution
CallInfo <- object$CallInfo
CallInfo$dotargs <- resolve.defaults(list(...), CallInfo$dotargs)
#
if(!is.null(window)) {
# insert fake response in new window
if(is.null(newdata)) newdata <- list()
window <- as.owin(window)
newdata[[CallInfo$responsename]] <- ppp(numeric(0), numeric(0),
window=window)
}
# process new data
newData <- slr.prepare(CallInfo, environment(CallInfo$formula), newdata,
clip=is.null(window))
newdf <- newData$df
newW <- newData$W
newloga <- newdf$logpixelarea
# compute link values
linkvalues <- predict(FIT, newdata=newdf, type="link")
#
linkinv <- family(FIT)$linkinv
switch(type,
probabilities={
values <- linkinv(linkvalues)
},
link={
values <- linkvalues
},
intensity={
# this calculation applies whether an offset was included or not
if(link == "cloglog") {
values <- exp(linkvalues - newloga)
} else {
probs <- linkinv(linkvalues)
values <- -log(1-probs)/exp(newloga)
}
}
)
v <- rep(NA, nrow(newdf))
which <- complete.cases(newdf)
v[which] <- values
out <- im(v, xcol=newW$xcol, yrow=newW$yrow, unitname=unitname(W))
return(out)
}
}
plot.slrm <- function(x, ..., type="intensity") {
xname <- deparse(substitute(x))
y <- predict(x, type=type)
do.call("plot.im", resolve.defaults(list(x=y), list(...), list(main=xname)))
}
formula.slrm <- function(x, ...) {
f <- x$CallInfo$formula
return(f)
}
terms.slrm <- function(x, ...) {
terms(formula(x), ...)
}
labels.slrm <- function(object, ...) {
# extract fitted trend coefficients
co <- coef(object)
# model terms
tt <- terms(object)
lab <- attr(tt, "term.labels")
if(length(lab) == 0)
return(character(0))
# model matrix
mm <- model.matrix(object)
ass <- attr(mm, "assign")
# 'ass' associates coefficients with model terms
# except ass == 0 for the Intercept
coef.ok <- is.finite(co)
relevant <- (ass > 0)
okterms <- unique(ass[coef.ok & relevant])
return(lab[okterms])
}
extractAIC.slrm <- function (fit, scale = 0, k = 2, ...)
{
edf <- length(coef(fit))
aic <- AIC(fit)
c(edf, aic + (k - 2) * edf)
}
model.matrix.slrm <- function(object,..., keepNA=TRUE) {
FIT <- object$Fit$FIT
mm <- model.matrix(FIT, ...)
if(!keepNA)
return(mm)
df <- object$Data$df
if(nrow(mm) == nrow(df))
return(mm)
which <- complete.cases(df)
mmplus <- matrix(NA, nrow(df), ncol(mm))
mmplus[which, ] <- mm
return(mmplus)
}
update.slrm <- function(object, ..., evaluate=TRUE, env=parent.frame()) {
e <- update.default(object, ..., evaluate=FALSE)
if(evaluate)
e <- eval(e, envir=env)
return(e)
}
anova.slrm <- function(object, ..., test=NULL) {
objex <- append(list(object), list(...))
if(!all(unlist(lapply(objex, is.slrm))))
stop("Some arguments are not of class slrm")
fitz <- lapply(objex, function(z){z$Fit$FIT})
do.call("anova.glm", append(fitz, list(test=test)))
}
vcov.slrm <- function(object, ...) {
vcov(object$Fit$FIT)
}
unitname.slrm <- function(x) {
return(unitname(x$Data$response))
}
"unitname<-.slrm" <- function(x, value) {
unitname(x$Data$response) <- value
return(x)
}
is.stationary.slrm <- function(x) {
fo <- formula(x)
trend <- fo[c(1,3)]
return(identical.formulae(trend, ~1))
}
is.poisson.slrm <- function(x) { TRUE }
simulate.slrm <- function(object, nsim=1, seed=NULL, ...,
window=NULL, covariates=NULL,
verbose=TRUE) {
# .... copied from simulate.lm ....
if (!exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE))
runif(1)
if (is.null(seed))
RNGstate <- get(".Random.seed", envir = .GlobalEnv)
else {
R.seed <- get(".Random.seed", envir = .GlobalEnv)
set.seed(seed)
RNGstate <- structure(seed, kind = as.list(RNGkind()))
on.exit(assign(".Random.seed", R.seed, envir = .GlobalEnv))
}
# determine simulation window and compute intensity
if(!is.null(window))
stopifnot(is.owin(window))
lambda <- predict(object, type="intensity", newdata=covariates, window=window)
# max lambda (for efficiency)
summ <- summary(lambda)
lmax <- summ$max + 0.05 * diff(summ$range)
# run
out <- list()
if(verbose && (nsim > 1))
cat(paste("Generating", nsim, "simulations... "))
for(i in 1:nsim) {
out[[i]] <- rpoispp(lambda, lmax=lmax)
if(verbose) progressreport(i, nsim)
}
# pack up
out <- as.listof(out)
names(out) <- paste("Simulation", 1:nsim)
attr(out, "seed") <- RNGstate
return(out)
}