######################################################################
#  heterogenous_poissonpp.txt
#
#
#  Deeper examination of inhomogenous poisson pp
#  Requires spatstat package and splines package
#
######################################################################

par (mfrow=c(1,1))
library(spatstat)
library(splines)

######################################################################
# First, we do example of linear increasing mean in x and y
######################################################################

#underlying function
func=function(x,y) { (x+y)/20 }

#generate grid locations
x=c(1:20)
y=c(1:20)
z=matrix(0,20,20)
for (i in 1:20){
  for (j in 1:20) {
     z[i,j]=func(i,j)
  } 
}
#view true lambda function
par(mfrow=c(2,2))
image(x,y,z, main="True Underlying Lambda Function")
persp(x,y,z,main="True Underlying Lambda Function")

#generate random points
linearpp=rpoispp(func,win=owin(c(0,20),c(0,20)))
plot(linearpp, cols="blue", pch=19,main="Randomly Generated Data from Lambda Function")
plot(density(linearpp),main="Kernel Smoothed Estimate of Lambda Function")



###  Fit Flat CSR Model and look at residuals and model fit
model1=ppm(linearpp,~1)
print.ppm(model1)
#Get model diagnostics values
residuals(model1)
par(mfrow=c(1,1))
diagnose.ppm(model1, cols="red", pch=19)



#Chi-Square test of goodness of fit
CHSQ=quadrat.test(model1,nx=5,ny=5)
CHSQ
plot(linearpp, cols="blue",pch=19, main="Chi-Square Test Assuming CSR")
plot(CHSQ, add=TRUE, col="red")

#Kolmogorov-Smirnov Test
test1=cdf.test(model1, "x")
plot(test1, lwd=3, cex=.8)
test1=cdf.test(model1, "y")
plot(test1)


###  Now Fit Model in X and Y and look at residuals and model fit
model2=ppm(linearpp,~polynom(x,y,1))
print.ppm(model2)
#Get model diagnostics values
par(mfrow=c(1,1))
diagnose.ppm(model2, cols="red", pch=19)

#Chi-Square test of goodness of fit
CHSQ2=quadrat.test(model2,nx=5,ny=5)
CHSQ2
plot(linearpp, cols="blue",pch=19, main="Chi-Square Test for Linearly Increasing Model")
plot(CHSQ2, add=TRUE, col="red")

#Kolmogorov-Smirnov Test
test2=cdf.test(model2, "y")
plot(test2, lwd=3, cex=.8)
test2=cdf.test(model2, "x")
plot(test2)



###  Now Overfit model with quadratic and look at residuals and model fit
model3=ppm(linearpp,~polynom(x,y,2))
print.ppm(model3)
#Get model diagnostics values
par(mfrow=c(1,1))
diagnose.ppm(model2)

#Chi-Square test of goodness of fit
CHSQ2=quadrat.test(model3,nx=5,ny=5)
CHSQ2
plot(linearpp, cols="blue",pch=19)
plot(CHSQ2, add=TRUE, col="red")

#Kolmogorov-Smirnov Test
test2=cdf.test(model3, "y")
plot(test2)
test2=cdf.test(model3, "x")
plot(test2)




#Compare AIC of all three models
AIC(model1)
AIC(model2)
AIC(model3)

#Likelihood Ratio Test
anova(model1,model2, test="LRT")
anova(model2,model3, test="LRT")





######################################################################
#Next, the conbination of TWO multivariate normal distributions
######################################################################

#generate grid locations
x <- 1:20
y <- 1:20
z <- matrix(0,20,20)

#set means and inverses of covariance matrices
mu1 <- c(3,3)
mu2 <- c(16,14)
sigmasq1 <- 6
sigmasq2 <- 12
Siginv <- matrix(0,2,2)
Siginv[1,1] <- 1/sigmasq1
Siginv[2,2] <- 1/sigmasq2

Siginv2 <- matrix(0,2,2)
Siginv2[1,1] <- 1/60
Siginv2[2,2] <- 1/35

#calculate density function height combination on a 20 x 20 grid
#z[i,j] has intensity at each location.

