rm(list=ls())
library(INLA)
library(PBSmapping)

# Read in the data
masterDat = read.csv("masterDat.csv")

# use PBS mapping to convert lat /lon to UTM in zone 10
old.coords = cbind(masterDat$BEST_LON_DD, masterDat$BEST_LAT_DD)
coords = old.coords
utm = as.data.frame(coords)
utm$PID = 1
utm$POS = seq(1,dim(utm)[1])
names(utm)[1:2] = c("X","Y")
attr(utm, "zone") <- 10
attr(utm, "projection") ="LL"
utm = convUL(xydata = utm, km = TRUE)
coords = cbind(utm$X,utm$Y)
n <- nrow(coords) # n is number of sets

# Create Year and Site as indexing variables
YEAR = masterDat$YEAR - min(masterDat$YEAR)+1
k = max(YEAR)
SITE = seq(1,n)

####################################################################
# Convert data from data frame to matrix form, and split data into 
# positive / presence-absence sets
####################################################################
y = matrix(NA, n,k)
y2 = matrix(NA, n,k)
yrmat = matrix(0, n,k) 
eff = matrix(NA, n,k)
dep = matrix(NA,n,k)
mon = matrix(NA,n,k)
covar = matrix(NA,n,k)
for(i in 1:length(YEAR)) {
	# There's more efficient ways to do this in R, but for clarity, we're
	# 1) looping over observations
	# 2) splitting the data into 0s/1s and the positive observations
	# 3) recording effort, area swept in meters^2
	# 4) including covariates -- depth and temperature
	# 5) including a design matrix of year effects so that they can be inluded as fixed effects
	y[SITE[i],YEAR[i]] = ifelse(masterDat$Eulachon[i] > 0, masterDat$Eulachon[i], NA)
	y2[SITE[i],YEAR[i]] = masterDat$E01[i]
	eff[SITE[i],YEAR[i]] = masterDat$AREA_SWEPT_MSQ[i]/1000
  	dep[SITE[i],YEAR[i]] = masterDat$BEST_DEPTH_M[i]
  	covar[SITE[i],YEAR[i]] = masterDat$Gtemp_c[i]
  	yrmat[SITE[i],YEAR[i]] = 1 # design matrix for factor year effects
}


# This block of code creates the mesh
bnd = inla.nonconvex.hull(coords, convex=40)		
mesh1 = inla.mesh.2d(loc=coords, boundary = bnd, cutoff=50, offset=c(50,50), max.edge=c(40,40) )	  
# Make a simple plot of the mesh with the actual trawl survey locations ('coords')
plot(mesh1)
points(coords,col="red",cex=0.3)

# Create the data frame for the presence absence 0-1 binomial model
dat <- data.frame(y=as.vector((y2)), time=rep(1:k, each=n), xcoo=rep(coords[,1], k),ycoo=rep(coords[,2], k), haultime = as.vector(eff), depth = as.vector(dep),depth2 = as.vector(dep^2), covar=as.vector(covar), covar2 = as.vector(covar^2), y03 = yrmat[,1], y04 = yrmat[,2], y05 = yrmat[,3], y06 = yrmat[,4],y07 = yrmat[,5], y08 = yrmat[,6], y09 = yrmat[,7], y10 = yrmat[,8], y11 = yrmat[,9], y12 = yrmat[,10])

# Create the data frame for the positive model (we'll model log(y) ~ normal)
datpos <- data.frame(y=as.vector((log(y))), time=rep(1:k, each=n), xcoo=rep(coords[,1], k),ycoo=rep(coords[,2], k), haultime = as.vector(eff), depth = as.vector(dep),depth2 = as.vector(dep^2), covar=as.vector(covar), covar2 = as.vector(covar^2), y03 = yrmat[,1], y04 = yrmat[,2], y05 = yrmat[,3], y06 = yrmat[,4],y07 = yrmat[,5], y08 = yrmat[,6], y09 = yrmat[,7], y10 = yrmat[,8], y11 = yrmat[,9], y12 = yrmat[,10])   

####################################################################
# Prepare data / covariate objects for INLA 
####################################################################
# Create an inla.spde2 model object
spde=inla.spde2.matern(mesh1, alpha=3/2)
iset <- inla.spde.make.index('i', n.spde=spde$n.spde, n.group=k)
   
# Make covariate matrix
X.1 = dat[,6:19]
Covar.names <- colnames(X.1)
XX.list <- as.list(X.1)
effect.list <- list()						
effect.list[[1]] <- c(iset, list(Intercept=1))
for (i in 1:ncol(X.1)) effect.list[[i+1]] <- XX.list[[i]]
names(effect.list) <- c("1", Covar.names)

# Create 
A <- inla.spde.make.A(mesh=mesh1, loc=cbind(dat$xcoo, dat$ycoo),group = dat$time)
A.list = list()
A.list[[1]] = A
for (i in 1:ncol(X.1)) A.list[[i+1]] <- 1

####################################################################
# Do model fitting 
####################################################################
# This specifies the actual formula for the model. Syntax similar to lm() or glm(), but see INLA manual for details
formula1 = as.formula(paste0("y ~ -1 +",  paste(Covar.names, collapse="+"), "+ f(i, model=spde, group=i.group, control.group=list(model='ar1'))"))		# field evolves with AR1 by year

# Fit presence-absence model
sdat <- inla.stack(tag='stdata', data=list(y=dat$y), A=A.list, effects=effect.list)
res1 <- inla(formula1, family = "binomial", data=inla.stack.data(sdat),control.predictor=list(compute=TRUE, A=inla.stack.A(sdat)), verbose = TRUE, debug=TRUE, keep=FALSE,control.compute = list(dic=TRUE, cpo=TRUE,config=TRUE), control.fixed = list(correlation.matrix=TRUE))

# Fit positive model  
sdat <- inla.stack(tag='stdata', data=list(y=datpos$y), A=A.list, effects=effect.list)
res1pos <- inla(formula1, family = "gaussian", offset=log(dat$haultime), data=inla.stack.data(sdat),control.predictor=list(compute=TRUE, A=inla.stack.A(sdat)), verbose = TRUE, debug=TRUE, keep=FALSE,control.compute = list(dic=TRUE, cpo=TRUE,config=TRUE), control.fixed = list(correlation.matrix=TRUE))
 
 
####################################################################
# Illustrate how to get summary statistics
####################################################################
fitted.pos = res1pos$summary.fitted.values$mean
fitted.pres = res1$summary.fitted.values$mean

####################################################################
# Illustrate how to generate samples from posterior
####################################################################
N=20000
pred.01 <- inla.posterior.sample(N, res1)       
pred.pos <- inla.posterior.sample(N, res1pos)