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

  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)
    adj=unlist(shape_nb)
    sumNumNeigh=length(unlist(shape_nb))

  8. 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"")], )

  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)
    adj=unlist(shape_nb)
    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)