I’m finally biting the bullet and trying to do some basic spatial modeling with my analyses. My 1st test case is a single species (bobolink) as it occurs in my Ch 2 study region (Badlands and Prairies, Potholes and Hardwood Transition). I have created points in the middle of each of my BBS route segments and first off, I need to read them into R.
library(rgdal)
library(spdep)
library(maptools)
library(INLA)
ch2points <- readOGR("CH2POINTS.shp",layer="CH2POINTS")
Earlier in the code, I made a subset of the data frame for my species of interest, and I only want to keep points where they occur.
species_subset <- ch2points[ch2points@data$route %in% vars$route_piece,]
This code will use an adjacency matrix which I’ll generate by tessellating polygons around each point. First, I need to make an identifier to match up the points with the polygons.
species_subset@data$ID <- 1:dim(species_subset@data)[1]
A nice function to create Voronoi polygons is already out there! So, I copied and pasted into my code and it works great. Once the polygons were created (v), I created the adjacency matrix, and generated the identifier to match up the point and polygon layers (but since the order was the same, I ultimately just copied the identifier column I wanted to the polygon layer). Then I translated the adjacency file to an INLA format.
ch2_adj <- poly2nb(v)
v$ID <- 1:length(ch2_adj)
v$route_piece <- species_subset$route
nb2INLA("ch2.adj", ch2_adj)
In order to link the polygons to the original data frame, I matched them by identifiers.
vars$ID <- species_subset$ID[match(vars$route_piece,species_subset$route)]
Now to the real heart of the matter…
model = abundance ~ f(ID,model="besag",graph="ch2.adj") + f(Year,model="rw2") + yearfix + f(rteobs, model="iid") + firstyear
To make the work space as clean as possible, I removed everything but what the model needs…
rm(list= ls()[!(ls() %in% c('vars','model'))])
…before running it! These models are substantially more intense than my previous models, which is to be expected. Still, this model took over 2 hours to run on my laptop!