clftest.R
#
# clftest.R
#
# $Revision: 1.9 $ $Date: 2012/05/15 09:07:53 $
#
# Monte Carlo tests for CSR (etc)
#
clf.test <- function(X, ..., rinterval=NULL, use.theo=FALSE) {
Xname <- deparse(substitute(X))
envelopeTest(X, ..., power=2,
use.theo=use.theo, rinterval=rinterval, Xname=Xname)
}
mad.test <- function(X, ..., rinterval=NULL, use.theo=FALSE) {
Xname <- deparse(substitute(X))
envelopeTest(X, ..., power=Inf,
use.theo=use.theo, rinterval=rinterval, Xname=Xname)
}
envelopeTest <- function(X, ...,
power=1, rinterval=NULL,
use.theo=FALSE, Xname=NULL, verbose=TRUE,
internal=NULL) {
if(is.null(Xname)) Xname <- deparse(substitute(X))
check.1.real(power)
explain.ifnot(power >= 0)
# compute or extract simulated functions
if(use.theo) {
# using theoretical function as reference.
# ensure envelope object includes theoretical function.
internal <- resolve.defaults(internal, list(csr=TRUE))
}
X <- envelope(X, ..., savefuns=TRUE, Yname=Xname,
internal=internal, verbose=verbose)
Y <- attr(X, "simfuns")
# extract values
r <- with(X, .x)
obs <- with(X, .y)
sim <- as.matrix(as.data.frame(Y))[, -1]
nsim <- ncol(sim)
# choose function as reference
has.theo <- ("theo" %in% names(X))
if(use.theo && !has.theo)
warning("No theoretical function available; use.theo ignored")
if(use.theo && has.theo) {
reference <- with(X, theo)
used.theo <- TRUE
} else {
# compute sample mean of simulations *and* observed
reference <- apply(cbind(sim, obs), 1, mean, na.rm=TRUE)
used.theo <- FALSE
}
# determine interval of r values for computation
if(!is.null(rinterval)) {
stopifnot(is.numeric(rinterval))
stopifnot(length(rinterval) == 2)
stopifnot(rinterval[1] < rinterval[2])
if(max(r) < rinterval[2]) {
oldrinterval <- rinterval
rinterval <- intersect.ranges(rinterval, range(r))
if(verbose)
warning(paste("The interval", prange(oldrinterval),
"is too large for the available data;",
"it has been trimmed to", prange(rinterval)))
}
ok <- (rinterval[1] <= r & r <= rinterval[2])
obs <- obs[ok]
sim <- sim[ok, ]
reference <- reference[ok]
} else rinterval <- range(r)
# compute test statistic
if(is.infinite(power)) {
# MAD
devdata <- max(abs(obs-reference))
names(devdata) <- "mad"
devsim <- apply(abs(sim-reference), 2, max)
testname <- "Maximum absolute deviation test"
} else {
a <- diff(rinterval) * (if(used.theo) 1 else (nsim/(nsim - 1))^power)
if(power == 2) {
# Cramer-von Mises
devdata <- a * mean((obs - reference)^2)
names(devdata) <- "u"
devsim <- a * colMeans((sim - reference)^2)
testname <- "Cressie-Loosmore-Ford test"
} else if(power == 1) {
# integral absolute deviation
devdata <- a * mean(abs(obs - reference))
names(devdata) <- "L1"
devsim <- a * colMeans(abs(sim - reference))
testname <- "Integral absolute deviation test"
} else {
# general p
devdata <- a * mean((abs(obs - reference)^power))
names(devdata) <- "Lp"
devsim <- a * colMeans((abs(sim - reference)^power))
testname <- paste("Integrated", ordinal(power), "Power Deviation test")
}
}
# compute rank
datarank <- sum(devdata < devsim) + sum(devdata == devsim)/2 + 1
pvalue <- datarank/(nsim+1)
# bookkeeping
statistic <- data.frame(devdata, rank=datarank)
colnames(statistic)[1] <- names(devdata)
e <- attr(X, "einfo")
nullmodel <-
if(identical(e$csr, TRUE)) "CSR" else
if(!is.null(e$simtype)) {
switch(e$simtype,
csr = "CSR",
rmh = paste("fitted",
if(identical(e$pois, TRUE)) "Poisson" else "Gibbs",
"model"),
kppm = "fitted cluster model",
expr = "model simulated by evaluating expression",
list = "model simulated by drawing patterns from a list",
"unrecognised model")
} else "unrecognised model"
fname <- deparse(attr(X, "ylab"))
testname <- c(paste(testname, "of", nullmodel),
paste("Monte Carlo test based on", nsim, "simulations"),
paste("Summary function:", fname),
paste("Reference function:",
if(used.theo) "theoretical" else "sample mean"))
result <- structure(list(statistic = statistic,
p.value = pvalue,
method = testname,
data.name = e$Yname),
class = "htest")
return(result)
}