- Exercise 9.5 - Download and extract zip folder into your preferred location
- Set working directory to the extracted folder in R under File - Change dir...
-
First we need to load the packages needed for the exercise
library(R2WinBUGS)
library(sp)
library(spdep)
library(rgdal) -
Now open the script "SnowyR2winbugs.R" and run code directly from the script
- Import dataset created previously
df <- read.table("SnowyFinal.txt", sep="\t")
head(df)
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
levels(df$NewAge)<-list(Yearling=c("Calf","Yearling"),Adult=c("2-5","3","6+","Adult"))
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 - 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
#Create vectors of the x and y points
#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)
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) - 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)
adj=unlist(shape_nb)
sumNumNeigh=length(unlist(shape_nb)) - Run correlation on covariates to prevent modeling similar covariates
rcorr.adjust(df[,c("Status","water","developed","wetland","forest","shrub", "grass","crop","HSI",
"dem","slope","aspect","ndviMay09","ndviMay25", "ndviJune10","ndviJuly12","ndviJuly28","ndviAug29","ndviSep14", "ndviSep30"")], ) - 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 - 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)
adj=unlist(shape_nb)
sumNumNeigh=length(unlist(shape_nb)) - 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) - 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 - 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)