for (i in 1:20) {
  for (j in 1:20) {
    z[i,j] <- (1/(2*pi))*exp((-1/2)*
              ( (c(x[i],y[j]) - mu1)%*%Siginv%*%(c(x[i],y[j]) - mu1) ) )+
              (1/(2*pi))*exp((-1/2)*
              ( (c(x[i],y[j]) - mu2)%*%Siginv2%*%(c(x[i],y[j]) - mu2) ) )
  }
}


#set plot picture
par(mfcol=c(1,3),pty='s',mar=c(2,2,2,2))

#Goal is to generate 200 points over this region - so we generate 400 points
#from a Poisson distribution with intensity = max intensity over entire region

#generate random locations on 20x20 grid
xrand <- runif(500,min=0,max=20)
yrand <- runif(500,min=0,max=20)

#generate random values between 0 and max intensity
test <- runif(500,min=,max=max(z))

#zvalues are approximate intensity values at 300 random generated points
zval <- 1:500
for (i in 1:500) {
   zval[i] <- z[trunc(xrand[i])+1,trunc(yrand[i])+1]
  }

#finally, keep 200 first points such that random number between 0 and max intensity 
#is less than curve height at that point

plot1x <- xrand[test<zval][1:200]
plot1y <- yrand[test<zval][1:200]

#make plot
par(mfrow=c(1,1), pty="s")
plot(plot1x,plot1y,xlim=c(0,20),ylim=c(0,20),xlab="u",ylab="v",cex.lab=1.5,pch=19,col='red',
     main="Two Multivariate Normals")
contour(x,y,z,lwd=1,lty=3,add=T,drawlabels=F)



#turn into a ppp object
multivar<-ppp(plot1x,plot1y,xrange=c(0,20),yrange=c(0,20))

#a 3-d view of the 'true' underlying functions
persp(z,theta=45,phi=45,xlab="u",ylab="v",zlab="lambda",cex.lab=1.5,main="Sigma=3")


#a kernel smoothing estimates
dens3<-density.ppp(multivar,sigma=3,dimyx=c(40,40))
image(dens3,theta=45,phi=45,xlab="u",ylab="v",zlab="lambda",cex.lab=1.5,main="Sigma=3")


####################################################
#  Now try parametric density estimation



model1=ppm(multivar,trend=~polynom(x,y,3))
print.ppm(model1)

predmod1=predict.ppm(model1)
persp(predmod1,theta=45,phi=45,xlab="u",ylab="v",zlab="lambda",cex.lab=1.5,main="Sigma=3")


#plot true and two estimates next to each other
par(mfrow=c(1,3))
persp(z,theta=45,phi=45,xlab="u",ylab="v",zlab="lambda",cex.lab=1.5,main="True Intensity")
persp(dens3,theta=45,phi=45,xlab="u",ylab="v",zlab="lambda",cex.lab=1.5,main="Kernel Smooth")
persp(predmod1,theta=45,phi=45,xlab="u",ylab="v",zlab="lambda",cex.lab=1.5,main="Inhom Poisson")


#residual analysis : see Baddeley,A.,Turner,R.,Moller,J.andHazelton,M.(2005)Residualanalysisforspatialpointprocesses.JournaloftheRoyalStatisticalSociety,SeriesB67,617–666.
par(mfrow=c(1,1))
diagnose.ppm(model1)


#Chi-Square test of goodness of fit
CHSQ2=quadrat.test(model1,nx=5,ny=5)
CHSQ2
plot(linearpp, cols="blue",pch=19)
plot(CHSQ2, add=TRUE, col="red")

#Kolmogorov-Smirnov Test
test2=cdf.test(model1, "y")
plot(test2, lwd=3,cex=.8)
test2=cdf.test(model1, "x")
plot(test2)






#############################################################
###
##  Situation with known quadratic increasing intensity function
###
############################################################


z2=matrix(0,20,20)
for (i in 1:20) {
  for (j in 1:20) {
    z2[i,j] <- (i^2+j^2)/400
  }
}

