#######################################
#######################################
#
#  MarkedSPP.R.txt
#
#  Marked Spatial Point Patterns
#
#  Bivariate and Multi-Group Relationships
#
#######################################
######################################

library(spatstat)


#########################################################
##  Try two random CSR processes and see what happens
#########################################################

#100 random points on square (CSR) in two groups
x<-runif(200,1,100)
y<-runif(200,1,100)
CSR.ppp<-ppp(x,y,xrange=c(1,100),yrange=c(1,100), marks=factor(c(rep(1,100), rep(2,100)), levels=c(1,2), labels=c("A","B")))

#get summary information
summary(CSR.ppp)

#see separate plots
plot(split(CSR.ppp),chars=16,cols="red", main="Plot by Group")

#make plot of all points
plot(CSR.ppp, cols=c(2:3), chars=c(16:17), main="Plot for both Groups")

#NOW.  Start to investigate intensity
plot(density(split(CSR.ppp)),ribbon=FALSE, main="Kernel Intensity Estimates")

#Relative Risk - This is showing the relative risk of seeing a white pine (Case)
#versus cedar (Control - since it was listed first).  Given not as a ratio of density 
#functions but rather as relative probability (i.e. if relative risk is 5, then 
#probability is .5/.1 or 5/6.  If Relative Risk were 1/2, then probabiltiy would be
#.2/.6 or 1/3.  Can get back to relative risk from plot using
#Relative Risk = Prob / (1-Prob)

plot(relrisk(CSR.ppp), main="Relative Probability of each Group")

#As we would expect, no signs of departures from CSR
#G cross function
plot(Gcross(CSR.ppp,"A","B"), lwd=3)
env=envelope(CSR.ppp,fun=Gcross,nsim=99)
plot(env)

#K and L cross function
plot(Kcross(CSR.ppp,"A","B"), lwd=3)
plot(Lcross(CSR.ppp,"A","B"), lwd=3)
env=envelope(CSR.ppp,fun=Kcross,nsim=99)
plot(env)
env=envelope(CSR.ppp,fun=Lcross,nsim=99)
plot(env)


#Flip A and B (i.e. change reference group)
plot(Gcross(CSR.ppp,"B","A"), lwd=3)
env=envelope(CSR.ppp,fun=Gcross, i="B",j="A",nsim=99)
plot(env)

#K and L cross function
plot(Kcross(CSR.ppp,"B","A"), lwd=3)
plot(Lcross(CSR.ppp,"B","A"), lwd=3)
env=envelope(CSR.ppp,fun=Kcross, i="B", j="A" ,nsim=99)
plot(env)
env=envelope(CSR.ppp,fun=Lcross, i="B", j="A" ,nsim=99)
plot(env)


#########################################################
##  Next do a random CSR AND a Regular Grid
#########################################################

#100 random points on square (CSR) AND a Grid
x<-c(runif(100,1,10),rep(c(1:10),10))
y<-c(runif(100,1,10),sort(rep(c(1:10),10)))

CSRGrid.ppp<-ppp(x,y,xrange=c(.5,10.5),yrange=c(.5,10.5), marks=factor(c(rep(1,100), rep(2,100)), levels=c(1,2), labels=c("CSR","GRID")))

#get summary information
summary(CSRGrid.ppp)

#see separate plots
plot(split(CSRGrid.ppp),chars=16,cols="red", main="Plot by Group")

#make plot of all points
plot(CSRGrid.ppp, cols=c(2:3), chars=c(16:17), main="Plot for both Groups")

#NOW.  Start to investigate intensity
plot(density(split(CSRGrid.ppp)),ribbon=FALSE, main="Kernel Intensity Estimates")

#relative risk
plot(relrisk(CSRGrid.ppp), main="Relative Probability of each Group")

#See page 190 of PDF spatial point course for discussion of crossdist (cross distances)
# and nncross (nearest neighbor distance of across dataset).

#this shows grid is closer than we'd expect after distance of .3
plot(Gcross(CSRGrid.ppp,"CSR","GRID"), lwd=3)
env=envelope(CSRGrid.ppp,fun=Gcross,nsim=99)
plot(env)

