Abstract
The objective of this project was to leverage ten years of Melanoplus sanguinipes (Msang) field survey data to assess the response of nymph recruitment under projected climate conditions through the year 2040. To achieve this goal, a multi-level, joint modeling framework was developed that individually assessed nymph and adult life stages while concurrently accounting for observation bias connected to preferential sampling. To ensure that density-dependence between Msang life stages was realistically captured, the model incorporated shared spatiotemporal effects that permitted adult–nymph interactions to inform demographic estimates.
0.1 Open Science Framework (OSF)
Example code, simulated grasshopper occurrences, and additional information are available at OSF: https://osf.io/ny4b5/
0.2 Citation:
Humphreys JM, Srygley RB, Branson DH. Geographic Variation in Migratory Grasshopper Recruitment under Projected Climate Change. Geographies. 2022; 2(1):12-30. https://doi.org/10.3390/geographies2010003
Simplified grasshopper lifecycle.
1 Preliminaries
1.1 Load Packages
suppressMessages(library(INLA))
suppressMessages(library(rgdal))
suppressMessages(library(reshape2))
suppressMessages(library(raster))
suppressMessages(library(rasterVis))
suppressMessages(library(rgeos))
suppressMessages(library(maptools))
suppressMessages(library(mapproj))
suppressMessages(library(spdep))
suppressMessages(library(ggplot2))
suppressMessages(library(ggthemes))
suppressMessages(library(GISTools))
suppressMessages(library(lattice))
suppressMessages(library(gridExtra))
suppressMessages(library(dplyr))
suppressMessages(library(deldir))
suppressMessages(library(sp))
suppressMessages(library(viridis))
suppressMessages(library(pals))
suppressMessages(library(kableExtra))
suppressMessages(library(stringi))1.2 Define Spatial Domain
Study area jurisdictional boundaries.
LL84 = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
nProj = "+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0
+x_0=0 +y_0=0 +k=1 +units=km +nadgrids=@null +no_defs"
States = map("state", fill = TRUE, plot = FALSE)
IDs = sapply(strsplit(States$names, ":"), function(x) x[1])
States = map2SpatialPolygons(States, IDs = IDs, proj4string = CRS(LL84))
pid = sapply(slot(States, "polygons"), function(x) slot(x, "ID"))
p.df = data.frame(ID = 1:length(States), row.names = pid)
States = SpatialPolygonsDataFrame(States, p.df)
States$NAME = rownames(States@data)
West.states = subset(States, NAME == "wyoming")
West.states$NAME = stri_trans_totitle(West.states$NAME)
West.states$Miss = "West"
Domain = West.states1.3 Load Simulated Data
Key data attributes include:
1. Set indicating if location is a grasshopper (GH), grid location from raster (grid), or spatial node (Node) from triangulated mesh reconstructed below.
2. Class Age class of grasshopper or designation as grid or node as decribed in line above.
3. Year Calendar year of grasshopper survey. A value of 0 (zero) indicates data is a node or grid location as described above.
4. Exp.w Geographic exposure, or area of natural neighborhood surronding location.
all.data = read.csv("simulated_gh.csv", stringsAsFactors = FALSE)
dim(all.data)[1] 185915 29unique(all.data$Set)[1] "GH" "grid" "Node"unique(all.data$Class)[1] "Adult" "Nymph" "grid" "Node" unique(all.data$YEAR)[1] 2020 2021 0range(all.data$Exp.w, na.rm = T)[1] 0.0000 174.5256all.data = SpatialPointsDataFrame(all.data[, c("Long", "Lat")], all.data)
proj4string(all.data) = nProj
GH.pnts.p = subset(all.data, Set == "GH")
GH.pnts.pclass : SpatialPointsDataFrame
features : 596
extent : -12361.26, -11585.93, 5015.391, 5620.163 (xmin, xmax, ymin, ymax)
crs : +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=km +nadgrids=@null +wktext +no_defs
variables : 29
names : ID, Set, YEAR, Long, Lat, Count, Class, Exp.w, Domain, NN, cur_clim_LD, clim_LD_2040_126, clim_LD_2040_245, clim_LD_2040_370, clim_LD_2040_585, ...
min values : 8, GH, 2020, -12361.2554370416, 5015.39106268664, 1, Adult, 0, in, 1.56953592522768e-06, -2.86047489912884, -3.87705961002293, -3.77291540223685, -3.81862904696371, -4.04823840040748, ...
max values : 1767, GH, 2021, -11585.9311691453, 5620.16327383931, 30, Nymph, 0, in, 78.1600000219902, 2.10521616116018, 0.900618845429641, 0.734917451059067, 0.802405526310221, 0.888729793142446, ... 1.4 View Observations
State of Wyoming wrt the United States
States.ll = fortify(States)
Domain.ll = fortify(Domain)
ewbrks = seq(-120, -80, 10)
nsbrks = seq(25, 50, 5)
ewlbls = unlist(lapply(ewbrks, function(x) ifelse(x < 0, paste(x, "°W"), ifelse(x >
0, paste(x, "°W"), x))))
nslbls = unlist(lapply(nsbrks, function(x) ifelse(x < 0, paste(x, "°S"), ifelse(x >
0, paste(x, "°N"), x))))
ggplot(States.ll, aes(long, lat, group = group)) + geom_polygon(fill = "gray95",
col = "gray60", size = 0.5, alpha = 0.5) + geom_polygon(data = Domain.ll, aes(long,
lat, group = group), col = "gray80", fill = "gray60", size = 0.5) + xlab(" ") +
ylab(" ") + scale_x_continuous(breaks = ewbrks, labels = ewlbls, expand = c(0,
0)) + scale_y_continuous(breaks = nsbrks, labels = nslbls, expand = c(0, 0)) +
coord_equal() + theme(panel.grid.minor = element_line(colour = "gray90", size = 0.5),
panel.grid.major = element_blank(), panel.background = element_blank(), plot.background = element_blank(),
panel.border = element_blank(), strip.text = element_text(face = "bold", size = 18),
strip.background = element_blank(), axis.text.x = element_text(face = "bold",
size = 22), axis.text.y = element_text(face = "bold", size = 22), legend.title = element_text(size = 17,
face = "bold"), legend.text = element_text(size = 11, face = "bold"), legend.position = "bottom",
plot.title = element_blank())
Panels depict location and relative counts of grasshoppers by year.
GH.plt.ll = spTransform(GH.pnts.p, LL84)
GH.plt.df = GH.plt.ll@data %>%
mutate(Long = GH.plt.ll@coords[, 1], Lat = GH.plt.ll@coords[, 2])
GH.plt.df$Size = Scale(GH.plt.df$Count) * 8
GH.plt.df$Size = GH.plt.df$Size + 0.5
ggplot(Domain.ll, aes(long, lat, group = group, col = Class, size = Size)) + geom_polygon(fill = "gray90",
col = "darkgray", size = 0.01) + geom_point(data = GH.plt.df, aes(Long, Lat,
group = NA), pch = 1, stroke = 0.5, size = GH.plt.df$Size) + xlab(" ") + ylab(" ") +
facet_wrap(~YEAR, ncol = 5) + coord_equal() + theme(panel.grid.minor = element_blank(),
panel.grid.major = element_blank(), panel.background = element_blank(), plot.background = element_blank(),
panel.border = element_blank(), strip.text = element_text(face = "bold", size = 18),
strip.background = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(),
axis.ticks.x = element_blank(), axis.ticks.y = element_blank(), legend.title = element_text(size = 17,
face = "bold"), legend.text = element_text(size = 11, face = "bold"), legend.position = "bottom",
axis.title.x = element_blank(), axis.title.y = element_blank(), plot.title = element_blank())
2 Spatial Mesh
Example spatial triangulation with buffer around exterior. Simulated grasshopper points overlain.
Domain.p = spTransform(Domain, nProj)
Hull = gConvexHull(Domain.p)
Hull = gBuffer(Hull, width = 10)
dom.bnds = inla.sp2segment(Hull)
mesh.dom = inla.mesh.2d(boundary = dom.bnds, loc = GH.pnts.p, cutoff = 10, max.edge = 15,
offset = c(10, 25), min.angle = 30)
mesh.dom$n[1] 6441plot(mesh.dom, main = " ")
plot(Domain, add = T, border = "red", lwd = 2)
plot(GH.pnts.p, col = "red", add = T)
3 Model Construction
3.1 Node Time Series
Duplicate mesh nodes for each year in analysis. A separate copy of mesh nodes are created for each year with grasshopper observations.
gh.N = length(which(all.data@data$Set == "GH"))
node.1 = all.data@data %>%
filter(Set == "Node")
node.1 = node.1[sample(nrow(node.1), gh.N), ]
dim(node.1)[1] 596 29myYears = seq(2020, 2021, 1)
for (i in 1:length(myYears)) {
node.tmp = node.1
node.tmp$YEAR = myYears[i]
if (i == 1) {
time.nodes = node.tmp
} else {
time.nodes = rbind(time.nodes, node.tmp)
}
}
dim(time.nodes)[1] 1192 29dim(time.nodes)[1]/mesh.dom$n[1] 0.1850644unique(time.nodes$YEAR)[1] 2020 2021GH.A = all.data@data %>%
filter(Class == "Adult")
GH.N = all.data@data %>%
filter(Class == "Nymph")
FE.A = rbind(GH.A, time.nodes)
FE.N = rbind(GH.N, time.nodes)3.2 Organize Data
Combing and organizing data as list() objects.
3.2.1 Adult Stage
FE.df = FE.A #Copy for adults
k = length(levels(factor(factor(FE.df$YEAR)))) #Time steps
FE.df$Year.int = as.integer(factor(FE.df$YEAR)) #Year as integer
locs = cbind(FE.df$Long, FE.df$Lat) #projected locations
A.mat = inla.spde.make.A(mesh.dom, #Matrix joining matching to spatial mesh
alpha = 2,
loc=locs,
group = FE.df$Year.int) #indexed by time
A.mat2 = A.mat #copy matrix
spde0 = inla.spde2.pcmatern(mesh.dom, alpha = 2, #Spatial prior
prior.range=c(600, 0.01), # a 1% probability that spatial range is withing 600km
prior.sigma=c(1, 0.01),
constr = TRUE)
Field.1 = inla.spde.make.index("Field.1",
spde0$n.spde,
n.group=k) #Index adults spatial fields by time (Pattern level)
Field.2 = inla.spde.make.index("Field.2",
spde0$n.spde,
n.group=k) #Copy of above index adults spatial fields by time (process level)
FE.df$Exp.sc = FE.df$Exp.w #Exposure: geographic areal extents (natural neighborhoods) for mesh nodes
FE.df$NN2 = round(FE.df$NN/10, 3) #scale
range(FE.df$NN2)[1] 0.000 10.597FE.lst = list(c(Field.1, #spatial index (pattern level)
list(intercept1 = 1)), #dedicated intercept
list(NN.1 = FE.df[,"NN2"],
Steps = FE.df[,"Year.int"],
ID = FE.df[,"ID"]))
FE.df$OBS = ifelse(FE.df$Set == "GH", 1, 0) #detection = 1, non detection = 0
App.stk = inla.stack(data = list(Y = cbind(FE.df$OBS, NA, NA, NA), #four columns response matrix (2 adult levels, 2 nymph levels); adult pattern level here
e = FE.df$Exp.sc), #exposure
A = list(A.mat2, 1), #projection matrix
effects = FE.lst, #variables
tag = "App.0") #label
FE.lst2 = list(c(Field.2, #Process Level spatial index
list(intercept2 = 1)), #dedicated intercept
list(hist.clim1 = FE.df[,"cur_clim_LD"], #climate
pred.soil1 = FE.df[,"soils_LD"], #soils
pred.gHab1 = FE.df[,"gHab_LD"], #vegetation
NN = FE.df[,"NN2"],
Year = FE.df[,"YEAR"],
Year.int = FE.df[,"Year.int"],
Step = FE.df[,"Year.int"],
IDx = FE.df[,"ID"]))
A.resp.stk = inla.stack(data = list(Y = cbind(NA, FE.df$Count, NA, NA), #2nd column = adults counts
e = FE.df$Exp.sc),
A = list(A.mat, 1),
effects = FE.lst2,
tag = "Aresp.0")
Adult.stk = inla.stack(App.stk, A.resp.stk) #combine adult pattern and process levels 3.3 Nymph Stage
Same code used to setup adult levels, but now specific to nymphs.
FE.df = FE.N #Copy for nymphs
k = length(levels(factor(factor(FE.df$YEAR)))) #time steps
FE.df$Year.int = as.integer(factor(FE.df$YEAR)) #year ids
locs = cbind(FE.df$Long, FE.df$Lat) #nymph locations
A.nym = inla.spde.make.A(mesh.dom, #spatial projection matrix
alpha = 2,
loc=locs,
group = FE.df$Year.int) #index by year
A.nym2 = A.nym #copy
Field.1n = inla.spde.make.index("Field.1n",
spde0$n.spde,
n.group=k)
Field.2n = inla.spde.make.index("Field.2n",
spde0$n.spde,
n.group=k)
Field.3a = inla.spde.make.index("Field.3a", #extra index to copy adult level
spde0$n.spde,
n.group=k)
#Nymph Pattern Level
FE.df$Exp.sc = FE.df$Exp.w #exposure
FE.df$NN2 = round(FE.df$NN/10, 3) #scaling
FE.lst = list(c(Field.1n, #Pattern level spatial index for nymphs
list(intercept3 = 1)),
list(NN.1 = FE.df[,"NN2"],
Steps = FE.df[,"Year.int"],
ID = FE.df[,"ID"]))
FE.df$OBS = ifelse(FE.df$Set == "GH", 1, 0) #detection = 1, non detection = 0
Npp.stk = inla.stack(data = list(Y = cbind(NA, NA, FE.df$OBS, NA), # 3rd column is nymph pattern
e = FE.df$Exp.sc), #areal exposure
A = list(A.nym2, 1), #spatial matrix
effects = FE.lst, #variables
tag = "Npp.0") #label
#Nymph Process Level
FE.lst2 = list(c(Field.2n, #Nymph spatial index
Field.3a, #spatial inde to copy adult level
list(intercept4 = 1)),
list(hist.clim2 = FE.df[,"cur_clim_LD"],
pred.soil2 = FE.df[,"soils_LD"],
pred.gHab2 = FE.df[,"gHab_LD"],
NN = FE.df[,"NN2"],
Year = FE.df[,"YEAR"],
Year.int = FE.df[,"Year.int"],
Step = FE.df[,"Year.int"],
IDx = FE.df[,"ID"]))
Nresp.stk = inla.stack(data = list(Y = cbind(NA, NA, NA, FE.df$Count), #4th column is nymph counts
e = FE.df$Exp.sc),
A = list(A.nym, 1),
effects = FE.lst2,
tag = "Nresp.0")
Nym.stk = inla.stack(Npp.stk, Nresp.stk) # combine nymph pattern and process levels 3.4 Join Adult and Nymph Stages
big.stk = inla.stack(Adult.stk, Nym.stk)4 Formula
hc1 = list(theta = list(prior = 'normal', param = c(0, 5))) #normal prior
ctr.g = list(model = 'ar1', hyper = hc1) #order-1 autoregressive prior
Frm1 = Y ~ -1 + intercept1 + #adult pattern intercept
intercept2 + #adult process intercept
intercept3 + #nymph pattern intercept
intercept4 + #nymph process intercept
f(Field.1, #adult pattern spatial index
model=spde0, #spatial prior as defined above
group = Field.1.group, #group by year
control.group=ctr.g) + #prior
f(Field.2, #copy adult pattern to adult process level
copy = "Field.1", #copy this field
group = Field.2.group,
fixed = FALSE,
hyper = hc1) + #noraml prior
f(Field.1n, #nymph pattern
model=spde0,
group = Field.1n.group,
control.group=ctr.g) +
f(Field.2n, #copy nymph pattern to nymph process level
copy = "Field.1n",
group = Field.2n.group,
fixed = FALSE,
hyper = hc1) +
f(Field.3a, # copy adult pattern to nymph process level
copy = "Field.2",
group = Field.3a.group,
fixed = FALSE,
hyper = hc1) +
hist.clim1 + pred.soil1 + pred.gHab1 + #environmental variables specific to adults
hist.clim2 + pred.soil2 + pred.gHab2 #environmental variables specific to nymphs 5 Excute Model
# theta is mean hyperparametrs from prior model run, inclusion speeds up reruns
theta1 = c(3.653, 3.552, 3.101, -0.893, 3.456, 3.246, -0.806, 3.521, 5.517, 5.622, -0.241)
Model.gh = inla(Frm1, #model formula
num.threads = 12, #cores per node to use (turn off for default desktops)
data = inla.stack.data(big.stk), #combined adult and nymph data
family = c("gaussian", "poisson", "gaussian", "poisson"), #pattern levels with gaussian, process level counts are poisson
verbose = FALSE, #show outputs as model runs
E = inla.stack.data(big.stk)$e, #Areal exposure
control.fixed = list(prec = 1, prec.intercept=1), #proper intercepts (fixed priors)
control.predictor = list(
A = inla.stack.A(big.stk), #matrices
compute = TRUE, # include fitted results
link = 1), # link functions
control.mode = list(restart = TRUE, theta = theta1), # use theta above to speed up run
control.inla = list(strategy="adaptive", #strategy to reduce run time
int.strategy = "eb"), # use modes to speed up this example
control.compute=list(dic = TRUE, cpo = TRUE, waic = TRUE)) #mdel indices to estimate