#plot of true intensity function lambda(s)
persp(z2,theta=45,phi=45,xlab="u",ylab="v",zlab="lambda",cex.lab=1.5,main="Sigma=3")


#define function
lamb <- function(x,y) { (x^2+y^2)/400 }

#use the rpoispp() to generate the random points
sloped <- rpoispp(lamb, lmax=100,win=owin(c(1,20),c(1,20)))
plot(sloped $x,sloped $y,xlim=c(0,20),ylim=c(0,20),xlab="u",ylab="v",cex.lab=1.5,pch=19,col='red')
contour(x,y,z2,lwd=1,lty=3,add=T,drawlabels=F)

#kernel density estimate
dens3<-density.ppp(temp,sigma=3,dimyx=c(20,20))
persp(dens3,theta=45,phi=45,xlab="u",ylab="v",zlab="lambda",cex.lab=1.5,main="Sigma=3")

#model - I messed around with different corerctions, but doesn't seem to matter much
model1=ppm(sloped ,trend=~polynom(x,y,2), correction = "border")
print.ppm(model1)
predmod1=predict.ppm(model1)
persp(predmod1,theta=45,phi=45,xlab="u",ylab="v",zlab="lambda",cex.lab=1.5,main="Sigma=3")

plot.msr(residuals(model1))
diagnose.ppm(model1)


#see what happens if fit flat surface
model2=ppm(sloped,trend=~1, correction = "border")
print.ppm(model2)
predmod2=predict.ppm(model2)
persp(predmod2,theta=45,phi=45,xlab="u",ylab="v",zlab="lambda",cex.lab=1.5,main="Sigma=3")

plot.msr(residuals(model2))
diagnose.ppm(model2)



#############################################################
#   
#  Poisson Model with Covariate - BEI DATA
#   
#############################################################

#attach dataset for tropical rain forest
data(bei)

#Z is the information on slope (gradient)
Z <- bei.extra$grad
plot(Z, main="BEI Tree Locations vs. Slope")
plot(bei, add=TRUE, pch="+")

########################################
#First fit as linear trend in x and y

par(mfrow=c(1,1))
model1=ppm(bei, ~polynom(x,y,1))
model1
plot(model1, how="image", se=FALSE)
diagnose.ppm(model1)

#Chi-Square test of goodness of fit
CHSQ=quadrat.test(model1,nx=5,ny=5)
CHSQ
plot(bei, cols="blue",pch=19, cex=.5)
plot(CHSQ, add=TRUE, col="red")

#Kolmogorov-Smirnov Test
test2=cdf.test(model1, "y")
plot(test2, lwd=3,cex=.8)


########################################
#Now Fit polynomial in covariates - also includes underlying linear trend

grad=bei.extra$grad
elev=bei.extra$elev

par(mfrow=c(1,1), pty='m')
model2=ppm(bei ~ polynom(grad,elev,2)+ x+y , data=bei.extra)
plot(model2)
diagnose.ppm(model2)
summary(model2)

AIC(model1)
AIC(model2)

#Chi-Square test of goodness of fit
CHSQ2=quadrat.test(model2,nx=10,ny=5)
CHSQ2
plot(bei, cols="blue",pch=19, cex=.5)
plot(CHSQ2, add=TRUE, col="red")

#Kolmogorov-Smirnov Test
test2=cdf.test(model2, "y")
plot(test2, lwd=3,cex=.8)

#Kolmogorov-Smirnov Test
test2=cdf.test(model2, "x")
plot(test2, lwd=3,cex=.8)






#############################################################
#############################################################
#   
#  K and L Functions for Inhomogenous Poisson PP
#   
#############################################################
#############################################################

#get functions for nice pictures
source("http://reuningscherer.net/fes781/rscripts/SpatialPPFuncs.R.txt")


#compare homogenouse L-function to inhomogenous L function
par(mfrow=c(3,2), pty="s")

plot(multivar$x,multivar$y,xlim=c(0,20),ylim=c(0,20),xlab="u",ylab="v",cex.lab=1.5,pch=19,col='red',
     main="Two Multivariate Normals")
