Nothing
#'
#' Header for all (concatenated) test files
#'
#' Require spatstat.model
#' Obtain environment variable controlling tests.
#'
#' $Revision: 1.5 $ $Date: 2020/04/30 05:31:37 $
require(spatstat.model)
FULLTEST <- (nchar(Sys.getenv("SPATSTAT_TEST", unset="")) > 0)
ALWAYS <- TRUE
cat(paste("--------- Executing",
if(FULLTEST) "** ALL **" else "**RESTRICTED** subset of",
"test code -----------\n"))
#
# tests/ppmBadData.R
#
# $Revision: 1.6 $ $Date: 2020/04/30 05:23:52 $
# Testing robustness of ppm and support functions
# when data are rubbish
local({
if(ALWAYS) {
## from Rolf: very large proportion of data is NA
SEED <- 42
K <- 101
A <- 500
X <- seq(0, A, length=K)
G <- expand.grid(x=X, y=X)
FOO <- function(x,y) { sin(x)^2 + cos(y)^2 }
M1 <- im(matrix(FOO(G$x, G$y), K, K), xcol=X, yrow=X)
M <- im(matrix(FOO(G$x, G$y), K, K))
BAR <- function(x) { exp(-6.618913 + 5.855337 * x - 8.432483 * x^2) }
V <- im(BAR(M$v), xcol=X, yrow=X)
# V <- eval.im(exp(-6.618913 + 5.855337 * M - 8.432483 * M^2))
set.seed(SEED)
Y <- rpoispp(V)
fY <- ppm(Y ~cv + I(cv^2), data=list(cv=M), correction="translate")
diagnose.ppm(fY)
lurking(fY, covariate=as.im(function(x,y){x}, square(A)), type="raw")
}
if(ALWAYS) {
## from Andrew Bevan: numerical overflow, ill-conditioned Fisher information
SEED <- 42
nongranite<- owin(poly = list(x = c(0, 8500, 7000, 6400, 6400, 6700, 7000, 7200, 7300, 8000, 8100, 8800, 9500, 10000, 10000, 0), y = c(0, 0, 2000, 3800, 4000, 5000, 6500, 7400, 7500, 8000, 8100, 9000, 9500, 9600, 10000, 10000)))
## Trend on raster grid
rain <- as.im(X=function(x,y) { x^2 + y^2 }, W=nongranite, dimyx=100)
## Generate a point pattern via a Lennard-Jones process
set.seed(SEED)
mod4<- rmhmodel(cif="lennard",
par=list(beta=1, sigma=250, epsilon=2.2),
trend=rain, w=nongranite)
ljtr<- rmh(mod4, start=list(n.start=80), control=list(p=1, nrep=1e5))
## Fit a point process model to the pattern with rain as a covariate
## NOTE INCORRECT TREND FORMULA
ljtrmod <- ppm(ljtr, trend= ~ Z, interaction=NULL, data=list(Z=rain))
ss <- summary(ljtrmod)
}
if(FULLTEST) {
## From Ege
## Degenerate but non-null argument 'covariates'
xx <- list()
names(xx) <- character(0)
fit <- ppm(cells ~x, covariates = xx)
st <- summary(fit)
}
})
#' tests/ppmclass.R
#'
#' Class support for ppm
#'
#' $Revision: 1.9 $ $Date: 2022/03/07 03:26:09 $
if(FULLTEST) {
local({
#' (1) print.ppm, summary.ppm, print.summary.ppm
Z <- as.im(function(x,y){x}, Window(cells))
fitZ <- ppm(cells ~ Z)
print(fitZ)
print(summary(fitZ))
#' logistic
fitl <- ppm(swedishpines ~ x+y, method="logi")
print(fitl)
print(summary(fitl))
#' Model with covariate arguments
f <- function(x,y,b) { x+b }
fitf <- ppm(cells ~ f, covfunargs=list(b=1))
print(fitf)
print(summary(fitf))
#' Invalid model
fitN <- ppm(redwood ~ 1, Strauss(0.1))
print(fitN)
print(summary(fitN))
#' standard errors in output
fat <- ppm(cells ~ x, Strauss(0.12))
op <- spatstat.options(print.ppm.SE='always')
print(fat)
spatstat.options(print.ppm.SE='never')
print(fat)
print(fitZ)
spatstat.options(op)
## (2) plot.ppm
plot(fitZ)
plot(fat, trend=FALSE, cif=FALSE, se=FALSE)
## (3) emend.ppm
fitZe <- emend(fitZ, trace=TRUE)
ZZ <- Z
fitZZ <- ppm(cells ~ Z + ZZ)
fitZZe <- emend(fitZZ, trace=TRUE)
fitOK <- ppm(redwood ~1, Strauss(0.1), emend=TRUE)
print(fitOK)
fitNot <- ppm(redwood ~1, Strauss(0.1))
fitSlow <- emend(fitNot, trace=TRUE)
print(fitSlow)
op <- spatstat.options(project.fast=TRUE)
fitFast <- emend(fitNot, trace=TRUE)
print(fitFast)
fitZZe <- emend(fitZZ, trace=TRUE)
spatstat.options(op)
#' (4) methods for other generics
logLik(fitZ, absolute=TRUE)
unitname(fitZ)
unitname(fat) <- c("metre", "metres")
is.expandable(fitf)
fit0 <- update(fitZ, . ~ 1)
anova(fit0, fitZ, override=TRUE)
interactionfamilyname(fat)
interactionorder(fat)
hardcoredist(fat)
#' (5) miscellaneous
## example from Robert Aue - handling offsets
X <- demohyper$Points[[1]]
GH <- Hybrid(G=Geyer(r=0.1, sat=3), H=Hardcore(0.01))
fit <- ppm(X ~ 1, GH)
valid.ppm(fit)
#' hard core distance of hybrid
hardcoredist(fit)
#' interaction order of hybrid
interactionorder(fit)
#' case of boundingbox
boundingbox(cells, ppm(cells ~ 1))
})
reset.spatstat.options()
}
#
# tests/ppmgam.R
#
# Test ppm with use.gam=TRUE
#
# $Revision: 1.4 $ $Date: 2020/04/30 05:23:52 $
#
if(FULLTEST) {
local({
fit <- ppm(nztrees ~s(x,y), use.gam=TRUE)
mm <- model.matrix(fit)
mf <- model.frame(fit)
v <- vcov(fit)
prd <- predict(fit)
})
}
#'
#' tests/ppmlogi.R
#'
#' Tests of ppm(method='logi')
#' and related code (predict, leverage etc)
#'
#' $Revision: 1.15 $ $Date: 2020/04/30 05:23:52 $
#'
local({
if(FULLTEST) {
fit <- ppm(cells ~x, method="logi")
f <- fitted(fit)
p <- predict(fit)
u <- summary(fit)
fitS <- ppm(cells ~x, Strauss(0.12), method="logi")
fS <- fitted(fitS)
pS <- predict(fitS)
uS <- summary(fitS)
print(uS)
plot(leverage(fit))
plot(influence(fit))
plot(dfbetas(fit))
plot(leverage(fitS))
plot(influence(fitS))
plot(dfbetas(fitS))
}
if(FULLTEST) {
#' same with hard core - A1 is singular
fitH <- ppm(cells ~x, Strauss(0.08), method="logi")
print(fitH)
fH <- fitted(fitH)
pH <- predict(fitH)
uH <- summary(fitH)
print(uH)
plot(leverage(fitH))
plot(influence(fitH))
plot(dfbetas(fitH))
}
if(FULLTEST) {
#' logistic fit to data frame of covariates
z <- c(rep(TRUE, 5), rep(FALSE, 5))
df <- data.frame(A=z + 2* runif(10),
B=runif(10))
Y <- quadscheme.logi(runifpoint(5), runifpoint(5))
fut <- ppm(Y ~ A+B, data=df, method="logi")
sf <- summary(fut)
print(sf)
}
if(FULLTEST) {
#' vblogit code, just to check that it runs.
fee <- ppm(cells ~ x, method="VBlogi", nd=21)
print(fee)
summary(fee)
logLik(fee)
AIC(fee)
extractAIC(fee)
Z <- predict(fee)
summary(Z)
print(fee$internal$glmfit) # print.vblogit
}
})
#
# tests/ppmmarkorder.R
#
# $Revision: 1.4 $ $Date: 2020/04/30 05:23:52 $
#
# Test that predict.ppm, plot.ppm and plot.fitin
# tolerate marks with levels that are not in alpha order
#
if(ALWAYS) { # locale-dependent?
local({
X <- amacrine
levels(marks(X)) <- c("ZZZ", "AAA")
fit <- ppm(X ~marks, MultiStrauss(c("ZZZ","AAA"), matrix(0.06, 2, 2)))
aa <- predict(fit, type="trend")
bb <- predict(fit, type="cif")
plot(fit)
plot(fitin(fit))
})
}
#
# tests/ppmscope.R
#
# Test things that might corrupt the internal format of ppm objects
#
# $Revision: 1.7 $ $Date: 2022/01/19 09:18:20 $
#
if(ALWAYS) { # dependent on R version?
local({
## (1) Scoping problem that can arise when ppm splits the data
fit <- ppm(bei ~elev, data=bei.extra)
mm <- model.matrix(fit)
## (2) Fast update mechanism
fit1 <- ppm(cells ~x+y, Strauss(0.07))
fit2 <- update(fit1, ~y)
fit3 <- update(fit2, ~x)
## (3) New formula-based syntax
attach(bei.extra)
slfit <- ppm(bei ~ grad)
sl2fit <- update(slfit, ~grad + I(grad^2))
slfitup <- update(slfit, use.internal=TRUE)
sl2fitup <- update(sl2fit, use.internal=TRUE)
## (4) anova.ppm
fut1 <- ppm(cells ~ 1, Strauss(0.1))
futx <- ppm(cells ~ x, Strauss(0.1))
anova(fut1, test="Chi")
anova(futx, test="Chi")
fut1a <- ppm(cells ~ 1, Strauss(0.1), rbord=0)
anova(fut1a, futx, test="Chi")
fut1d <- ppm(cells ~ 1, Strauss(0.1), nd=23)
anova(fut1d, futx, test="Chi")
## This now works!
futxyg <- ppm(cells ~ x + s(y), Strauss(0.1), use.gam=TRUE)
anova(futx, futxyg)
## marked
fatP <- ppm(amacrine ~ marks)
fatM <- ppm(amacrine ~ marks, MultiStrauss(matrix(0.07, 2, 2)))
anova(fatP, fatM, test="Chi")
## (5) expansion of "." in update.ppm
fitb <- ppm(bei ~ . , data=bei.extra)
step(fitb)
})
}
grep#
# tests/ppmtricks.R
#
# Test backdoor exits, hidden options, internals and tricks in ppm
#
# $Revision: 1.19 $ $Date: 2020/04/30 05:23:52 $
#
local({
## (1) skip.border
if(ALWAYS) { # needed below
fit <- ppm(cells, ~1, Strauss(0.1), skip.border=TRUE)
}
## (2) subset arguments of different kinds
if(FULLTEST) {
fut <- ppm(cells ~ x, subset=(x > 0.5))
fot <- ppm(cells ~ x, subset=(x > 0.5), method="logi")
W <- owin(c(0.4, 0.8), c(0.2, 0.7))
fut <- ppm(cells ~ x, subset=W)
fot <- ppm(cells ~ x, subset=W, method="logi")
V <- as.im(inside.owin, Window(cells), w=W)
fet <- ppm(cells ~ x, subset=V)
fet <- ppm(cells ~ x, subset=V, method="logi")
}
## (3) profilepl -> ppm
## uses 'skip.border' and 'precomputed'
## also tests scoping for covariates
if(FULLTEST) {
splants <- split(ants)
mess <- splants[["Messor"]]
cats <- splants[["Cataglyphis"]]
ss <- data.frame(r=seq(60,120,by=20),hc=29/6)
dM <- distmap(mess,dimyx=256)
mungf <- profilepl(ss, StraussHard, cats ~ dM)
mungp <- profilepl(ss, StraussHard, trend=~dM, Q=cats)
}
## (4) splitting large quadschemes into blocks
if(FULLTEST) {
mop <- spatstat.options(maxmatrix=5000)
qr <- quadBlockSizes(quadscheme(cells))
pr <- predict(ppm(cells ~ x, AreaInter(0.05)))
spatstat.options(mop)
qr <- quadBlockSizes(quadscheme(cells))
}
## (5) shortcuts in summary.ppm
## and corresponding behaviour of print.summary.ppm
if(FULLTEST) {
print(summary(fit, quick=TRUE))
print(summary(fit, quick="entries"))
print(summary(fit, quick="no prediction"))
print(summary(fit, quick="no variances"))
}
## (6) suffstat.R
if(ALWAYS) {
fitP <- update(fit, Poisson())
suffstat.poisson(fitP, cells)
fit0 <- killinteraction(fit)
suffstat.poisson(fit0, cells)
}
## (7) various support for class ppm
if(FULLTEST) {
fut <- kppm(redwood ~ x)
A <- quad.ppm(fut)
Z <- as.im(function(x,y){x}, Window(cells))
fitZ <- ppm(cells ~ Z)
U <- getppmOriginalCovariates(fitZ)
}
## (8) support for class profilepl
if(FULLTEST) {
rr <- data.frame(r=seq(0.05, 0.15, by=0.02))
ps <- profilepl(rr, Strauss, cells)
## plot(ps) ## covered in plot.profilepl.Rd
simulate(ps, nrep=1e4)
parameters(ps)
fitin(ps)
predict(ps, type="cif")
}
## (9) class 'plotppm'
if(FULLTEST) {
fut <- ppm(amacrine ~ marks + polynom(x,y,2), Strauss(0.07))
p <- plot(fut, plot.it=FALSE)
print(p)
plot(p, how="contour")
plot(p, how="persp")
}
## (10) ppm -> mpl.engine -> mpl.prepare
if(ALWAYS) { # includes C code
fit <- ppm(cells, NULL)
fit <- ppm(cells ~ x, clipwin=square(0.7))
fit <- ppm(cells ~ x, subset=square(0.7))
DG <- as.im(function(x,y){x+y < 1}, square(1))
fit <- ppm(cells ~ x, subset=DG)
fit <- ppm(cells ~ x, GLM=glm)
fit <- ppm(cells ~ x, famille=quasi(link='log', variance='mu'))
fit <- ppm(cells ~ x, Hardcore(0.07), skip.border=TRUE, splitInf=TRUE)
}
## (11) unidentifiable model (triggers an error in ppm)
if(FULLTEST) {
Q <- quadscheme(cells)
M <- mpl.prepare(Q, cells, as.ppp(Q), trend=~1, covariates=NULL,
interaction=Hardcore(0.3), correction="none")
}
})
reset.spatstat.options()
#
# tests/prediction.R
#
# Things that might go wrong with predict()
#
# $Revision: 1.20 $ $Date: 2020/04/30 05:41:59 $
#
local({
if(ALWAYS) {
## test of 'covfunargs' - platform dependent?
f <- function(x,y,a){ y - a }
fit <- ppm(cells ~x + f, covariates=list(f=f), covfunargs=list(a=1/2))
p <- predict(fit)
## prediction involving 0 * NA
qc <- quadscheme(cells, nd=10)
r <- minnndist(as.ppp(qc))/10
fit <- ppm(qc ~ 1, Strauss(r)) # model has NA for interaction coefficient
p1 <- predict(fit)
p2 <- predict(fit, type="cif", ngrid=10)
stopifnot(all(is.finite(as.matrix(p1))))
stopifnot(all(is.finite(as.matrix(p2))))
}
if(FULLTEST) {
## test of 'new.coef' mechanism
fut <- ppm(cells ~ x, Strauss(0.15), rbord=0)
p0 <- predict(fut, type="cif")
pe <- predict(fut, type="cif", new.coef=coef(fut))
pn <- predict(fut, type="cif", new.coef=unname(coef(fut)))
if(max(abs(pe-p0)) > 0.01)
stop("new.coef mechanism is broken!")
if(max(abs(pn-p0)) > 0.01)
stop("new.coef mechanism gives wrong answer, for unnamed vectors")
#' adaptcoef
a <- c(A=1,B=2,Z=42)
b <- c(B=41,A=0)
ab <- adaptcoef(a, b, drop=TRUE)
}
if(FULLTEST) {
## tests of relrisk.ppm
fut <- ppm(amacrine ~ x * marks)
a <- relrisk(fut, control=2, relative=TRUE)
a <- relrisk(fut, se=TRUE)
a <- relrisk(fut, relative=TRUE, se=TRUE)
fut <- ppm(sporophores ~ marks + x)
a <- relrisk(fut, control=2, relative=TRUE)
a <- relrisk(fut, se=TRUE)
a <- relrisk(fut, relative=TRUE, se=TRUE)
## untested cases of predict.ppm
fit0 <- ppm(cells)
a <- predict(fit0, interval="confidence")
a <- predict(fit0, interval="confidence", type="count")
fit <- ppm(cells ~ x)
b <- predict(fit, se=TRUE, locations=cells)
b <- predict(fit, se=TRUE, interval="confidence")
b <- predict(fit, type="count", se=TRUE)
b <- predict(fit, type="count", window=square(0.5), se=TRUE)
b <- predict(fit, type="count", window=quadrats(cells, 3), se=TRUE)
d <- predict(fit, type="count", interval="prediction", se=TRUE)
d <- predict(fit, type="count", interval="confidence", se=TRUE)
d <- predict(fit, interval="confidence", se=TRUE)
foot <- ppm(cells ~ x, StraussHard(0.12))
d <- predict(foot, ignore.hardcore=TRUE)
dX <- predict(foot, ignore.hardcore=TRUE, locations=cells)
## superseded usages
b <- predict(fit, type="se", getoutofjail=TRUE)
b <- predict(fit, type="se", locations=cells) # warning
b <- predict(fit, total=TRUE)
b <- predict(fit, total=square(0.5))
b <- predict(fit, total=quadrats(cells, 3))
## supporting code
u <- model.se.image(fit, square(0.5))
u <- model.se.image(fit, square(0.5), what="cv")
u <- model.se.image(fit, square(0.5), what="ce")
co <- c(Intercept=5, slope=3, kink=2)
re <- c("Intercept", "slope")
a <- fill.coefs(co, re) # warning
b <- fill.coefs(co, rev(names(co)))
d <- fill.coefs(co, letters[1:3])
## model matrix etc
v <- model.frame(ppm(cells))
fut <- ppm(cells ~ x, Strauss(0.1))
v <- model.matrix(fut, subset=(x<0.5), keepNA=FALSE)
df <- data.frame(x=runif(10), y=runif(10),
Interaction=sample(0:1, 10, TRUE))
m10 <- PPMmodelmatrix(fut, data=df)
mmm <- PPMmodelmatrix(fut, Q=quad.ppm(fut))
#' effectfun for Gibbs
effectfun(fut, "x")
effectfun(fut, "x", se.fit=TRUE)
#' implicit covariate when there is only one
effectfun(fut)
effectfun(fut, se.fit=TRUE)
#' given covariate
dlin <- distfun(copper$SouthLines)
copfit <- ppm(copper$SouthPoints ~ dlin, Geyer(1,1))
effectfun(copfit, "dlin")
effectfun(copfit)
#' covariate that is not used in model
effectfun(fut, "y", x=0)
futS <- ppm(cells ~ 1, Strauss(0.1))
effectfun(futS, "x")
effectfun(futS, "y")
#' factor covariate
fot <- ppm(amacrine~x+marks)
effectfun(fot, "marks", x=0.5, se.fit=TRUE)
#' covariate retained but not used
W <- Window(swedishpines)
a <- solist(A=funxy(function(x,y){x < 20}, W),
B=funxy(function(x,y){factor(x < 20)}, W))
fvt <- ppm(swedishpines ~ A, data=a, allcovar=TRUE)
effectfun(fvt, "A", se.fit=TRUE)
effectfun(fvt, "B", A=TRUE, se.fit=TRUE)
## ppm with covariate values in data frame
X <- rpoispp(42)
Q <- quadscheme(X)
weirdfunction <- function(x,y){ 10 * x^2 + 5 * sin(10 * y) }
Zvalues <- weirdfunction(x.quad(Q), y.quad(Q))
fot <- ppm(Q ~ y + Z, data=data.frame(Z=Zvalues))
effectfun(fot, "y", Z=0)
effectfun(fot, "Z", y=0)
#' multitype
modX <- ppm(amacrine ~ polynom(x,2))
effectfun(modX)
effectfun(modX, "x")
modXM <- ppm(amacrine ~ marks*polynom(x,2))
effectfun(modXM, "x", marks="on")
modXYM <- ppm(amacrine ~ marks*polynom(x,y,2))
effectfun(modXYM, "x", y=0, marks="on")
df <- as.data.frame(simulate(modXM, drop=TRUE))
df$marks <- as.character(df$marks)
dfpr <- predict(modXM, locations=df)
}
})
#
# tests/project.ppm.R
#
# $Revision: 1.7 $ $Date: 2020/04/30 05:41:59 $
#
# Tests of projection mechanism
#
local({
chk <- function(m) {
if(!valid.ppm(m)) stop("Projected model was still not valid")
return(invisible(NULL))
}
if(FULLTEST) {
## a very unidentifiable model
fit <- ppm(cells ~Z, Strauss(1e-06), covariates=list(Z=0))
chk(emend(fit))
## multitype
r <- matrix(1e-06, 2, 2)
fit2 <- ppm(amacrine ~1, MultiStrauss(types=c("off", "on"), radii=r))
chk(emend(fit2))
## complicated multitype
fit3 <- ppm(amacrine ~1, MultiStraussHard(types=c("off", "on"),
iradii=r, hradii=r/5))
chk(emend(fit3))
#' code coverage
op <- spatstat.options(project.fast=TRUE)
fut <- emend(fit, trace=TRUE)
chk(fut)
spatstat.options(op)
#' hierarchical
ra <- r
r[2,1] <- NA
fit4 <- ppm(amacrine ~1, HierStrauss(types=c("off", "on"), radii=r))
chk(emend(fit4))
#' complicated hierarchical
fit5 <- ppm(amacrine ~1, HierStraussHard(types=c("off", "on"),
iradii=r, hradii=r/5))
chk(emend(fit5))
## hybrids
r0 <- min(nndist(redwood))
ra <- 1.25 * r0
rb <- 0.8 * r0
f1 <- ppm(redwood ~1, Hybrid(A=Strauss(ra), B=Geyer(0.1, 2)), project=TRUE)
chk(f1)
f2 <- ppm(redwood ~1, Hybrid(A=Strauss(rb), B=Geyer(0.1, 2)), project=TRUE)
chk(f2)
f3 <- ppm(redwood ~1, Hybrid(A=Strauss(ra), B=Strauss(0.1)), project=TRUE)
chk(f3)
f4 <- ppm(redwood ~1, Hybrid(A=Strauss(rb), B=Strauss(0.1)), project=TRUE)
chk(f4)
f5 <- ppm(redwood ~1, Hybrid(A=Hardcore(rb), B=Strauss(0.1)), project=TRUE)
chk(f5)
f6 <- ppm(redwood ~1, Hybrid(A=Hardcore(rb), B=Geyer(0.1, 2)), project=TRUE)
chk(f6)
f7 <- ppm(redwood ~1, Hybrid(A=Geyer(rb, 1), B=Strauss(0.1)), project=TRUE)
chk(f7)
}
})
reset.spatstat.options()
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.