#
# summary.ppm.R
#
# summary() method for class "ppm"
#
# $Revision: 1.51 $ $Date: 2012/01/31 07:52:12 $
#
# summary.ppm()
# print.summary.ppm()
#
summary.ppm <- function(object, ..., quick=FALSE) {
verifyclass(object, "ppm")
x <- object
y <- list()
####### Extract main data components #########################
QUAD <- object$Q
DATA <- QUAD$data
TREND <- x$trend
INTERACT <- x$interaction
if(is.null(INTERACT)) INTERACT <- Poisson()
####### Check version #########################
mpl.ver <- versionstring.ppm(object)
int.ver <- versionstring.interact(INTERACT)
current <- versionstring.spatstat()
virgin <- min(package_version(c(mpl.ver, int.ver)))
y$antiquated <- antiquated <- (virgin <= package_version("1.5"))
y$old <- old <- (virgin < majorminorversion(current))
y$version <- as.character(virgin)
####### Determine type of model ############################
y$entries <- list()
y$no.trend <- identical.formulae(TREND, NULL) || identical.formulae(TREND, ~1)
y$trendvar <- trendvar <- variablesinformula(TREND)
y$stationary <- y$no.trend || all(trendvar == "marks")
y$poisson <- is.poisson.interact(INTERACT)
y$marked <- is.marked.ppp(DATA)
y$multitype <- is.multitype.ppp(DATA)
y$marktype <- if(y$multitype) "multitype" else
if(y$marked) "marked" else "unmarked"
if(y$marked) y$entries$marks <- marks(DATA)
y$name <- paste(
if(y$stationary) "Stationary " else "Nonstationary ",
if(y$poisson) {
if(y$multitype) "multitype "
else if(y$marked) "marked "
else ""
},
INTERACT$name,
sep="")
###### Fitting algorithm ########################################
y$method <- x$method
y$problems <- x$problems
y$fitter <- x$fitter
if(y$fitter %in% c("glm", "gam"))
y$converged <- x$internal$glmfit$converged
###### Coefficients were changed after fit? #####################
y$projected <- identical(x$projected, TRUE)
y$changedcoef <- y$projected || (y$method != "mpl")
###### Extract fitted model coefficients #########################
y$entries$coef <- COEFS <- x$coef
y$entries$Vnames <- Vnames <- x$internal$Vnames
y$entries$IsOffset <- x$internal$IsOffset
###### Extract fitted interaction and summarise #################
FITIN <- fitin(x)
y$interaction <- summary(FITIN)
# Exit here if quick=TRUE
if(identical(quick, TRUE)) {
class(y) <- "summary.ppm"
return(y)
}
###### Does it have external covariates? ####################
if(!antiquated) {
covars <- x$covariates
y$has.covars <- !is.null(covars) && (length(covars) > 0)
y$covnames <- names(covars)
used <- (y$trendvar %in% names(covars))
y$covars.used <- y$trendvar[used]
y$uses.covars <- any(used)
# describe covariates
covtype <- function(x) {
if(is.im(x)) "im" else
if(inherits(x, "distfun")) "distfun" else
if(is.function(x)) "function" else
if(is.owin(x)) "owin" else
if(is.numeric(x) && length(x) == 1) "number" else "unknown"
}
y$covar.type <- ctype <- unlist(lapply(covars, covtype))
y$covar.descrip <- y$covar.type
# are there any functions?
y$has.funcs <- any(isfun <- (ctype == "function"))
# do covariates depend on additional arguments?
if(y$has.funcs) {
y$covfunargs <- x$covfunargs
funs <- covars[isfun]
fdescrip <- function(f) {
alist <- paste(names(formals(f)), collapse=", ")
paste("function(", alist, ")", sep="")
}
y$covar.descrip[isfun] <- unlist(lapply(funs, fdescrip))
# find any extra arguments (after args 1 & 2) explicitly named
xargs <- function(f) {
ar <- names(formals(f))[-(1:2)]
return(ar[ar != "..."])
}
fargs <- lapply(funs, xargs)
nxargs <- unlist(lapply(fargs, length))
y$has.xargs <- any(nxargs > 0)
if(y$has.xargs) {
# identify which function arguments are fixed in the call
fmap <- data.frame(Covariate=rep(names(funs), nxargs),
Argument=unlist(fargs))
fmap$Given <- (fmap$Argument %in% names(y$covfunargs))
y$xargmap <- fmap
}
}
} else {
# Antiquated format
# Interpret the function call instead
callexpr <- parse(text=x$call)
callargs <- names(as.list(callexpr[[1]]))
# Data frame of covariates was called 'data' in versions up to 1.4-x
y$has.covars <- !is.null(callargs) && !is.na(pmatch("data", callargs))
# conservative guess
y$uses.covars <- y$has.covars
y$covfunargs <- NULL
}
###### Arguments in call ####################################
y$args <- x[c("call", "correction", "rbord")]
####### Main data components #########################
y$entries <- append(list(quad=QUAD,
data=DATA,
interaction=INTERACT),
y$entries)
####### Summarise data ############################
y$data <- summary.ppp(DATA, checkdup=FALSE)
y$quad <- summary.quad(QUAD, checkdup=FALSE)
if(is.character(quick) && (quick == "no prediction"))
return(y)
###### Trend component #########################
y$trend <- list()
y$trend$name <- if(y$poisson) "Intensity" else "Trend"
y$trend$formula <- if(y$no.trend) NULL else TREND
if(y$poisson && y$no.trend) {
# uniform Poisson process
y$trend$value <- lambda <- exp(COEFS[[1]])
y$trend$label <- switch(y$marktype,
unmarked="Uniform intensity",
multitype="Uniform intensity for each mark level",
marked="Uniform intensity in product space",
"")
} else if(y$stationary) {
# stationary
switch(y$marktype,
unmarked={
# stationary non-poisson non-marked
y$trend$label <- "First order term"
y$trend$value <- c(beta=exp(COEFS[[1]]))
},
multitype={
# stationary, multitype
mrk <- marks(DATA)
y$trend$label <-
if(y$poisson) "Intensities" else "First order terms"
# Use predict.ppm to evaluate the fitted intensities
lev <- factor(levels(mrk), levels=levels(mrk))
nlev <- length(lev)
marx <- list(x=rep(0, nlev), y=rep(0, nlev), marks=lev)
betas <- predict(x, locations=marx, type="trend")
names(betas) <- paste("beta_", as.character(lev), sep="")
y$trend$value <- betas
},
marked={
# stationary, marked
y$trend$label <- "Fitted intensity coefficients"
y$trend$value <- blankcoefnames(COEFS)
})
} else {
# not stationary
y$trend$label <- "Fitted coefficients for trend formula"
# extract trend terms without trying to understand them much
if(is.null(Vnames))
trendbits <- COEFS
else {
agree <- outer(names(COEFS), Vnames, "==")
whichbits <- apply(!agree, 1, all)
trendbits <- COEFS[whichbits]
}
y$trend$value <- blankcoefnames(trendbits)
}
# ----- parameters with SE --------------------------
if(y$poisson && (length(COEFS) > 0)) {
# compute standard errors
se <- x$internal$se
if(is.null(se)) {
vc <- vcov(x, matrix.action="silent")
if(!is.null(vc))
se <- sqrt(diag(vc))
}
if(!is.null(se)) {
two <- qnorm(0.975)
lo <- COEFS - two * se
hi <- COEFS + two * se
pval <- 2 * pnorm(abs(COEFS)/se, lower.tail=FALSE)
psig <- cut(pval, c(0,0.001, 0.01, 0.05, 1),
labels=c("***", "**", "*", " "))
notapplic <- names(COEFS) %in% c("(Intercept)", "log(lambda)")
psig[notapplic] <- " "
# table of coefficient estimates with SE and 95% CI
y$coefs.SE.CI <- data.frame(Estimate=COEFS, S.E.=se, Ztest=psig,
CI95.lo=lo, CI95.hi=hi)
}
}
class(y) <- "summary.ppm"
return(y)
}
coef.summary.ppm <- function(object, ...) {
object$coefs.SE.CI
}
print.summary.ppm <- function(x, ...) {
if(x$old)
warning("Model was fitted by an older version of spatstat")
if(is.null(x$args)) {
# this is the quick version
cat(paste(x$name, "\n"))
return(invisible(NULL))
}
# otherwise - full details
cat("Point process model\n")
methodchosen <-
if(is.null(x$method))
"unspecified method"
else
switch(x$method,
mpl={
if(!x$poisson) {
"maximum pseudolikelihood (Berman-Turner approximation)"
} else {
# Poisson process
if(x$fitter == "exact") "maximum likelihood" else
"maximum likelihood (Berman-Turner approximation)"
}
},
ho="Huang-Ogata method (approximate maximum likelihood)",
paste("unrecognised method", sQuote(x$method)))
cat(paste("Fitting method:", methodchosen, "\n"))
howfitted <- switch(x$fitter,
exact= "analytically",
gam = "using gam()",
glm = "using glm()",
ho = NULL,
paste("using unrecognised fitter", sQuote(x$fitter)))
if(!is.null(howfitted)) cat(paste("Model was fitted", howfitted, "\n"))
if(x$fitter %in% c("glm", "gam")) {
if(x$converged) cat("Algorithm converged\n")
else cat("*** Algorithm did not converge ***\n")
}
if(x$projected)
cat("Fit was projected to obtain a valid point process model\n")
cat("Call:\n")
print(x$args$call)
if(x$old)
cat(paste("** Executed by old spatstat version", x$version, " **\n"))
cat(paste("Edge correction:", dQuote(x$args$correction), "\n"))
if(x$args$correction == "border")
cat(paste("\t[border correction distance r =", x$args$rbord,"]\n"))
cat("\n----------------------------------------------------\n")
# print summary of quadrature scheme
print(x$quad)
cat("\n----------------------------------------------------\n")
cat("FITTED MODEL:\n\n")
# This bit is currently identical to print.ppm()
# except for a bit more fanfare
# and the inclusion of the 'gory details' bit
notrend <- x$no.trend
stationary <- x$stationary
poisson <- x$poisson
markeddata <- x$marked
multitype <- x$multitype
markedpoisson <- poisson && markeddata
# ----------- Print model type -------------------
cat(x$name)
cat("\n")
if(markeddata) mrk <- x$entries$marks
if(multitype) {
cat("Possible marks: \n")
cat(paste(levels(mrk)))
}
# ----- trend --------------------------
cat(paste("\n\n ---- ", x$trend$name, ": ----\n\n", sep=""))
if(!notrend) {
cat("Trend formula: ")
print(x$trend$formula)
if(x$uses.covars)
cat(paste("Model depends on external",
ngettext(length(x$covars.used), "covariate", "covariates"),
commasep(sQuote(x$covars.used)), "\n"))
}
if(x$has.covars) {
if(notrend || !x$uses.covars)
cat("Model object contains external covariates\n")
if(!is.null(cd <- x$covar.descrip)) {
# print description of each covariate
cat("\nCovariates provided:\n")
namescd <- names(cd)
for(i in seq_along(cd))
cat(paste("\t", namescd[i], ": ", cd[i], "\n", sep=""))
}
if(!is.null(cfa <- x$covfunargs) && length(cfa) > 0) {
cat("Covariate function arguments (covfunargs) provided:\n")
namescfa <- names(cfa)
for(i in seq_along(cfa)) {
cat(paste(namescfa[i], "= "))
cfai <- cfa[[i]]
if(is.numeric(cfai) && length(cfai) == 1) {
cat(paste(cfai, "\n"))
} else print(cfa[[i]])
}
}
}
cat(paste("\n", x$trend$label, ":\n", sep=""))
tv <- x$trend$value
if(!is.list(tv))
print(tv)
else
for(i in seq_along(tv))
print(tv[[i]])
# table of coefficient estimates with SE and 95% CI
if(!is.null(cose <- x$coefs.SE.CI)) {
cat("\n")
print(cose)
}
# ---- Interaction ----------------------------
if(!poisson) {
cat("\n\n ---- Interaction: -----\n\n")
print(x$interaction)
}
####### Gory details ###################################
cat("\n\n----------- gory details -----\n")
COEFS <- x$entries$coef
cat("\nFitted regular parameters (theta): \n")
print(COEFS)
cat("\nFitted exp(theta): \n")
print(exp(unlist(COEFS)))
##### Warnings issued #######
probs <- x$problems
if(!is.null(probs) && is.list(probs) && (length(probs) > 0))
lapply(probs,
function(a) {
if(is.list(a) && !is.null(p <- a$print))
cat(paste("Problem:\n", p, "\n\n"))
})
return(invisible(NULL))
}
no.trend.ppm <- function(x) {
summary.ppm(x, quick=TRUE)$no.trend
}
is.stationary <- function(x) {
UseMethod("is.stationary")
}
is.poisson <- function(x) {
UseMethod("is.poisson")
}
is.stationary.ppm <- function(x) {
TREND <- x$trend
if(is.null(TREND) || identical.formulae(TREND, ~1))
return(TRUE)
if(all(variablesinformula(TREND) == "marks"))
return(TRUE)
return(FALSE)
}
is.poisson.ppm <- function(x) {
stopifnot(is.ppm(x))
y <- x$interaction
if(is.null(y)) y <- Poisson()
is.poisson.interact(y)
}
is.marked.ppm <- function(X, ...) {
summary.ppm(X, quick=TRUE)$marked
}
is.multitype.ppm <- function(X, ...) {
summary.ppm(X, quick=TRUE)$multitype
}
blankcoefnames <- function(x) {
# remove name labels from ppm coefficients
# First decide whether there are 'labels within labels'
unlabelled <- unlist(lapply(x,
function(z) { is.null(names(z)) } ))
if(all(unlabelled))
value <- unlist(x)
else {
value <- list()
for(i in seq_along(x))
value[[i]] <- if(unlabelled[i]) unlist(x[i]) else x[[i]]
}
return(value)
}