contour(x,y,z,lwd=1,lty=3,add=T,drawlabels=F)

plot(sloped$x,sloped$y,xlim=c(0,20),ylim=c(0,20),xlab="u",ylab="v",cex.lab=1.5,pch=19,col='red')
contour(x,y,z2,lwd=1,lty=3,add=T,drawlabels=F)

L(envelope(sloped,func=Kest),Kest(sloped),1,1,"L-plot, Homogenous PP")
L(envelope(multivar,func=Kest),Kest(multivar),1,1,"L-plot, Homogenous PP")

L(envelope(sloped,func=Kinhom),Kinhom(sloped),1,1,"L-plot, Inhomogenous PP")
L(envelope(multivar,func=Kinhom),Kinhom(multivar, sigma=2),1,1,"L-plot, Inhomogenous PP")







##############################################################################
###
###  Tim's Tree Data
###
##############################################################################

###################################
#Tim's MAINE tree data
###################################

# 1. import population and store population size and total stem weight

maine <- read.fwf("http://reuningscherer.net/fes781/data/maine.txt", skip=44,
     col.names=c("treecode","spp","crown", "dbh", "ht", "stemwt", "x", "y"),
     widths=c(7,3,3,8,8,8,8,8), strip.white=T)
maine$ba <- pi*(maine$dbh/2/12)^2

  maine$dbh.cm <- maine$dbh * 2.54;  maine$height.m <- maine$ht * .3048
  maine$ba.m2 <- maine$ba * .3048**2

  maine$tree.label <- 1:nrow(maine); maine$tree.count <- rep(1, nrow(maine))
  maine$stemwt.kg <- maine$stemwt / 2.2046
attach(maine)


N <- nrow(maine)       # store population size
N                      # display population size

spruce <- maine[spp == 1,]  ; spruce$Nsp <- nrow(spruce)
hemlock <- maine[spp == 2,]  ; hemlock$Nhe <- nrow(hemlock)
redmaple <- maine[spp == 3,]  ; redmaple$Nrm <- nrow(redmaple)
cedar <- maine[spp == 4,]  ; cedar$Nce <- nrow(cedar)
balsamfir <- maine[spp == 5,]  ; balsamfir$Nbf <- nrow(balsamfir)
paperbirch <- maine[spp == 6,]  ; paperbirch$Npb <- nrow(paperbirch)
whitepine <- maine[spp == 7,]  ; whitepine$Nwp <- nrow(whitepine)


#make two ppp objects

whitepine.ppp<-ppp(whitepine$x,whitepine$y,xrange=range(whitepine$x),yrange=range(whitepine$y))
paperbirch.ppp<-ppp(paperbirch$x,paperbirch$y,xrange=range(paperbirch$x),yrange=range(paperbirch$y))

plot(whitepine.ppp)
plot(paperbirch.ppp)
#plot kernel smoothed plots

par(pty="s", mfrow=c(1,2))
dens3paper<-density.ppp(paperbirch.ppp,sigma=70,dimyx=c(20,20))
contour(dens3paper,xlab="u",ylab="v",main="Paper Birch Sigma=70",lwd=1.5,col='blue')
dens3white<-density.ppp(whitepine.ppp,sigma=70,dimyx=c(20,20))
contour(dens3white,xlab="u",ylab="v",main="White Pine Sigma=70",lwd=1.5,col='blue')


#######what we saw before : get K and L functions, assuming HOMOGENOUS PP
paperbirch.env<-envelope(paperbirch.ppp)
paperbirch.K<-Kest(paperbirch.ppp)

whitepine.env<-envelope(whitepine.ppp)
whitepine.K<-Kest(whitepine.ppp)

#plot two species
par(mfrow=c(3,2))
par(pty="s")
plot(paperbirch.ppp$x,paperbirch.ppp$y,pch=19,col='red',xlab="X",ylab="Y",main="Paper Birch")
plot(whitepine.ppp$x,whitepine.ppp$y,pch=19,col='red',xlab="X",ylab="Y",main="White Pine")

