lppm.R
#
# lppm.R
#
# Point process models on a linear network
#
# $Revision: 1.44 $ $Date: 2018/05/02 02:04:45 $
#
lppm <- function(X, ...) {
UseMethod("lppm")
}
lppm.formula <- function(X, interaction=NULL, ..., data=NULL) {
## remember call
callstring <- paste(short.deparse(sys.call()), collapse = "")
cl <- match.call()
########### INTERPRET FORMULA ##############################
if(!inherits(X, "formula"))
stop(paste("Argument 'X' should be a formula"))
formula <- X
if(spatstat.options("expand.polynom"))
formula <- expand.polynom(formula)
## check formula has LHS and RHS. Extract them
if(length(formula) < 3)
stop(paste("Formula must have a left hand side"))
Yexpr <- formula[[2L]]
trend <- formula[c(1L,3L)]
## FIT #######################################
thecall <- call("lppm", X=Yexpr, trend=trend,
data=data, interaction=interaction)
ncall <- length(thecall)
argh <- list(...)
nargh <- length(argh)
if(nargh > 0) {
thecall[ncall + 1:nargh] <- argh
names(thecall)[ncall + 1:nargh] <- names(argh)
}
callenv <- list2env(as.list(data), parent=parent.frame())
result <- eval(thecall, envir=callenv)
result$call <- cl
result$callstring <- callstring
return(result)
}
lppm.lpp <- function(X, ..., eps=NULL, nd=1000, random=FALSE) {
Xname <- short.deparse(substitute(X))
callstring <- paste(short.deparse(sys.call()), collapse = "")
cl <- match.call()
nama <- names(list(...))
resv <- c("method", "forcefit")
if(any(clash <- resv %in% nama))
warning(paste(ngettext(sum(clash), "Argument", "Arguments"),
commasep(sQuote(resv[clash])),
"must not be used"))
stopifnot(inherits(X, "lpp"))
Q <- linequad(X, eps=eps, nd=nd, random=random)
fit <- ppm(Q, ..., method="mpl", forcefit=TRUE)
if(!is.poisson.ppm(fit))
warning("Non-Poisson models currently use Euclidean distance")
out <- list(X=X, fit=fit, Xname=Xname, call=cl, callstring=callstring)
class(out) <- "lppm"
return(out)
}
is.lppm <- function(x) { inherits(x, "lppm") }
# undocumented
as.ppm.lppm <- function(object) { object$fit }
fitted.lppm <- function(object, ..., dataonly=FALSE, new.coef=NULL,
leaveoneout=FALSE) {
pfit <- object$fit
v <- fitted(pfit, dataonly=dataonly, new.coef=new.coef,
leaveoneout=leaveoneout)
return(v)
}
predict.lppm <- function(object, ...,
type="trend", locations=NULL,
new.coef=NULL) {
type <- pickoption("type", type,
c(trend="trend", cif="cif", lambda="cif"))
X <- object$X
fit <- object$fit
L <- as.linnet(X)
if(!is.null(locations)) {
#' locations given; return a vector/matrix of predicted values
if(is.lpp(locations)) locations <- as.ppp(locations)
values <- predict(fit, locations=locations, type=type, new.coef=new.coef)
return(values)
}
# locations not given; want a pixel image
# pixellate the lines
Llines <- as.psp(L)
linemask <- as.mask.psp(Llines, ...)
lineimage <- as.im(linemask)
# extract pixel centres
xx <- rasterx.mask(linemask)
yy <- rastery.mask(linemask)
mm <- linemask$m
xx <- as.vector(xx[mm])
yy <- as.vector(yy[mm])
pixelcentres <- ppp(xx, yy, window=as.rectangle(linemask), check=FALSE)
pixdf <- data.frame(xc=xx, yc=yy)
# project pixel centres onto lines
p2s <- project2segment(pixelcentres, Llines)
projloc <- as.data.frame(p2s$Xproj)
projmap <- as.data.frame(p2s[c("mapXY", "tp")])
projdata <- cbind(pixdf, projloc, projmap)
# predict at the projected points
if(!is.multitype(fit)) {
values <- predict(fit, locations=projloc, type=type, new.coef=new.coef)
# map to nearest pixels
Z <- lineimage
Z[pixelcentres] <- values
# attach exact line position data
df <- cbind(projdata, values)
out <- linim(L, Z, df=df, restrict=FALSE)
} else {
# predict for each type
lev <- levels(marks(data.ppm(fit)))
out <- list()
for(k in seq(length(lev))) {
markk <- factor(lev[k], levels=lev)
locnk <- cbind(projloc, data.frame(marks=markk))
values <- predict(fit, locations=locnk, type=type, new.coef=new.coef)
Z <- lineimage
Z[pixelcentres] <- values
df <- cbind(projdata, values)
out[[k]] <- linim(L, Z, df=df, restrict=FALSE)
}
out <- as.solist(out)
names(out) <- as.character(lev)
}
return(out)
}
coef.lppm <- function(object, ...) {
coef(object$fit)
}
print.lppm <- function(x, ...) {
splat("Point process model on linear network")
print(x$fit)
terselevel <- spatstat.options('terse')
if(waxlyrical('extras', terselevel))
splat("Original data:", x$Xname)
if(waxlyrical('gory', terselevel))
print(as.linnet(x))
return(invisible(NULL))
}
summary.lppm <- function(object, ...) {
splat("Point process model on linear network")
print(summary(object$fit))
terselevel <- spatstat.options('terse')
if(waxlyrical('extras', terselevel))
splat("Original data:", object$Xname)
if(waxlyrical('gory', terselevel))
print(summary(as.linnet(object)))
return(invisible(NULL))
}
plot.lppm <- function(x, ..., type="trend") {
xname <- short.deparse(substitute(x))
y <- predict(x, type=type)
do.call(plot, resolve.defaults(list(y),
list(...),
list(main=xname)))
}
anova.lppm <- function(object, ..., test=NULL) {
stuff <- list(object=object, ...)
if(!is.na(hit <- match("override", names(stuff)))) {
warning("Argument 'override' is outdated and was ignored")
stuff <- stuff[-hit]
}
#' extract ppm objects where appropriate
mod <- sapply(stuff, is.lppm)
stuff[mod] <- lapply(stuff[mod], getElement, name="fit")
#' analysis of deviance or adjusted composite deviance
do.call(anova.ppm, append(stuff, list(test=test)))
}
update.lppm <- function(object, ...) {
stopifnot(inherits(object, "lppm"))
X <- object$X
fit <- object$fit
Xname <- object$Xname
callframe <- environment(formula(fit))
aargh <- list(...)
islpp <- sapply(aargh, is.lpp)
if(any(islpp)) {
# trap point pattern argument & convert to quadscheme
ii <- which(islpp)
if((npp <- length(ii)) > 1)
stop(paste("Arguments not understood:", npp, "lpp objects given"))
X <- aargh[[ii]]
aargh[[ii]] <- linequad(X)
}
isfmla <- sapply(aargh, inherits, what="formula")
if(any(isfmla)) {
# trap formula pattern argument, update it, evaluate LHS if required
jj <- which(isfmla)
if((nf <- length(jj)) > 1)
stop(paste("Arguments not understood:", nf, "formulae given"))
fmla <- aargh[[jj]]
fmla <- update(formula(object), fmla)
if(!is.null(lhs <- lhs.of.formula(fmla))) {
X <- eval(lhs, envir=list2env(list("."=X), parent=callframe))
Qpos <- if(any(islpp)) ii else (length(aargh) + 1L)
aargh[[Qpos]] <- linequad(X)
}
aargh[[jj]] <- rhs.of.formula(fmla)
}
newfit <- do.call(update.ppm,
append(list(fit), aargh),
envir=callframe)
if(!is.poisson.ppm(newfit))
warning("Non-Poisson models currently use Euclidean distance")
out <- list(X=X, fit=newfit, Xname=Xname)
class(out) <- "lppm"
return(out)
}
terms.lppm <- function(x, ...) {
terms(x$fit, ...)
}
logLik.lppm <- function(object, ...) {
logLik(object$fit, ...)
}
deviance.lppm <- function(object, ...) {
as.numeric(-2 * logLik(object, ...))
}
pseudoR2.lppm <- function(object, ..., keepoffset=TRUE) {
dres <- deviance(object, ..., warn=FALSE)
nullfmla <- . ~ 1
if(keepoffset && has.offset.term(object)) {
off <- attr(model.depends(object), "offset")
offterms <- row.names(off)[apply(off, 1, any)]
if(length(offterms)) {
nullrhs <- paste(offterms, collapse=" + ")
nullfmla <- as.formula(paste(". ~ ", nullrhs))
}
}
nullmod <- update(object, nullfmla, forcefit=TRUE)
dnul <- deviance(nullmod, warn=FALSE)
return(1 - dres/dnul)
}
formula.lppm <- function(x, ...) {
formula(x$fit, ...)
}
extractAIC.lppm <- function(fit, ...) {
extractAIC(fit$fit, ...)
}
as.owin.lppm <- function(W, ..., fatal=TRUE) {
stopifnot(inherits(W, "lppm"))
as.owin(as.linnet(W), ..., fatal=fatal)
}
Window.lppm <- function(X, ...) { as.owin(X) }
data.lppm <- function(object) { object$X }
model.images.lppm <- local({
model.images.lppm <- function(object, L=as.linnet(object), ...) {
stopifnot(inherits(object, "lppm"))
stopifnot(inherits(L, "linnet"))
m <- model.images(object$fit, W=as.rectangle(L), ...)
if(length(m)) {
## restrict images to L
type <- if(is.hyperframe(m)) "hyperframe" else
if(is.imlist(m)) "imlist" else
if(is.list(m) && all(sapply(m, is.im))) "imlist" else
stop("Internal error: model.images not understood", call.=FALSE)
switch(type,
imlist = {
## list of images: convert to list of linims
ZL <- netmask(L, template=m[[1L]])
m <- tolinims(m, L=L, imL=ZL)
},
hyperframe = {
## hyperframe, each column being a list of images
## extract columns
rownam <- row.names(m)
m <- as.list(m)
ZL <- netmask(L, template=m[[1L]][[1L]])
mm <- lapply(m, tolinims, L=L, imL=ZL)
m <- do.call(hyperframe, mm)
row.names(m) <- rownam
})
}
return(m)
}
netmask <- function(L, template) {
as.im(as.mask.psp(as.psp(L), xy=as.mask(template)))
}
tolinim <- function(x, L, imL) linim(L, eval.im(x * imL), restrict=FALSE)
tolinims <- function(x, L, imL) solapply(x, tolinim, L=L, imL=imL)
model.images.lppm
})
model.matrix.lppm <- function(object,
data=model.frame(object, na.action=NULL),
..., keepNA=TRUE) {
stopifnot(is.lppm(object))
if(missing(data)) data <- NULL
model.matrix(object$fit, data=data, ..., keepNA=keepNA)
}
model.frame.lppm <- function(formula, ...) {
stopifnot(inherits(formula, "lppm"))
model.frame(formula$fit, ...)
}
domain.lppm <- as.linnet.lppm <- function(X, ...) {
as.linnet(X$X, ...)
}
nobs.lppm <- function(object, ...) {
npoints(object$X)
}
is.poisson.lppm <- function(x) { is.poisson(x$fit) }
is.stationary.lppm <- function(x) { is.stationary(x$fit) }
is.multitype.lppm <- function(X, ...) { is.multitype(X$fit) }
is.marked.lppm <- function(X, ...) { is.marked(X$fit) }
vcov.lppm <- function(object, ...) {
if(!is.poisson(object))
stop("vcov.lppm is only implemented for Poisson models")
vcov(object$fit, ...)
}
valid.lppm <- function(object, ...) {
valid(object$fit, ...)
}
emend.lppm <- function(object, ...) {
object$fit <- emend(object$fit, ...)
return(object)
}