# PMD code for HW 4

1)
graves <- read.table('data/graves.txt', header=T, as.is=T)

a) 
graves1.ppp <- ppp(graves$x, graves$y, window=owin(c(4000,10500),c(3000,10500)))

graves1.env <- envelope(graves1.ppp, Lest, nsim=199, nrank=5)

temp <- data.frame(graves1.env)
graves1.env[,2:5] <- temp[,2:5]-temp$r
plot(graves1.env, legendargs=list(cex=0.7))

b) 
gravepoly <- read.table('data/gravespoly.txt', header=T, as.is=T)
graves2.ppp <- ppp(graves$x, graves$y, poly=gravepoly)

graves2.env <- envelope(graves2.ppp, Lest, nsim=199, nrank=5)

temp <- data.frame(graves2.env)
graves2.env[,2:5] <- temp[,2:5]-temp$r
plot(graves2.env, legendargs=list(cex=0.7), legendpos='topright')

c) 
plot(graves1.ppp)

d) 
graves2.Genv <- envelope(graves2.ppp, Gest, nsim=199, nrank=5)

plot(graves2.Genv, legendargs=list(cex=0.7))

e)
graves2.Fenv <- envelope(graves2.ppp, Fest, nsim=199, nrank=5)

plot(graves2.Fenv, legendargs=list(cex=0.7))

f) 
graves2.genv <- envelope(graves2.ppp, pcf, nsim=199, nrank=5)

temp <- data.frame(graves2.genv)
graves2.genv[,2:5] <- log(temp[,2:5])

# could also use:
graves2.genv <- eval.fv(log(graves2.genv))

plot(graves2.genv, legendargs=list(cex=0.7), main='log g(x)', 
  legendpos='topright')
locator()
# to visually find places where obs crosses outside envelope

g) 
temp <- Kest(graves2.ppp, r=c(0,175))
#  need to include r=0 because spatstat demands it
data.frame(temp)
# temp is an fv object.  Doesn't print nicely.  
lambda <- npoints(graves.ppp)/area.owin(graves.ppp)
lambda*(temp$iso[2] - pi*175^2)
lambda*pi*175^2
lambda*temp$iso[2]

2) Analysis of grave marks
a) 
graves$aff <- factor(graves$aff)
# make sure aff is a factor

graves3.ppp <- ppp(graves$x, graves$y, poly=gravepoly, marks=graves$aff)

plot(density(split(graves3.ppp)))

b,c)
r <- seq(0,1650,50)
nsim <- 199

Kdiff <- Kdiff2 <- matrix(0, nrow=nsim+1, ncol=length(r))
K00 <- Kcross(graves3.ppp, i='0', j='0', r=r)
K10 <- Kcross(graves3.ppp, i='1', j='0', r=r)
K11 <- Kcross(graves3.ppp, i='1', j='1', r=r)

Kdiff[1,] <- K11$iso-K00$iso
Kdiff2[1,] <- K11$iso-K10$iso

for (i in 1:nsim) {
  temp <- rlabel(graves3.ppp)

  K00 <- Kcross(temp, i='0', j='0', r=r)
  K10 <- Kcross(temp, i='1', j='0', r=r)
  K11 <- Kcross(temp, i='1', j='1', r=r)

  Kdiff[i+1,] <- K11$iso-K00$iso
  Kdiff2[i+1,] <- K11$iso-K10$iso
  }

Kd.env <- apply(Kdiff,2,quantile,c(0.025, 0.975))
Kd2.env <- apply(Kdiff2,2,quantile,c(0.025, 0.975))

matplot(r, cbind(Kdiff[1,],t(Kd.env)), lty=c(1,3,3), type='l')
title('K11 - K00')
matplot(r, cbind(Kdiff2[1,],t(Kd2.env)), lty=c(1,3,3), type='l')
title('K11 - K10')

d) 
library(dixon)
dixon(as.matrix(graves[,c('x','y','aff')]))

3) madeup data
To read and setup the data:
madeup <- read.table('madeup.txt', header=T, as.is=T)
madeup.ppp <- ppp(madeup$x, madeup$y, window=owin(c(0,20), c(0,18)))

To look at intensity:
plot(density(madeup.ppp))

To model trend in intensity
madeup.tr <- ppm(madeup.ppp, ~x + y)
print(madeup.tr)

To look at clustering or regularity:
madeup.Lenv <- envelope(madeup.ppp, Lest, nsim=199, nrank=5)
temp <- data.frame(madeup.Lenv)
madeup.Lenv[,2:5] <- temp[,2:5] - temp$r
plot(madeup.Lenv, legendpos='bottomright', main='madeup data', ylab='L(r) - r')

madeup.Genv <- envelope(madeup.ppp, Gest, nsim=199, nrank=5)
madeup.Fenv <- envelope(madeup.ppp, Fest, nsim=199, nrank=5)