#######################################
#######################################
#
# KernelRatios.R.txt
#
# Look at ratio of kernel density estimates and K-functions
#######################################
######################################

#Get spatstat package
library(spatstat)

###################################
#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)



########################
#This is for Tim's Maine Data for paper birch and white pine
######################


paperbirch.ppp<-ppp(paperbirch$x,paperbirch$y,xrange=range(paperbirch$x),yrange=range(paperbirch$y))
paperbirch.env<-envelope(paperbirch.ppp)
paperbirch.K<-Kest(paperbirch.ppp)

whitepine.ppp<-ppp(whitepine$x,whitepine$y,xrange=range(whitepine$x),yrange=range(whitepine$y))
whitepine.env<-envelope(whitepine.ppp)
whitepine.K<-Kest(whitepine.ppp)

#kernel density estimates for paper birch
par(mfrow=c(2,2))
par(pty="s")
plot(paperbirch.ppp$x,paperbirch.ppp$y,pch=19,col='red',xlab="X",ylab="Y",main="Paper Birch")
dens1paper<-density.ppp(paperbirch.ppp,sigma=30,dimyx=c(20,20))
dens2paper<-density.ppp(paperbirch.ppp,sigma=50,dimyx=c(20,20))
dens3paper<-density.ppp(paperbirch.ppp,sigma=70,dimyx=c(20,20))
contour(dens1paper,xlab="u",ylab="v",main="Paper Birch Sigma=30",lwd=1.5,col='blue')
contour(dens2paper,xlab="u",ylab="v",main="Paper Birch Sigma=50",lwd=1.5,col='blue')
contour(dens3paper,xlab="u",ylab="v",main="Paper Birch Sigma=70",lwd=1.5,col='blue')

#kernel density estimates for white pine
par(mfrow=c(2,2))
par(pty="s")
plot(whitepine.ppp$x,whitepine.ppp$y,pch=19,col='red',xlab="X",ylab="Y",main="White Pine")
dens1white<-density.ppp(whitepine.ppp,sigma=30,dimyx=c(20,20))
dens2white<-density.ppp(whitepine.ppp,sigma=50,dimyx=c(20,20))
dens3white<-density.ppp(whitepine.ppp,sigma=70,dimyx=c(20,20))
contour(dens1white,xlab="u",ylab="v",main="White Pine Sigma=30",lwd=1.5,col='blue')
contour(dens2white,xlab="u",ylab="v",main="White Pine Sigma=50",lwd=1.5,col='blue')
contour(dens3white,xlab="u",ylab="v",main="White Pine Sigma=70",lwd=1.5,col='blue')

#Ratio of Density Estimates

par(mfrow=c(2,2))
densratio1<-dens3paper
densratio2<-dens3white
#lines below change from INTENSITY estimates to DENSITY estimates and then 
#takes ratio of one to the other  (both ratios are created - paper vs white and
#white vs paper

densratio1$v<-(dens3paper$v/sum(dens3paper$v))/(dens3white$v/sum(dens3white$v))
densratio2$v<-(dens3white$v/sum(dens3white$v))/(dens3paper$v/sum(dens3paper$v))
contour(densratio1,xlab="X",ylab="Y",main="Ratio of PB to WP, Sigma=70",lwd=1.5,col='blue')
contour(densratio2,xlab="X",ylab="Y",main="Ratio of WP to PB, Sigma=70",lwd=1.5,col='blue')
persp(densratio1,theta=45,phi=45,xlab="X",ylab="Y",zlab="lambda",cex.lab=1.5,main="Sigma=70")
persp(densratio2,theta=45,phi=45,xlab="X",ylab="Y",zlab="lambda",cex.lab=1.5,main="Sigma=70")


########################
#This is for 2 CSR from same region
######################

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

x1<-runif(100,0,10)
y1<-runif(100,0,10)
x2<-runif(50,0,10)
y2<-runif(50,0,10)
plot(x1,y1,pch=19,col='red',xlab="",ylab="",main="Cases (Red) and Controls (Blue)")
points(x2,y2,col='blue')
CSR1<-ppp(x1,y1,xrange=range(x1),yrange=range(y1))
CSR2<-ppp(x2,y2,xrange=range(x2),yrange=range(y2))

#kernel density estimates 
dens1<-density.ppp(CSR1,sigma=2,dimyx=c(20,20))
dens2<-density.ppp(CSR2,sigma=2,dimyx=c(20,20))
contour(dens1,xlab="X",ylab="Y",main="Cases, Sigma=2",lwd=1.5,col='blue')
contour(dens2,xlab="X",ylab="Y",main="Controls, Sigma=2",lwd=1.5,col='blue')

#Ratio of Density Estimates
densratio1<-dens1
densratio1$v<-(dens1$v/sum(dens1$v))/(dens2$v/sum(dens2$v))
contour(densratio1,xlab="X",ylab="Y",main="Ratio of Cases to Controls, Sigma=2",lwd=1.5,col='blue')

persp(densratio1,theta=45,phi=45,xlab="X",ylab="Y",zlab="lambda",cex.lab=1.5,main="Sigma=2")