# 9.5 Data preparation for R2WinBUGS

1. Exercise 9.5 - Download and extract zip folder into your preferred location
2. Set working directory to the extracted folder in R under File - Change dir...
3. First we need to load the packages needed for the exercise
library(R2WinBUGS)
library(sp)
library(spdep)
library(rgdal)

4. Now open the script "SnowyR2winbugs.R" and run code directly from the script

5. Import dataset created previously

df <- read.table("SnowyFinal.txt", sep="\t")
str(df)
df\$IDvar <- substr(df\$id, 2,5)
df\$IDvar <- as.integer(df\$IDvar)
#Clean up and remove missing data for disease status, sex, and age
summary(df\$Status)
df <- subset(df, df\$Status !="Unknown")
df\$Status <- factor(df\$Status)
summary(df\$Status)
#Recode Sex classes and remove NAs
summary(df\$Sex)
df <- subset(df, df\$Sex !="")
df\$Sex <- factor(df\$Sex)
df\$Sex2 <- as.character(df\$Sex)
df\$Sex2[df\$Sex2 == "Male"] <- "1"
df\$Sex2[df\$Sex2 == "Female"] <- "0"
df\$Sex2
#Recode Age classes and remove NAs
summary(df\$Age)
df <- subset(df, df\$Age !="")
df\$Age <- factor(df\$Age)
table(df\$Age)
#Age classes
# 2-5 3 6+ Adult Calf Yearling
# 188 1 79 5 36 30
#Combine ages classes
df\$NewAge <- df\$Age
summary(df\$NewAge)
#Now convert age to numeric with Yearling=0 (baseline) and Adult=1
df\$Age2 <- df\$NewAge
df\$Age2 <- as.character(df\$Age2)
df\$Age2[df\$Age2 == "Adult"] <- "1"
df\$Age2[df\$Age2 == "Yearling"] <- "0"
df\$Age2
df\$KillYear <- as.factor(df\$KillYear)

summary(df\$KillYear)
#2009 2011 2012
23     15     1

6. Now we are going to re-create our spatial grid used in the previous code. Alternatively, we could export the grid as a shapefile and import it here but this may be preferable if sampling grid cell IDs match up.

#Make a spatial data frame of locations after removing outliers
coords<-data.frame(x = df\$X_Coordina, y = df\$Y_Coordina)
utm.crs<-"+proj=utm +zone=12 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"
utm.spdf <- SpatialPointsDataFrame(coords= coords, data = df, proj4string = CRS(utm.crs))

