Fest.R
#
# Fest.R
#
# Computes estimates of the empty space function
#
# $Revision: 4.43 $ $Date: 2016/10/04 01:04:13 $
#
Fhazard <- function(X, ...) {
Z <- Fest(X, ...)
if(!any(names(Z) == "km"))
stop("Kaplan-Meier estimator 'km' is required for hazard rate")
## strip off Poisson F
Z <- Z[, (colnames(Z) != "theo")]
## relabel the fv object
Z <- rebadge.fv(Z,
new.ylab=quote(h(r)),
new.fname="h",
tags=c("hazard", "theohaz"),
new.tags=c("hazard", "theo"),
new.labl=c("hat(%s)[km](r)", "%s[pois](r)"),
new.desc=c(
"Kaplan-Meier estimate of %s",
"theoretical Poisson %s"),
new.dotnames=c("hazard", "theo"),
new.preferred="hazard")
## strip off unwanted bits
Z <- Z[, c("r", "hazard", "theo")]
return(Z)
}
Fest <- function(X, ..., eps = NULL, r=NULL, breaks=NULL,
correction=c("rs", "km", "cs"),
domain=NULL) {
verifyclass(X, "ppp")
if(!is.null(domain))
stopifnot(is.subset.owin(domain, Window(X)))
rorbgiven <- !is.null(r) || !is.null(breaks)
## Intensity estimate
W <- X$window
npts <- npoints(X)
lambda <- npts/area(W)
## First discretise
dwin <- as.mask(W, eps=eps)
dX <- ppp(X$x, X$y, window=dwin, check=FALSE)
## histogram breakpoints
rmaxdefault <- rmax.rule("F", dwin, lambda)
breaks <- handle.r.b.args(r, breaks, dwin, eps,
rmaxdefault=rmaxdefault)
rvals <- breaks$r
rmax <- breaks$max
if(rorbgiven) check.finespacing(rvals,
if(is.null(eps)) NULL else eps/4,
dwin,
rmaxdefault=rmaxdefault,
action="fatal",
rname="r",
context="in Fest(X, r)")
## choose correction(s)
# correction.given <- !missing(correction) && !is.null(correction)
if(is.null(correction)) {
correction <- c("rs", "km", "cs")
} else correction <- pickoption("correction", correction,
c(none="none",
border="rs",
rs="rs",
KM="km",
km="km",
Kaplan="km",
cs="cs",
ChiuStoyan="cs",
Hanisch="cs",
han="cs",
best="km"),
multi=TRUE)
## initialise fv object
df <- data.frame(r=rvals, theo=1-exp(-lambda * pi * rvals^2))
Z <- fv(df, "r", substitute(F(r), NULL), "theo", . ~ r,
c(0,rmax),
c("r", "%s[pois](r)"),
c("distance argument r", "theoretical Poisson %s"),
fname="F")
nr <- length(rvals)
zeroes <- numeric(nr)
## compute distances and censoring distances
if(X$window$type == "rectangle") {
## original data were in a rectangle
## output of exactdt() is sufficient
e <- exactdt(dX)
dist <- e$d
bdry <- e$b
if(!is.null(domain)) {
ok <- inside.owin(raster.xy(e$w), , domain)
dist <- dist[ok]
bdry <- bdry[ok]
}
} else {
## window is irregular..
# Distance transform & boundary distance for all pixels
e <- exactdt(dX)
b <- bdist.pixels(dX$window, style="matrix")
## select only those pixels inside mask
mm <- dwin$m
if(!is.null(domain)) {
ok <- inside.owin(raster.xy(e$w), , domain)
mm <- as.vector(mm) & ok
}
dist <- e$d[mm]
bdry <- b[mm]
}
## censoring indicators
d <- (dist <= bdry)
## observed distances
o <- pmin.int(dist, bdry)
## start calculating estimates of F
if("none" %in% correction) {
## UNCORRECTED e.d.f. of empty space distances
if(npts == 0)
edf <- zeroes
else {
hh <- hist(dist[dist <= rmax],breaks=breaks$val,plot=FALSE)$counts
edf <- cumsum(hh)/length(dist)
}
Z <- bind.fv(Z, data.frame(raw=edf), "hat(%s)[raw](r)",
"uncorrected estimate of %s", "raw")
}
if("cs" %in% correction) {
## Chiu-Stoyan correction
if(npts == 0)
cs <- zeroes
else {
## uncensored distances
x <- dist[d]
## weights
a <- eroded.areas(W, rvals)
## calculate Hanisch estimator
h <- hist(x[x <= rmax], breaks=breaks$val, plot=FALSE)$counts
H <- cumsum(h/a)
cs <- H/max(H[is.finite(H)])
}
## add to fv object
Z <- bind.fv(Z, data.frame(cs=cs),
"hat(%s)[cs](r)",
"Chiu-Stoyan estimate of %s",
"cs")
}
if(any(correction %in% c("rs", "km"))) {
## calculate Kaplan-Meier and/or border corrected (Reduced Sample) estimators
want.rs <- "rs" %in% correction
want.km <- "km" %in% correction
selection <- c(want.rs, want.km, want.km, want.km)
tags <- c("rs", "km", "hazard", "theohaz")[selection]
labels <- c("hat(%s)[bord](r)", "hat(%s)[km](r)",
"hat(h)[km](r)", "h[pois](r)")[selection]
descr <- c("border corrected estimate of %s",
"Kaplan-Meier estimate of %s",
"Kaplan-Meier estimate of hazard function h(r)",
"theoretical Poisson hazard h(r)")[selection]
if(npts == 0) {
result <- as.data.frame(matrix(0, nr, length(tags)))
names(result) <- tags
} else {
result <- km.rs.opt(o, bdry, d, breaks, KM=want.km, RS=want.rs)
result$theohaz <- 2 * pi * lambda * rvals
result <- as.data.frame(result[tags])
}
## add to fv object
Z <- bind.fv(Z, result,
labels, descr, if(want.km) "km" else "rs")
}
## wrap up
unitname(Z) <- unitname(X)
## remove 'hazard' from the dotnames
nama <- names(Z)
fvnames(Z, ".") <- rev(setdiff(nama, c("r", "hazard", "theohaz")))
## determine recommended plot range
attr(Z, "alim") <- with(Z, range(.x[is.finite(.y) & .y <= 0.9]))
return(Z)
}