diagnose.R
if(dev.cur() <= 1) {
dd <- getOption("device")
if(is.character(dd)) dd <- get(dd)
dd()
}
oldpar <- par(ask = interactive() &&
(.Device %in% c("X11", "GTK", "windows", "Macintosh")))
par(mfrow=c(1,1))
oldoptions <- options(warn = -1)
#
#######################################################
#
X <- rpoispp(function(x,y) { 1000 * exp(- 4 * x)}, 1000)
plot(X, main="Inhomogeneous Poisson pattern")
fit.hom <- ppm(X, ~1, Poisson())
fit.inhom <- ppm(X, ~x, Poisson())
diagnose.ppm(fit.inhom, which="marks", type="Pearson",
main="Mark plot\nCircles for positive residual mass\nColour for negative residual density")
par(mfrow=c(1,2))
diagnose.ppm(fit.hom, which="marks",
main=c("Wrong model (homogeneous Poisson)", "raw residuals"))
diagnose.ppm(fit.inhom, which="marks",
main=c("Right model (inhomogeneous Poisson)", "raw residuals"))
par(mfrow=c(1,1))
diagnose.ppm(fit.inhom, which="smooth", main="Smoothed residual field")
par(mfrow=c(1,2))
diagnose.ppm(fit.hom, which="smooth",
main=c("Wrong model (homogeneous Poisson)",
"Smoothed residual field"))
diagnose.ppm(fit.inhom, which="smooth",
main=c("Right model (inhomogeneous Poisson)",
"Smoothed residual field"))
par(mfrow=c(1,1))
diagnose.ppm(fit.inhom, which="x")
par(mfrow=c(1,2))
diagnose.ppm(fit.hom, which="x",
main=c("Wrong model (homogeneous Poisson)",
"lurking variable plot for x coordinate"))
diagnose.ppm(fit.inhom, which="x",
main=c("Right model (inhomogeneous Poisson)",
"lurking variable plot for x coordinate"))
par(mfrow=c(1,1))
diagnose.ppm(fit.hom, type="Pearson",main="standard diagnostic plots")
par(mfrow=c(1,2))
diagnose.ppm(fit.hom, main="Wrong model (homogeneous Poisson)")
diagnose.ppm(fit.inhom, main="Right model (inhomogeneous Poisson)")
par(mfrow=c(1,1))
#
#######################################################
# Q - Q PLOTS
#
qqplot.ppm(fit.hom, 40)
#conclusion: homogeneous Poisson model is not correct
title(main="Q-Q plot of smoothed residuals")
qqplot.ppm(fit.inhom, 40) # TAKES A WHILE...
title(main=c("Right model (inhomogeneous Poisson)",
"Q-Q plot of smoothed residuals"))
# conclusion: fitted inhomogeneous Poisson model looks OK
#
#######################################################
#
data(cells)
plot(cells)
fitPoisson <- ppm(cells, ~1, Poisson())
diagnose.ppm(fitPoisson,
main=c("CSR fitted to cells data",
"Raw residuals",
"No suggestion of departure from CSR"))
diagnose.ppm(fitPoisson, type="pearson",
main=c("CSR fitted to cells data",
"Pearson residuals",
"No suggestion of departure from CSR"))
# These diagnostic plots do NOT show evidence of departure from uniform Poisson
qqplot.ppm(fitPoisson, 40)
title(main=c("CSR fitted to cells data",
"Q-Q plot of smoothed raw residuals",
"Strong suggestion of departure from CSR"))
# Q-Q plot DOES show strong evidence of departure from uniform Poisson.
#
fitStrauss <- ppm(cells, ~1, Strauss(r=0.15))
diagnose.ppm(fitStrauss,
main=c("Strauss model fitted to cells data",
"Raw residuals"))
diagnose.ppm(fitStrauss, type="pearson",
main=c("Strauss model fitted to cells data",
"Pearson residuals"))
# next line takes a LOOONG time ...
qqplot.ppm(fitStrauss, 40, type="pearson")
title(main=c("Strauss model fitted to cells data",
"Q-Q plot of smoothed Pearson residuals",
"Suggests adequate fit"))
# Conclusion: Strauss model seems OK
#
#######################################################
#
data(nztrees)
plot(nztrees)
fit <- ppm(nztrees, ~1, Poisson())
diagnose.ppm(fit, type="pearson")
title(main=c("CSR fitted to NZ trees",
"Pearson residuals"))
diagnose.ppm(fit, type="pearson", cumulative=FALSE)
title(main=c("CSR fitted to NZ trees",
"Pearson residuals (non-cumulative)"))
lurking(fit, expression(x), type="pearson", cumulative=FALSE,
splineargs=list(spar=0.3))
# Sharp peak at right is suspicious
qqplot.ppm(fit, 40, type="pearson")
title(main=c("CSR fitted to NZ trees",
"Q-Q plot of smoothed Pearson residuals"))
# Slight suggestion of departure from Poisson at top right of pattern.
par(oldpar)
options(oldoptions)