#plot K and L functions
Kplot(paperbirch.env, paperbirch.K,1,40,"K-plot Paper Birch, Homogenous PP")
Kplot(whitepine.env, whitepine.K,1,40,"K-plot White Pine, Homogenous PP")
L(paperbirch.env,paperbirch.K,275,30,"L-Plot Paper Birch, Homogenous PP")
L(whitepine.env,whitepine.K,275,30,"L-Plot Paper Birch, Homogenous PP")


#######NOW - redo as HETEROGENOUS PP, i.e. use density function, THEN calculate K function, compare results

#plot the results species

par(mfrow=c(2,2))
par(pty="s")

SIGMA=100

plot(paperbirch.ppp$x,paperbirch.ppp$y,pch=19,col='red',xlab="X",ylab="Y",main="Paper Birch")
paperbirch.env2<-envelope(paperbirch.ppp, func=Kinhom, sigma=SIGMA)
paperbirch.K2<-Kinhom(paperbirch.ppp, sigma=SIGMA)
plot(density(paperbirch.ppp,sigma=SIGMA))
#plot L functions, with and without kernel smoothing
L(paperbirch.env,paperbirch.K,275,30,"Paper Birch, Homogenous PP")
L(paperbirch.env2,paperbirch.K2,275,30,"Paper Birch, Heterogenous PP")


plot(whitepine.ppp$x,whitepine.ppp$y,pch=19,col='red',xlab="X",ylab="Y",main="White Pine")
plot(density(whitepine.ppp,sigma=SIGMA))
whitepine.env2<-envelope(whitepine.ppp, func=Kinhom, sigma=SIGMA)
whitepine.K2<-Kinhom(whitepine.ppp, sigma=SIGMA)
#plot L functions, with and without kernel smoothing
L(whitepine.env,whitepine.K,275,30,"White Pine, Homogenous PP")
L(whitepine.env2,whitepine.K2,275,30,"White Pine, Heterogenous PP")


#try same thing with third degree polynomial for white pine
#I DONT THINK THIS IS WORKING!!!!!

par(mfrow=c(1,1))
model1=ppm(whitepine.ppp ,trend=~polynom(x,y,3), correction = "border")
diagnose.ppm(model1)
summary(model1)
whitepine.env3<-envelope(model1, func=Kinhom)
#SOMETHING WRONG WITH STATEMENT BELOW
whitepine.K3<-Kinhom(whitepine.ppp, lambda=predict(model1))
L(whitepine.env,whitepine.K,100,30,"White Pine, Homogenous PP")
L(whitepine.env3,whitepine.K3,100,30,"White Pine, Heterogenous PP")

#############################################################
#############################################################
#   
#  Poison Cluster Processes 
#   
#############################################################
#############################################################


# as an example, we generate a Poisson Cluster Process

library(MASS)  #need to generate multivariate normal


par(mfrow=c(1,1), pty="s")

#Clustered process - Poisson mean 10 centers, Poisson mean 5 points per center, multivariate norm dist of points
centers=rpoispp(.05, win=owin(c(0,20),c(0,20)))
plot(centers, pch=19)
clusters=rpois(centers$n,lambda=5)
mvnorms=mvrnorm(n=sum(clusters), mu=c(0,0), Sigma=matrix(c(.5,0,0,.5),nrow=2))
xvals=rep(centers$x,clusters)+mvnorms[,1]
yvals=rep(centers$y,clusters)+mvnorms[,2]
points(xvals, yvals, col='red', pch=19)
cluster=ppp(xvals,yvals, window=owin(range(xvals), range(yvals)))
plot(cluster, pch=19, main="Clustered Process",xlab='x',ylab='y')
clustmod=kppm(cluster, trend=~1)
summary(clustmod)
#compare k-functions, and show Gaussian kernel fit to clusters
plot(clustmod)
plot(predict(clustmod))


#EXAMPLE redwood data
data(redwood)
plot(redwood$x, redwood$y, pch=19)
mod2=kppm(redwood, trend=~1)
plot(mod2)