#interestingly, K and L look just fine
plot(Kcross(CSRGrid.ppp,"CSR","GRID"), lwd=3)
plot(Lcross(CSRGrid.ppp,"CSR","GRID"), lwd=3)
env=envelope(CSRGrid.ppp,fun=Kcross,nsim=99)
plot(env)
env=envelope(CSRGrid.ppp,fun=Lcross,nsim=99)
plot(env)


#notice below when we start at grid, looks just fine
plot(Gcross(CSRGrid.ppp,"GRID","CSR"), lwd=3)
env=envelope(CSRGrid.ppp,fun=Gcross, i="GRID",j="CSR",nsim=99)
plot(env)

plot(Kcross(CSRGrid.ppp,"GRID","CSR"), lwd=3)
plot(Lcross(CSRGrid.ppp,"GRID","CSR"), lwd=3)
env=envelope(CSRGrid.ppp,fun=Kcross, i="GRID", j="CSR",nsim=99)
plot(env)
env=envelope(CSRGrid.ppp,fun=Lcross, i="GRID", j="CSR", nsim=99)
plot(env)


#########################################################
##  Next do a random CSR AND a Clustered Process
#########################################################


#100 random points on square (CSR) AND a grid
x<-runif(100,1,10)
y<-runif(100,1,10)
#10 clusters of 10 points
xa<-runif(100)
ya<-runif(100)
movex<-runif(10,.5,9.5)
movex<-sort(rep(movex,10))+xa
movey<-runif(10,.5,9.5)
movey<-rep(movey,each=10)+ya
x=c(x, movex)
y=c(y,movey)

CSRClust.ppp<-ppp(x,y,xrange=c(.5,10.5),yrange=c(.5,10.5), marks=factor(c(rep(1,100), rep(2,100)), levels=c(1,2), labels=c("CSR","CLUSTER")))

#get summary information
summary(CSRClust.ppp)

#see separate plots
plot(split(CSRClust.ppp),chars=16,cols="red", main="Plot by Group")

#make plot of all points
plot(CSRClust.ppp, cols=c(2:3), chars=c(16:17), main="Plot for both Groups")

#NOW.  Start to investigate intensity
plot(density(split(CSRClust.ppp)),ribbon=FALSE, main="Kernel Intensity Estimates")

#relative risk
plot(relrisk(CSRClust.ppp), main="Relative Probability of each Group")

#See page 190 of PDF spatial point course for discussion of crossdist (cross distances)
# and nncross (nearest neighbor distance of across dataset).

plot(Gcross(CSRClust.ppp,"CSR","CLUSTER"), lwd=3)
env=envelope(CSRClust.ppp,fun=Gcross,nsim=99)
plot(env)
plot(Kcross(CSRClust.ppp,"CSR","CLUSTER"), lwd=3)
plot(Lcross(CSRClust.ppp,"CSR","CLUSTER"), lwd=3)
env=envelope(CSRClust.ppp,fun=Kcross,nsim=99)
plot(env)
env=envelope(CSRClust.ppp,fun=Lcross, nsim=99)
plot(env)

#
plot(Gcross(CSRClust.ppp,"CLUSTER","CSR"), lwd=3)
env=envelope(CSRClust.ppp,fun=Gcross, i="CLUSTER",j="CSR",nsim=99)
plot(env)

plot(Kcross(CSRClust.ppp,"CLUSTER","CSR"), lwd=3)
plot(Lcross(CSRClust.ppp,"CLUSTER","CSR"), lwd=3)
env=envelope(CSRClust.ppp,fun=Kcross, i="CLUSTER", j="CSR",nsim=99)
plot(env)
env=envelope(CSRClust.ppp,fun=Lcross, i="CLUSTER", j="CSR", nsim=99)
plot(env)



###################################
#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 species list
spec=factor(spp, levels=c(1:7),labels=c("spruce","hemlock","redmaple","cedar","balsamfir","paperbirch","whitepine"))

#make marked point process with FACTOR species
maine.ppp<-ppp(maine$x,maine$y,xrange=range(maine$x),yrange=range(maine$y),marks=spec)

#make two species example - white pine and cedar

