### NMST543 Spatial Statistics, exercise 1 ###
##############################################


library(spatstat)

############################
### Estimating intensity ###
############################

# Datasets of interest:
help(swedishpines)
plot(swedishpines)

?bei
plot(bei)
plot(bei, pch=20)
plot(bei, pch=".")

data(package="spatstat.data")

### Homogeneous intensity ###
#############################

# Under the assumption of constant intensity:
npoints(swedishpines)/area(Window(swedishpines))
intensity(swedishpines)

# We calculated mean number of points per square unit:
unitname(swedishpines) # 0.0074 trees per square decimeter
rescale(swedishpines, s=10)
intensity(rescale(swedishpines, s=10)) # 0.74 trees per square meter

# Important information about a dataset can be accessed easily:
summary(swedishpines)

# What other information is available?
attributes(swedishpines)
swedishpines$window
swedishpines$n
swedishpines$x
swedishpines$y

# More detailed technical information:
str(swedishpines)


### What if we suspect inhomogeneity? ###
#########################################

# One can e.g. divide observation window into smaller parts (quadrats)
# and count the number of points in each of them. If the pattern is
# homogeneous, the quadrat counts should be approximately the same.

q.swp <- quadratcount(swedishpines, nx=3, ny=3)
q.swp
plot(swedishpines, pch="+")
plot(q.swp, add=TRUE, col="red", cex=1.5, lty=2)

# Note that it is difficult to choose a reasonable number of quadrats
# (bias vs. variability trade-off). Using large quadrats means small
# relative error of the counts but also disregards the spatial variation
# in intensity within each quadrat.


### Quadrat counting test of homogeneity ###
############################################

# For the moment let us assume we observe a Poisson process.
# Null hypothesis: homogeneous Poisson process.
# Alternative hypothesis: inhomogeneous process.

# It is not difficult to develop a chi-square type of test of uniformity:
q.swp.test <- quadrat.test(swedishpines, 3, 3)
q.swp.test

# We plot the observed counts n_j, expected counts e_j and Pearson residuals
# r_j = (n_j - e_j) / sqrt(e_j)
plot(swedishpines, pch=20, cols="gray")
plot(q.swp.test, add=TRUE)


### Smoothing estimation of intensity function ###
##################################################

?density      # generic function
?density.ppp  # method for the objects of the type "ppp" (planar point pattern),
              # not necessary to write explicitely

# Gaussian smoothing kernel is used by default.

# Raw estimate (no edge-correction):
plot(density(swedishpines, edge=FALSE), zlim=c(0,0.015)) 
points(swedishpines, pch="+")

# Uniformly corrected estimate:
plot(density(swedishpines, edge=TRUE, diggle=FALSE), zlim=c(0,0.015)) 
points(swedishpines, pch="+")

# Diggle's correction:
plot(density(swedishpines, edge=TRUE, diggle=TRUE), zlim=c(0,0.015)) 
points(swedishpines, pch="+")

# Other plotting options:
den <- density(swedishpines)
plot(den)
persp(den)
persp(den,theta=10,phi=30)
contour(den, axes = FALSE)

# Selecting the smoothing bandwidth - which value is "the right one"?
plot(density(swedishpines, sigma=5), zlim=c(0,0.026))
plot(density(swedishpines, sigma=10), zlim=c(0,0.026))
plot(density(swedishpines, sigma=15), zlim=c(0,0.026))

# default choice
b <- attr(den,"sigma")
b

# likelihood cross-validation method (assumes inhomogeneous Poisson process):
b1 <- bw.ppl(swedishpines)
b1
plot(b1)
plot(density(swedishpines, sigma=b1))

# Diggle and Berman's mean square error cross-validation method
# (aims at reconstructing the random intensity surface of a Cox process
# of the Neyman-Scott type):
b2 <- bw.diggle(swedishpines)
b2
plot(b2)
plot(density(swedishpines, sigma=b2))


### Spatially adaptive estimation ###
#####################################

# Sometimes it may be useful to consider variable-bandwidth smoothers instead
# of the fixed-bandwidth smoother discussed above (details not covered here).

# If there is an abrupt change in true intensity one can take advantage of
# the Dirichlet-Voronoi tesselation generated by the point pattern:
plot(dirichlet(swedishpines))
points(swedishpines, pch=20, col="red")

# Intensity estimates are the reciprocals of the area of the individual tiles:
plot(adaptive.density(swedishpines, f=1))

# Or one could use only fraction of the points to generate the tesselation
# and use the rest for counting the number of points in each tile:
plot(adaptive.density(swedishpines, f=0.1))

# And repeat several times and average the results:
plot(adaptive.density(swedishpines, f=0.1, nrep=30))

# Now for a dataset with a sudden change in intensity:
sim.pois.pattern <- rpoispp(lambda=function(x,y){20+300*(x>1/2)})
plot(sim.pois.pattern, pch=20)
plot(adaptive.density(sim.pois.pattern, f=0.1, nrep=30))
plot(adaptive.density(sim.pois.pattern, f=0.1, nrep=100))


### Dependence of intensity on a covariate ###
##############################################

# Intensity function may depend on some covariates:
elev <- bei.extra$elev
slope <- bei.extra$grad
par(mfrow = c(3,1), mar=c(1,1,1,1))
plot(bei, pch=".")
plot(elev)
plot(slope)
dev.off()

# Quadrats determined by a covariate:
br <- quantile(elev, probs=(0:4)/4, type=2)
Zcut <- cut(elev, breaks=br, labels=1:4)
plot(Zcut)
textureplot(Zcut)
V <- tess(image=Zcut)

quadratcount(bei, tess=V)
quadrat.test(bei, tess=V)

# Functional dependence lambda(u) = rho( Z(u) ):
plot(rhohat(bei, elev))



### Inhomogeneity vs. interaction ###
#####################################

# It is difficult (impossible?) to distinguish whether the spatial arrangement
# of points is caused by inhomogeneous intensity or interpoint interactions.

# Even for a realization of a stationary process we obtain some estimate =>
# => be careful with interpretation of the results!
sim.thomas.pattern <- rThomas(kappa=10, sigma=0.02, mu=10)
sim.thomas.pattern
plot(density(sim.thomas.pattern))
points(sim.thomas.pattern, pch=20)



### Individual tasks:  ###
##########################

# 1) For the simulated dataset "sim.thomas.pattern" draw several estimates
# of the intensity function using different values of the bandwidth.
# Discuss what is the "correct" choice of bandwidth in this case.

# 2) Using the command "rMatClust" simulate a realization of the Matern cluster
# process in a unit square with on average 10 clusters, on average 10 points in
# a cluster and r = 0.05 (use help to find out the meaning of r). Plot the new
# realization next to the "sim.thomas.pattern". Can you describe the differences?