#'      SSAI course on spatial point patterns with spatstat
#'                     Perth, May 2017
#' 
#'                   Lecturer's R script 
#'                  Session 4: Fitting Poisson models
#' 
#'         Copyright (c) Adrian Baddeley and Ege Rubak 2017
#' 
library(spatstat)

plot(rpoispp(100))
plot(rpoispp(30, nsim=20), main.panel="", mar.panel=1, pch=16)
#' [Point process - 'ensemble']
#' [Poisson properties]

#' [Homogeneous Poisson process is completely specified by its intensity]
intensity(swedishpines)
#' Fit model to data
ppm(swedishpines)
#' No other arguments - assume homogeneous Poisson.
m <- ppm(swedishpines)
class(m)
coef(m)
simulate(m)
#'
lam <- function(x,y) { 400 * exp(-2*(x + y)) }
plot(funxy(lam, W = owin()))
plot(Y <- rpoispp(lam))
plot(rpoispp(lam, nsim=20), main.panel="", mar.panel=1, pch=16)
#' [Inhomogeneous Poisson process is completely specified by
#' its intensity function.]

#' Fit model. 
ppm(Y ~ x + y)
#' Formula specifies form of **log** intensity
#' [Write out formula before/after fit]
mxy <- ppm(Y ~ x + y)
coef(mxy)

X <- murchison$gold
L <- murchison$faults
X <- rescale(X, 1000, unitname="km")
L <- rescale(L, 1000, unitname="km")
plot(X, main="Murchison gold data", chars=16, cex=0.6, cols="blue")
plot(L, add=TRUE, col="red")
D <- distfun(L)
plot(D)
plot(L, add=TRUE, col="white")

fitD <- ppm(X ~ D)
fitD
coef(fitD)
#' Fitted coefficients [write formula]
plot(effectfun(fitD, "D"))
confint(fitD)
vcov(fitD)
plot(predict(fitD))
plot(simulate(fitD))

#' test for significance of effect
fit0 <- ppm(X ~ 1)
fit0
anova(fit0, fitD, test="LR")

#' Using the `data` argument of ppm (recommended)
murch <- solist(gold = X, faultdist = D)
ppm(gold ~ faultdist, data = murch)

#'
copP <- rotate(copper$SouthPoints, pi/2)
copL <- rotate(copper$SouthLines, pi/2)
plot(copP, main="Queensland copper data")
plot(copL, add=TRUE)
#'
copD <- distmap(copL)
cfD <- ppm(copP ~ copD)
cf0 <- ppm(copP ~ 1)
anova(cf0, cfD, test="LR")

#' model selection
cfD3 <- ppm(copP ~ polynom(copD, 3))
cfD3
step(cfD3)