pcspec=factor(c(rep(1,nrow(cedar)),rep(2,nrow(whitepine))),levels=c(1,2),labels=c("cedar","whitepine"))

pinecedar=rbind(cedar[,-16],whitepine[,-16])
pinecedar.ppp<-ppp(pinecedar$x,pinecedar$y,xrange=range(pinecedar$x),yrange=range(pinecedar$y),marks=pcspec)


#########################################################
##  First do two species and see how this works
#########################################################

#get summary information
summary(pinecedar.ppp)

#see separate plots
plot(split(pinecedar.ppp),chars=16,cols="red", main="Plot by Species")

#make plot of all points
plot(pinecedar.ppp, cols=c(2:3), chars=c(16:17), main="Plot for both Species")

#NOW.  Start to investigate intensity
plot(density(split(pinecedar.ppp)),ribbon=FALSE, main="Kernel Intensity Estimates")

#Relative Risk - This is showing the relative risk of seeing a white pine (Case)
#versus cedar (Control - since it was listed first).  Given not as a ratio of density 
#functions but rather as relative probability (i.e. if relative risk is 5, then 
#probability is .5/.1 or 5/6.  If Relative Risk were 1/2, then probabiltiy would be
#.2/.6 or 1/3.  Can get back to relative risk from plot using
#Relative Risk = Prob / (1-Prob)

plot(relrisk(pinecedar.ppp), main="Relative Probability of White Pine vs. Cedar")

RR=relrisk(pinecedar.ppp)
plot(eval.im(RR > 0.3), main="More than 30 percent White Pine")
points(whitepine$x,whitepine$y, pch=16)

#See page 190 of PDF spatial point course for discussion of crossdist (cross distances)
# and nncross (nearest neighbor distance of across dataset).

plot(Gcross(pinecedar.ppp,"cedar","whitepine"), lwd=3)
env=envelope(pinecedar.ppp,fun=Gcross,nsim=99)
plot(env)

plot(Kcross(pinecedar.ppp,"cedar","whitepine"), lwd=3)
plot(Lcross(pinecedar.ppp,"cedar","whitepine"), lwd=3)
env=envelope(pinecedar.ppp,fun=Kcross, i="cedar", j="whitepine",nsim=99)
plot(env)
env=envelope(pinecedar.ppp,fun=Lcross, i="cedar", j="whitepine", nsim=99)
plot(env)


plot(Gcross(pinecedar.ppp,"whitepine","cedar"), lwd=3)
env=envelope(pinecedar.ppp,fun=Gcross, i="whitepine",j="cedar",nsim=99)
plot(env)

plot(Kcross(pinecedar.ppp,"whitepine","cedar"), lwd=3)
plot(Lcross(pinecedar.ppp,"whitepine","cedar"), lwd=3)
env=envelope(pinecedar.ppp,fun=Kcross, i="whitepine", j="cedar",nsim=99)
plot(env)
env=envelope(pinecedar.ppp,fun=Lcross, i="whitepine", j="cedar", nsim=99)
plot(env)



#########################################################
##  NOW do all species
#########################################################


#get summary information
summary(maine.ppp)

#make plot of all points
plot(maine.ppp, cols=c(1:7), chars=16)

#see separate plots
plot(split(maine.ppp),chars=16,cols="red")

#NOW.  Start to investigate intensity
plot(density(split(maine.ppp)),ribbon=FALSE)

#relative risk
mainerr=relrisk(maine.ppp)
plot(mainerr, main="Maine Tree Relative Risk by Species")

#partition space by species
plot(which.max.im(mainerr), main="Most common species")


##  G/K/L from one to many

plot(Gdot(maine.ppp,"cedar"), lwd=3)
env=envelope(maine.ppp,fun=Gdot, i="cedar",nsim=99)
plot(env)

#reduced simulations to 20 since pretty computationally intensive
plot(Kdot(maine.ppp,"cedar"), lwd=3)
plot(Ldot(maine.ppp,"cedar"), lwd=3)
env=envelope(maine.ppp,fun=Kdot, i="cedar",nsim=20)
plot(env)
env=envelope(maine.ppp,fun=Ldot, i="cedar", nsim=20)
plot(env)