#
# summary.ppm.R
#
# summary() method for class "ppm"
#
# $Revision: 1.10 $ $Date: 2005/04/29 20:31:04 $
#
# summary.ppm()
# print.summary.ppm()
#
summary.ppm <- function(object, ..., quick=FALSE) {
verifyclass(object, "ppm")
x <- object
y <- list()
####### Check version #########################
ver <- object$version
antiquated <- is.null(ver) || !is.list(ver) ||
(ver$major == 1 && ver$minor < 5)
y$antiquated <- antiquated
####### Extract main data components #########################
QUAD <- object$Q
DATA <- QUAD$data
TREND <- x$trend
INTERACT <- x$interaction
if(is.null(INTERACT)) INTERACT <- Poisson()
####### Determine type of model ############################
y$no.trend <- identical.formulae(TREND, NULL) || identical.formulae(TREND, ~1)
y$stationary <- y$no.trend || identical.formulae(TREND, ~marks)
y$poisson <- is.null(INTERACT$family)
y$marked <- is.marked.ppp(DATA)
y$multitype <- y$marked && is.factor(DATA$marks)
if(y$marked) y$entries <- list(marks = DATA$marks)
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="")
if(is.logical(quick) && quick) {
class(y) <- "summary.ppm"
return(y)
}
###### Does it have external covariates? ####################
if(!antiquated) {
hc <- !is.null(x$covariates)
} 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
hc <- !is.null(callargs) && !is.na(pmatch("data", callargs))
}
y$has.covars <- hc
###### 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)
y$quad <- summary.quad(QUAD)
if(is.character(quick) && (quick == "no prediction"))
return(y)
###### Extract fitted model coefficients #########################
if(exists("is.R") && is.R())
theta <- x$coef # result of coef(glm(...))
else
theta <- x$theta # result of dummy.coef(glm(....))
y$entries$theta <- theta
# corresponding internal names of regressor variables
Vnames <- x$internal$Vnames
###### 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) {
lambda <- exp(theta[[1]])
if(!y$marked) {
y$trend$label <- "Uniform intensity"
y$trend$value <- lambda
} else {
y$trend$label <- "Uniform intensity for each mark level"
y$trend$value <- lambda
}
} else # process is at least one of: marked, nonstationary, non-poisson
if(y$stationary) {
if(!y$marked) {
# stationary non-poisson non-marked
y$trend$label <- "First order term"
y$trend$value <- c(beta=exp(theta[[1]]))
} else {
# stationary, marked
mrk <- DATA$marks
y$trend$label <-
if(y$poisson) "Fitted intensities" else "Fitted 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
}
} 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 <- theta
else {
agree <- outer(names(theta), Vnames, "==")
whichbits <- apply(!agree, 1, all)
trendbits <- theta[whichbits]
}
# decide whether there are 'labels within labels'
unlabelled <- unlist(lapply(trendbits,
function(x) { is.null(names(x)) } ))
if(all(unlabelled))
y$trend$value <- unlist(trendbits)
else {
y$trend$value <- list()
for(i in seq(trendbits))
y$trend$value[[i]] <-
if(unlabelled[i])
unlist(trendbits[i])
else
trendbits[[i]]
}
}
###### Interaction component #########################
if(!y$poisson) {
if(!is.null(INTERACT$interpret)) {
# invoke auto-interpretation feature
sensible <- (INTERACT$interpret)(x$coef, INTERACT)
header <- paste("Fitted", sensible$inames)
printable <- sensible$printable
} else {
# fallback
sensible <- NULL
header <- "Fitted interaction terms"
printable <- exp(unlist(theta[Vnames]))
}
y$interaction <- list(sensible=sensible,
header=header,
printable=printable)
}
class(y) <- "summary.ppm"
return(y)
}
print.summary.ppm <- function(x, ...) {
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")
cat("fitted by maximum pseudolikelihood\n")
cat("Call:\n")
print(x$args$call)
cat(paste("Edge correction: \'", x$args$correction, "\'\n", sep=""))
if(x$args$rbord > 0)
cat(paste("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
# names of interaction variables if any
Vnames <- x$Vnames
# their fitted coefficients
theta <- x$theta
# ----------- 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$has.covars)
cat("Model involves external covariates\n")
}
cat(paste("\n", x$trend$label, ":\n", sep=""))
tv <- x$trend$value
if(!is.list(tv))
print(tv)
else
for(i in seq(tv))
print(tv[[i]])
# ---- Interaction ----------------------------
if(!poisson) {
cat("\n\n ---- Interaction: -----\n\n")
print(x$entries$interaction)
cat(paste(x$interaction$header, ":\n", sep=""))
print(x$interaction$printable)
}
####### Gory details ###################################
cat("\n\n----------- gory details -----\n")
theta <- x$entries$theta
cat("\nFitted regular parameters (theta): \n")
print(theta)
cat("\nFitted exp(theta): \n")
print(exp(unlist(theta)))
return(invisible(NULL))
}
no.trend.ppm <- function(x) {
summary.ppm(x, quick=TRUE)$no.trend
}
is.stationary.ppm <- function(x) {
summary.ppm(x, quick=TRUE)$stationary
}
is.poisson.ppm <- function(x) {
summary.ppm(x, quick=TRUE)$poisson
}
is.marked.ppm <- function(X, ...) {
summary.ppm(X, quick=TRUE)$marked
}