#Now transform projections all to match DEM (i.e., Albers)
Albers.crs <-CRS("+proj=aea +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80
+towgs84=0,0,0,0,0,0,0 +units=m +no_defs")
df.spdf <-spTransform(utm.spdf, CRS=Albers.crs)

#NOTE: We added 3 cells around outer most samples to encompass all disease
#samples until can figure out how to expand grid polygons
snowy.df <- as.data.frame(df.spdf)
str(snowy.df)
minx <- (min(snowy.df\$x)-10830)
maxx <- (max(snowy.df\$x)+10830)
miny <- (min(snowy.df\$y)-10830)
maxy <- (max(snowy.df\$y)+10830)

#Create vectors of the x and y points
x <- seq(from = minx, to = maxx, by = 3610)
y <- seq(from = miny, to = maxy, by = 3610)

#Create a grid of all pairs of coordinates (as a data.frame)
xy <- expand.grid(x = x, y = y)

#Identifiy projection before creating Spatial Points Data Frame
Albers.crs2 <-"+proj=aea +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
#NOTE: Albers.crs2 is needed because SPDF needs different projection command than spTransform above
grid.pts<-SpatialPointsDataFrame(coords= xy, data=xy, proj4string = CRS(Albers.crs2))
#Need to define points as a grid to convert to Spatial Polygons below
gridded(grid.pts) <- TRUE
#Convert grid points to Spatial Polygons in essence converting to a shapefile
gridsp <- as(grid.pts, "SpatialPolygons")
#Now convert gridpts to Spatial Polygons Data Frame for added flexibility in manipulating layer
grid <- SpatialPolygonsDataFrame(gridsp, data=data.frame(id=row.names(gridsp), row.names=row.names(gridsp)))
class(grid)
plot(grid)
7. Compute adjacency matrix and adj, N, sumNumNeigh required by car.normal using sdep package
shape_nb <- poly2nb(grid)
NumCells= length(shape_nb)
num=sapply(shape_nb,length)
sumNumNeigh=length(unlist(shape_nb))

8. Run correlation on covariates to prevent modeling similar covariates
"dem","slope","aspect","ndviMay09","ndviMay25", "ndviJune10","ndviJuly12","ndviJuly28","ndviAug29","ndviSep14", "ndviSep30"")], type="pearson")

9. Prepare values for model inputs

Result<-df\$Status2
Grid_ID<-df\$ID
Sex<-df\$Sex2
Age<-df\$Age2
Dem <- df\$dem
Slope<-df\$slope
Aspect <- df\$aspect
HSI <- df\$HSI
Wat<-df\$water
Dev<-df\$developed
For<-df\$forest
Shru<-df\$shrub
Gras<-df\$grass
Crop<-df\$crop
Wet<-df\$wetland
Maymax <- df\$mayBigmax
Junmax <- df\$juneBigmax
Jul1max <- df\$jul1Bigmax
Jul2max <- df\$jul2Bigmax
Augmax <- df\$augBigmax
Sep1max <- df\$sep1Bigmax
Sep2max <- df\$sep2Bigmax
Maysum <- df\$mayBigsum
Julsum <- df\$julyBigsum
Sepsum <- df\$sepBigsum

10. Compute adjacency matrix and adj, N, sumNumNeigh required by car.normal function in GeoBUGS.

shape_nb <- poly2nb(grid)
NumCells= length(shape_nb)
num=sapply(shape_nb,length)
sumNumNeigh=length(unlist(shape_nb))

11. Define the model in BUGS language

sink("sublettepriors.bug")
cat("

model
{
#Priors for CAR model spatial random effects:
b.car[1:NumCells] ~ car.normal(adj[], weights[], num[], tau.car)
for (k in 1:sumNumNeigh)
{
weights[k] <- 1
}
for (j in 1:NumCells)
{
epsi[j] ~ dnorm(0,tau.epsi)
}

#Other priors
alpha ~ dflat()
beta1 ~ dnorm(0,1.0E-5)
beta2 ~ dnorm(0,1.0E-5)
beta3 ~ dnorm(0,1.0E-5)
beta4 ~ dnorm(0,1.0E-5)
beta5 ~ dnorm(0,1.0E-5)
beta6 ~ dnorm(0,1.0E-5)
beta7 ~ dnorm(0,1.0E-5)
beta8 ~ dnorm(0,1.0E-5)
beta9 ~ dnorm(0,1.0E-5)
tau.car ~ dgamma(1.0,1.0)
tau.epsi ~ dgamma(17.0393,4.1279)

sd.car<-sd(b.car[])
sd.epsi<-sd(epsi[])
lambda <- sd.car/(sd.car+sd.epsi)

for (i in 1 : 339)
{
Result[i] ~ dbern(pi[i])
logit(pi[i]) <- mu[i]
mu[i] <- alpha + beta1*Sex[i] + beta2*Age[i] + beta3*Dem[i] + be ta4*Slope[i] + beta5*Aspect[i] + beta6*HSI[i] +
beta7*Dev[i] + beta8*For[i] + beta9*Gras[i] + b.car[Grid_ID[i]] + epsi[Grid_ID[i]]
}}
# end model
", fill=TRUE)
sink()

# Bundle data
bugs.data <- list(Result=Result, Grid_ID=Grid_ID, NumCells=NumCells, sumNumNeigh=sumNumNeigh, num=num, adj=adj, Sex=Sex, Age=Age, Dem=Dem, Slope=Slope, Aspect=Aspect, HSI=HSI, Dev=Dev, For=For, Gras=Gras)

12. Load initial values

inits1<- list(alpha = 0, beta1 = 0, beta2 = 0, beta3 = 0, beta4 = 0, beta5 = 0, beta6 = 0, beta7 = 0, beta8 = 0, beta9 = 0)
inits2<- list(alpha = 0, beta1 = 0, beta2 = 0, beta3 = 0, beta4 = 0, beta5 = 0, beta6 = 0, beta7 = 0, beta8 = 0, beta9 = 0)
inits3<- list(alpha = 0, beta1 = 0, beta2 = 0, beta3 = 0, beta4 = 0, beta5 = 0, beta6 = 0, beta7 = 0, beta8 = 0, beta9 = 0)

inits<- list(inits1, inits2, inits3)

#Paramters to estimate and keep track of
parameters <- c("alpha","beta1","beta2", "beta3", "beta4", "beta5", "beta6", "beta7", "beta8", "beta9", "lambda")

# MCMC settings
niter <- 150000
nthin <- 20
nburn <- 10000
nchains <- 3

13. Locate WinBUGS by setting path below specifically for the computer used.

bugs.dir<-"C:\\Program Files\\WinBUGS14"

# Do the MCMC stuff from R
out <- bugs(data = bugs.data, inits = inits, parameters.to.save = parameters, model.file = "sublettepriors.bug", n.chains = nchains, n.thin=nthin, n.iter=niter, n.burnin=nburn, debug=TRUE, bugs.directory=bugs.dir)

print(out, 3)