nncorr.R
#
# nncorr.R
#
# $Revision: 1.7 $ $Date: 2011/10/11 10:45:24 $
#
nnmean <- function(X) {
stopifnot(is.ppp(X) && is.marked(X))
m <- numeric.columns(marks(X), logical=TRUE, others="na")
nv <- ncol(m)
nnid <- nnwhich(X)
ok <- (nndist(X) <= bdist.points(X))
if(!any(ok))
stop("Insufficient data")
numer <- unlist(lapply(as.data.frame(m[nnid[ok], ]), mean, na.rm=TRUE))
denom <- unlist(lapply(as.data.frame(m), mean, na.rm=TRUE))
ans <- rbind(unnormalised=numer,
normalised =numer/denom)
return(ans)
}
nnvario <- function(X) {
stopifnot(is.ppp(X) && is.marked(X))
m <- numeric.columns(marks(X), logical=TRUE, others="na")
f <- function(m1,m2) { ((m1-m2)^2)/2 }
ans <- nncorr(X %mark% m, f, denominator=diag(var(m)))
return(ans)
}
nncorr <- function(X, f = function(m1,m2) { m1 * m2}, ...,
use = "all.obs",
method = c("pearson", "kendall", "spearman"),
denominator=NULL) {
stopifnot(is.ppp(X) && is.marked(X))
m <- as.data.frame(marks(X))
nv <- ncol(m)
if(nv == 1) colnames(m) <- ""
#
if(missing(method) || is.null(method))
method <- "pearson"
#
if(missing(f)) f <- NULL
if(!is.null(f) && !is.function(f)) {
if(nv == 1) stop("f should be a function")
# could be a list of functions
if(!(is.list(f) && all(unlist(lapply(f, is.function)))))
stop("f should be a function or a list of functions")
if(length(f) != nv)
stop("Length of list f does not match number of mark variables")
}
# optional denominator(s)
if(!is.null(denominator) && !(length(denominator) %in% c(1, nv)))
stop("Denominator has incorrect length")
# multi-dimensional case
if(nv > 1) {
# replicate things
if(is.function(f)) f <- rep(list(f), nv)
if(length(denominator) <= 1) denominator <- rep(list(denominator), nv)
#
result <- matrix(NA, nrow=3, ncol=nv)
outnames <- c("unnormalised", "normalised", "correlation")
dimnames(result) <- list(outnames, colnames(m))
for(j in 1:nv) {
mj <- m[,j, drop=FALSE]
denj <- denominator[[j]]
nncj <- nncorr(X %mark% mj, f=f[[j]], use=use, method=method,
denominator=denj)
kj <- length(nncj)
result[1:kj,j] <- nncj
}
if(all(is.na(result[3, ]))) result <- result[1:2, ]
return(result)
}
# one-dimensional
m <- m[,1,drop=TRUE]
# select 'f' appropriately for X
chk <- check.testfun(f, X=X)
f <- chk$f
ftype <- chk$ftype
# denominator
Efmm <-
if(!is.null(denominator)) denominator else
switch(ftype,
mul={
mean(m)^2
},
equ={
sum(table(m)^2)/length(m)^2
},
general={
mean(outer(m, m, f, ...))
})
# border method
nn <- nnwhich(X)
ok <- (nndist(X) <= bdist.points(X))
if(!any(ok))
stop("Insufficient data")
mY <- m[nn[ok]]
mX <- m[ok]
Efmk <- switch(ftype,
mul = {
mean(mX * mY, ...)
},
equ = {
mean(mX == mY, ...)
},
general = {
mean(f(mX, mY, ...))
})
#
answer <- c(unnormalised=Efmk,
normalised=Efmk/Efmm)
if(ftype == "mul") {
classic <- cor(mX, mY, use=use, method=method)
answer <- c(answer, correlation=classic)
}
return(answer)
}