Overview
Taking steps to account for data collection biases, mediating effects, and variable confounding, we assessed the influence of natural climate oscillations on a 40-year record of grasshopper density in the Western US. Central to our analysis was employing spatially varying coefficients to model time and location-specific variation in grasshopper response to climate. Our results quantitatively demonstrated interannual changes in grasshopper density to be effected by seasonal El Nino/Southern Oscillation (ENSO) and Pacific Decadal Oscillation (PDO) variability and to exhibit spatial asynchrony and non-stationarity such that the relative influence of climate on grasshopper density varied through time and across geographic space.
This script provides code that demonstrates the general workflow and modeling framework used for analysis. Full models as described in the publication were executed using Integrated Laplace Approximation (Håvard Rue, Martino, and Chopin 2009; Lindgren and Rue 2015) and required approximately 10 days of run time on high-performance computers using 24 CPUs at 40GB memory, as well as, licensed PARDISO solver software (Verbosio et al. 2017; Kourounis, Fuchs, and Schenk 2018).
For full description of the study, please see the associated article:
Humphreys JM, Srygley RB, Lawton D, Hudson A, and Branson DH, 2022. Grasshoppers Exhibit Asynchrony and Spatial Non-Stationarity in Response to the El Nino/Southern and Pacific Decadal Oscillations (in review)
For data and code to run this analysis, please see the project Open Science Framework website: https://osf.io/dmyhf/
1 Libraries
Loading required packages for statistical analysis.
suppressMessages(library(rgdal))
suppressMessages(library(rgeos))
suppressMessages(library(maptools))
suppressMessages(library(mapproj))
suppressMessages(library(spdep))
suppressMessages(library(sp))
suppressMessages(library(ggplot2))
suppressMessages(library(GISTools))
suppressMessages(library(dplyr))
suppressMessages(library(viridis))
suppressMessages(library(kableExtra))
suppressMessages(library(dagitty))
suppressMessages(library(ggdag))
suppressMessages(library(INLA))2 Study Area
Defining the spatial domain used in the study.
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 = spTransform(States, nProj)
States$NAME = rownames(States@data)
West.states = subset(States, NAME == "arizona" | NAME == "california" | NAME == "colorado" |
NAME == "idaho" | NAME == "kansas" | NAME == "montana" | NAME == "nebraska" |
NAME == "nevada" | NAME == "new mexico" | NAME == "north dakota" | NAME == "oklahoma" |
NAME == "oregon" | NAME == "south dakota" | NAME == "texas" | NAME == "utah" |
NAME == "washington" | NAME == "wyoming")Western US
Study Area
West.fort = fortify(West.states)
States.fort = fortify(States)
ggplot(States.fort, aes(long, lat, group = group)) + geom_polygon(col = "gray85",
fill = "gray95", size = 0.1) + geom_polygon(data = West.fort, aes(long, lat,
group = group), fill = "gray85", col = "gray65", size = 1) + xlab(" ") + ylab(" ") +
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(),
legend.direction = "horizontal", legend.position = "bottom", strip.text = element_text(size = 12,
face = "bold"), strip.background = element_blank(), legend.key = element_blank(),
legend.key.size = unit(5, "line"), legend.key.width = unit(3, "line"), legend.text = element_text(size = 12,
face = "bold"), legend.title = element_text(size = 18, face = "bold"), axis.title.x = element_blank(),
axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(),
axis.ticks.x = element_blank(), axis.ticks.y = element_blank(), plot.title = element_blank())
3 Hexagon Grid
Constructing a regular, hexagon grid to cover the study domain.
States.buffer = gBuffer(West.states, width = 120, byid = FALSE)
size = 65 #cell center distance km
set.seed(1976)
hex_points = spsample(States.buffer, type = "hexagonal", cellsize = size)
hex_grid = HexPoints2SpatialPolygons(hex_points, dx = size)
hex_grid$States = over(hex_grid, West.states)[, "NAME"]
hex_grid = subset(hex_grid, is.na(States) == FALSE)
# Add a data frame
pid = sapply(slot(hex_grid, "polygons"), function(x) slot(x, "ID"))
p.df = data.frame(ID = 1:length(hex_grid), row.names = pid)
hex_grid = SpatialPolygonsDataFrame(hex_grid, p.df)
hex_grid$Area = gArea(hex_grid, byid = T)
hex_grid$Region = 1:dim(hex_grid@data)[1]View Grid
Buffer.fort = fortify(States.buffer)
hex_grid@data$id = rownames(hex_grid@data)
hex.fort = fortify(hex_grid, region = "id")
hex.fort = left_join(hex.fort, hex_grid@data, by = "id")
ggplot(Buffer.fort, aes(long, lat, group = group)) + geom_polygon(col = "gray85",
fill = "gray95", size = 0.1) + geom_polygon(data = hex.fort, aes(long, lat, group = group),
fill = "tan", alpha = 0.2, col = "gray50", size = 0.1) + geom_polygon(data = West.fort,
aes(long, lat, group = group), fill = "transparent", col = "gray65", size = 1) +
xlab(" ") + ylab(" ") + 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(), legend.direction = "horizontal", legend.position = "bottom",
strip.text = element_text(size = 12, face = "bold"), strip.background = element_blank(),
legend.key = element_blank(), legend.key.size = unit(5, "line"), legend.key.width = unit(3,
"line"), legend.text = element_text(size = 12, face = "bold"), legend.title = element_text(size = 18,
face = "bold"), axis.title.x = element_blank(), axis.title.y = element_blank(),
axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks.x = element_blank(),
axis.ticks.y = element_blank(), plot.title = element_blank())
4 Neighborhood Adjacency
Spatial neighborhood graph identifies hexagon cells with shared boundaries.
nb = poly2nb(hex_grid, queen = TRUE)
nb2INLA("J", nb)
J = inla.read.graph("J")View
plot(hex_grid, col = "transparent", border = "brown")
xy = sp::coordinates(hex_grid)
plot(nb, xy, col = "gray30", cex = 0.25, lwd = 0.1, add = TRUE)
5 Grasshopper Data
Simulated grasshopper densities and environmental data. This analysis will focus on assessing the ONI as a Spatially Varying Coefficient during the winter season. The full model concurrently models the ONI, PDO, temperature, and precipitation as SVCs over a 40 year period across all seasons. Code and and data are available at the OSF: https://osf.io/dmyhf/
data_df = read.csv("gh_demo_data.csv", stringsAsFactors = FALSE, header = TRUE, sep = ",")
datatable(data_df, filter = list(position = "top", clear = TRUE, plain = FALSE),
extensions = c("Buttons", "FixedColumns"), options = list(dom = "Bfrtip", buttons = c("copy",
"csv", "excel", "pdf"), scrollX = TRUE, fixedColumns = TRUE))View Simulated Grasshopper Densities
unique(data_df$Year)[1] 2018 2019 2020match.set = data_df %>%
dplyr::select(Year, Region, Density)
hex.2020 = hex.2019 = hex.2018 = hex.fort
hex.2020$Year = 2020
hex.2019$Year = 2019
hex.2018$Year = 2018
hex.yrs = rbind(hex.2020, hex.2019, hex.2018)
hex.yrs = left_join(hex.yrs, match.set, by = c("Year", "Region"))
ggplot(hex.yrs, aes(long, lat, group = group, fill = Density)) + geom_polygon(col = "gray85") +
scale_fill_viridis(name = "Simulated GH Density", discrete = F, option = "inferno",
direction = -1, na.value = "white") + geom_polygon(data = West.fort, aes(long,
lat, group = group), fill = "transparent", col = "gray65", size = 1) + xlab(" ") +
ylab(" ") + facet_wrap(~Year, ncol = 3) + 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(), legend.position = "bottom", strip.text = element_text(size = 25,
face = "bold"), strip.background = element_blank(), legend.key = element_blank(),
legend.key.size = unit(1, "line"), legend.key.width = unit(3, "line"), legend.text = element_text(size = 12,
face = "bold"), legend.title = element_text(size = 18, face = "bold"), axis.title.x = element_blank(),
axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(),
axis.ticks.x = element_blank(), axis.ticks.y = element_blank(), plot.title = element_blank()) +
guides(fill = guide_colourbar(title.position = "top", title.hjust = 0.5))
6 Directed Acyclic Graph (DAG)
Complete causal analysis is provided in the publication and at the OSF site (https://osf.io/dmyhf/), here we plot the DAG as a general refrence in interpreting data variables.
GH1 = dagitty("dag {
bb=\"0,0,1,1\"
Cult [pos=\"0.738,0.694\"]
D [pos=\"0.409,0.447\"]
E [pos=\"0.311,0.328\"]
ENSO [exposure,pos=\"0.158,0.456\"]
GH_obs [outcome,pos=\"0.843,0.472\"]
GH_true [latent,pos=\"0.643,0.468\"]
Inter [latent,pos=\"0.620,0.627\"]
Intra [latent,pos=\"0.497,0.684\"]
OC [pos=\"0.876,0.198\"]
P [pos=\"0.366,0.657\"]
PDO [exposure,pos=\"0.160,0.166\"]
PH [pos=\"0.845,0.066\"]
Pred [latent,pos=\"0.769,0.570\"]
S [latent,pos=\"0.768,0.238\"]
T [pos=\"0.369,0.086\"]
U [latent,pos=\"0.896,0.369\"]
V [pos=\"0.709,0.085\"]
Cult -> GH_true
D -> GH_true
D -> S
D -> V
E -> GH_true
E -> P
E -> S
E -> T
E -> V
ENSO -> P
ENSO -> T
GH_true -> GH_obs
Inter -> GH_true
Intra -> GH_true
OC -> S
P -> D
P -> GH_true
P -> S
P -> V
PDO -> P
PDO -> T
PH -> S
Pred -> GH_true
S -> GH_true
T -> D
T -> GH_true
T -> S
T -> V
U -> GH_obs
V -> GH_true
V -> S
}")View DAG
tidy_dag = tidy_dagitty(GH1)
ggdag_status(tidy_dag, layout = "nicely") + scale_color_manual(values = c("tan",
"gray85", "orange")) + geom_dag_text(col = "gray10") + theme_dag()
7 Statistical Model
The joint likelihood model included a binomial distribution for survey locations and a gamma distribution for density (grasshoppers per square meter).
\[\begin{align} \textit{y}_{1st} \sim \text{Binomial}(\pi_{st}) \quad \& \quad \textit{y}_{2st} \sim \text{Gamma}(\textit{a}_{st}, \textit{b}_{st}) \nonumber \\ \nonumber \end{align}\]
To bivariate response variable was coded such that human observation of grasshoppers (field survey) approximated sampling probability, with
\[\begin{equation} \textit{y}_{1st} = \begin{cases} 1, & \text{if field survey conducted at location (s) and time (t)} \\ 0, & \text{location not surveyed (s,t)} \nonumber \\ \end{cases} \end{equation}\]
By comparison, the response variable’s second column followed a Gamma distribution with the grasshopper density estimate used when the location was surveyed and all other values truncated (removed),
\[\begin{equation} \textit{y}_{2st} = \begin{cases} NA, & \text{if location was not surveyed (s,t)} \\ density, & \text{if a density estimate was reported with survey (s,t). } \nonumber \\ \end{cases} \end{equation}\]
Each part of the bivariate response was linked to a dedicated regression component. The observation process included an intercept (\(\alpha_{1}\)), a spatiotemporal random effect following an order-1 autoregressive prior (\(\zeta_{\textit{st}}^{1}\)), an effect to capture location specific variability (\(\varphi_{s}^1\)), and a shared spatiotemporal process (\(\textit{S}_{\textit{st}}\)) to control for correlation between the two model levels. The coefficient \(\lambda\) serves as a scaler to qunatify the strength of between tier correlation. This can be generalized as,
\[\begin{equation} \text{logit}(\pi_{st}) = \alpha^{1} + \zeta_{\textit{st}}^{1} + \varphi_{s}^1 + \lambda \textit{S}_{\textit{st}} \\ \end{equation}\]
The density level of the model similarly included an intercept, spatiotemporal effect for latencies, and a random effect to quantify variation at the level of the location. Additionally, sampling effort (Effort) was considered, other environmental factors (\(\beta_{2}X_{st}\)), and temporarily dynamic SVCs (\(\sum_{k}^{m} \textit{f} \hspace{0.1cm} (\xi_{k\textit{st}} \hspace{0.1cm} SVC_{kst}))\)).
\[\begin{equation} \text{log}(\mu_{st}) = \alpha^{2} + \beta_{1}\textit{Effort}_{st} + \beta_{2}X_{st} + \sum_{k}^{m} \textit{f} \hspace{0.1cm} (\xi_{k\textit{st}} \hspace{0.1cm} SVC_{kst}) + \zeta_{\textit{st}}^{2} + \varphi_{s}^2 + \textstyle\frac{1}{\lambda}\textit{S}_{\textit{st}} \\ \end{equation}\]
8 Organize Data
Create Year Index as integer.
data_df$Year1 = as.integer(as.factor(data_df$Year))
range(data_df$Year)[1] 2018 2020range(data_df$Year1)[1] 1 3Observation Level
Because not all variables are of the same format and dimensions, data is organized as a list object. For the observation part of the data, this is called the “observation stack” and shown as observation.stk.
obs.var = list(list(alpha.1 = rep(1, dim(data_df)[1])), #intercept
list(t.1 = data_df[,"Year1"], #time steps
zeta.1 = data_df[,"Region"], #spatiotemporal field
varphi.1 = data_df[,"Region"])) #individual variation
data_df$y.1 = ifelse(data_df$Density > 0, 1, 0) #Survey conducted = 1, not conducted = 0
data_df$y.1[is.na(data_df$y.1)] = 0
sum(data_df$y.1) [1] 3350observation.stk = inla.stack(data = list(Y = cbind(data_df$y.1, NA)), #y.1 id first part of the bivariate response
A = list(1,1),
effects = obs.var,
tag = "observation")Grasshopper Density Level
Similar to the observation stack, data for grasshopper density is organized as a list. These variable names/labels correspond to those presented in the publication and as shown in the Directed Acyclic Graph (DAG) used as a causal analysis heuristic.
dens.var = list(list(alpha.2 = rep(1, dim(data_df)[1])), #intercept
list(Effort = data_df[,"Effort"], #sampling effort
E = data_df[,"s_DEM"], #elevation
PH = data_df[,"s_PH"], #soil pH
OC = data_df[,"s_CARB"], #soil organic carbon
Cult = data_df[,"s_pCult"], #cultivated area
V = data_df[,"ndvi.win"], #vegetation
D = data_df[,"spei.win"], #drought
ONI = data_df[,"oni.win"], #winter ONI index from NOAA
PDO = data_df[,"pdo.win"], #winter PDO index
T = data_df[,"temp.win"], #winter temperature
P = data_df[,"ppt.win"], #winter precipitation
t.2a = data_df[,"Year1"], #time steps
t.2b = data_df[,"Year1"], #time steps (extra copy)
t.2c = data_df[,"Year1"], #time steps (extra copy)
zeta.2 = data_df[,"Region"], #spatiotemporal effect for latencies
S = data_df[,"Region"], #for copy of shared spatial field
xi = data_df[,"Region"], #Spatiaotemporal SVC index
varphi.2 = data_df[,"Region"])) #individual variation
data_df$y.2 = data_df$Density
range(data_df$y.2, na.rm=T) #un-surveyed locations are NA [1] 0.4180602 25.0836120density.stk = inla.stack(data = list(Y = cbind(NA, data_df$y.2)), #second part of bivariate response
A = list(1,1),
effects = dens.var,
tag = "density")Combine Observation and Density Levels
Data from both the observation and desnity parts of the model are combined.
Joint.stk = inla.stack(observation.stk, density.stk)9 Run Model
Check bivariate response Variable. Check 12 rows at random. Region lists cell ID numbers, Year is survey year, y.1 indicates if a
survey was conducted (1 = yes, 0 = no) and y.2 is grasshopper density (if surveyed).
set.seed(1)
data_df[sample(nrow(data_df), 12), ] %>%
dplyr::select(Region, Year, y.1, y.2) %>%
knitr::kable(caption = "Grasshopper densities at un-surveyed locations (y.1 = 0) are unknown (y.2 = NA).")| Region | Year | y.1 | y.2 | |
|---|---|---|---|---|
| 1017 | 1017 | 2018 | 1 | 0.8361204 |
| 4775 | 85 | 2020 | 0 | NA |
| 2177 | 2177 | 2018 | 1 | 13.3779264 |
| 5026 | 336 | 2020 | 0 | NA |
| 1533 | 1533 | 2018 | 0 | NA |
| 4567 | 2222 | 2019 | 1 | 5.4347826 |
| 2347 | 2 | 2019 | 0 | NA |
| 270 | 270 | 2018 | 1 | 2.5083612 |
| 4050 | 1705 | 2019 | 1 | 0.8361204 |
| 5307 | 617 | 2020 | 0 | NA |
| 3379 | 1034 | 2019 | 0 | NA |
| 4065 | 1720 | 2019 | 0 | NA |
Model Formula
norm.prior = list(theta=list(prior = "normal", param=c(0, 3)))
pcprior1 = list(prec = list(prior="pc.prec", param = c(1, 0.001))) #priors
Frm.j = Y ~ -1 + alpha.1 + #intercepts
alpha.2 +
f(zeta.1, #observation level spatiotemporal for latencies
model="besag", #Besag model
graph = J, #adjacency matrix
group = t.1, #temporal structure (i.e., separate realization for each time step)
scale.model=TRUE,
constr=TRUE, #constrain field to zero
control.group = list(model = "ar1")) + #autoregressive prior
f(zeta.2, #density level spatiotemporal for latencies
model="besag",
graph = J,
group = t.2a ,
scale.model=TRUE,
constr=TRUE,
control.group = list(model = "ar1")) +
f(S, #shares field, copy from observation level
copy = "zeta.1",
fixed = FALSE,
group = t.2b) +
f(xi, ONI, #spatiotemporal SVC weighted by the ONI index
model="besag",
graph = J,
group = t.2c,
scale.model=TRUE,
constr=FALSE, #do not constrain field to zero
control.group = list(model = "ar1")) +
f(varphi.1, #individual variation
model="iid", #independent and identically distributed
constr=TRUE,
hyper=norm.prior) +
f(varphi.2,
model="iid",
constr=TRUE,
hyper=norm.prior) +
E + Effort + PH + OC + Cult + V + D #environmental effects and effortRun Model
theta.1 = c(3.31254482, -3.07086277, 2.29563236, 0.45380020, 0.46629706, 2.48823476, 1.18624377, 3.10314798, 3.10314798, 0.04833974)
Model.1 = inla(Frm.j, #formula
num.threads = 12, #number of CPUs
data = inla.stack.data(Joint.stk), #data
family = c("binomial", "gamma"), #distributions
verbose = FALSE, #print progress to screen
control.fixed = list(prec = 1, prec.intercept=1), #intercept priors
control.predictor = list(
A = inla.stack.A(Joint.stk), #projection matrix (placeholder only here)
compute = TRUE, #estimate fitted values
link = 1), #default link
control.mode = list(restart = TRUE, theta = theta.1), #above means from previous run to speed-up
control.inla = list(strategy="adaptive", #strategy to speed up
int.strategy = "eb"), #use mode to reduce resources/run time
control.compute=list(dic = FALSE, cpo = FALSE, waic = FALSE)) #comparative indices to calculate10 Results
Although greatly reduced and simplified, the above model may take a long time to run on desktops and laptops (about 90min). A copy of the model outputs and results (“RES_demo.RData”) are provided with the data and code at the OSF site: https://osf.io/dmyhf/
Load previoulsy run model above:
load("RES_demo.RData")Pull Fitted Estimates
idat = inla.stack.index(Joint.stk, "density")$data
data_df$Fitted = exp(Model.1$summary.fitted.values$mean[idat])
cor(data_df$Density, data_df$Fitted, use = "pairwise.complete")[1] 0.995774As done prior to running model, check bivariate response variable; but this time, in relation to fitted estimates. Check 12 rows at random.
set.seed(1)
data_df[sample(nrow(data_df), 12), ] %>%
dplyr::select(Region, Year, y.1, y.2, Fitted) %>%
knitr::kable(caption = "Grasshopper fitted densities (un-surveyed locations were unknown (y.2 = NA).")| Region | Year | y.1 | y.2 | Fitted | |
|---|---|---|---|---|---|
| 1017 | 1017 | 2018 | 1 | 0.8361204 | 0.9070343 |
| 4775 | 85 | 2020 | 0 | NA | 2.6931593 |
| 2177 | 2177 | 2018 | 1 | 13.3779264 | 11.2886559 |
| 5026 | 336 | 2020 | 0 | NA | 1.3914606 |
| 1533 | 1533 | 2018 | 0 | NA | 8.4534501 |
| 4567 | 2222 | 2019 | 1 | 5.4347826 | 5.2637420 |
| 2347 | 2 | 2019 | 0 | NA | 1.6341613 |
| 270 | 270 | 2018 | 1 | 2.5083612 | 2.3980387 |
| 4050 | 1705 | 2019 | 1 | 0.8361204 | 0.8879924 |
| 5307 | 617 | 2020 | 0 | NA | 1.3211859 |
| 3379 | 1034 | 2019 | 0 | NA | 2.8267923 |
| 4065 | 1720 | 2019 | 0 | NA | 0.9167634 |
View random effect parameters (point estimates, mean)
kable(Model.1$summary.hyperpar[, c(1:3, 5)], caption = "Hyperparameters", digits = 2) %>%
kable_styling("striped", full_width = F) %>%
row_spec(0, font_size = 20) %>%
column_spec(1, bold = T)| mean | sd | 0.025quant | 0.975quant | |
|---|---|---|---|---|
| Precision parameter for the Gamma observations[2] | 25.68 | 0.27 | 25.13 | 26.22 |
| Precision for zeta.1 | 0.04 | 0.00 | 0.03 | 0.04 |
| GroupRho for zeta.1 | 0.82 | 0.01 | 0.80 | 0.84 |
| Precision for zeta.2 | 1.58 | 0.06 | 1.47 | 1.71 |
| GroupRho for zeta.2 | 0.27 | 0.05 | 0.19 | 0.38 |
| Precision for xi | 13.68 | 2.61 | 10.04 | 20.11 |
| GroupRho for xi | 0.93 | 0.04 | 0.84 | 0.98 |
| Precision for varphi.1 | 10.75 | 2.47 | 6.57 | 16.18 |
| Precision for varphi.2 | 26.35 | 4.39 | 18.69 | 35.93 |
| Beta for S | 0.04 | 0.01 | 0.03 | 0.06 |
10.0.1 Autocorrelation (\(\rho\))
Autocorrelation parameter for density as a demonstration.
names(Model.1$marginals.hyperpar) [1] "Precision parameter for the Gamma observations[2]"
[2] "Precision for zeta.1"
[3] "GroupRho for zeta.1"
[4] "Precision for zeta.2"
[5] "GroupRho for zeta.2"
[6] "Precision for xi"
[7] "GroupRho for xi"
[8] "Precision for varphi.1"
[9] "Precision for varphi.2"
[10] "Beta for S" ggplot(as.data.frame(Model.1$marginals.hyperpar[["GroupRho for zeta.2"]]), aes(x,y)) +
geom_line(size = 1, linetype = "solid") +
geom_vline(xintercept = 0,
linetype = "solid",
col = "black",
size = 0.5) +
ylab(" ") +
xlab(" ") +
ggtitle("Group Rho for Density \n(spatiotemporal correlation)") +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5),
legend.position = "bottom", #c(0.8, 0.8),
legend.direction = "vertical",
legend.key.width = unit(3,"line"),
legend.title=element_blank(),
legend.text=element_text(size=18),
strip.background = element_blank(),
title = element_text(face="bold", size=22, hjust=0.5),
strip.text = element_text(face="bold", size = 22),
axis.title.y = element_text(face="bold", size = 25),
axis.text.x = element_text(face="bold", size=20, vjust=0.5,
hjust=0.5, angle=0),
axis.text.y = element_text(face="bold", size=20),
axis.title.x = element_text(face="bold", size = 25)) 
10.0.2 Shared Field Interaction (\(\lambda\))
A positive relationship.
ggplot(as.data.frame(Model.1$marginals.hyperpar[["Beta for S"]]), aes(x,y)) +
geom_line(size = 1, linetype = "solid") +
geom_vline(xintercept = 0,
linetype = "solid",
col = "black",
size = 0.5) +
ylab(" ") +
xlab(" ") +
ggtitle("Observtaion-Density Interaction") +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5),
legend.position = "bottom", #c(0.8, 0.8),
legend.direction = "vertical",
legend.key.width = unit(3,"line"),
legend.title=element_blank(),
legend.text=element_text(size=18),
strip.background = element_blank(),
title = element_text(face="bold", size=22, hjust=0.5),
strip.text = element_text(face="bold", size = 22),
axis.title.y = element_text(face="bold", size = 25),
axis.text.x = element_text(face="bold", size=20, vjust=0.5,
hjust=0.5, angle=0),
axis.text.y = element_text(face="bold", size=20),
axis.title.x = element_text(face="bold", size = 25)) 
10.0.3 Individual levels variation (\(\varphi_{s}^2\))
density.iid = Model.1$summary.random$varphi.2[, c(1:4, 6)]
names(density.iid) = c("ID", "Mean", "sd", "Q0.025", "Q0.975")
ggplot(density.iid, aes(ID, Mean)) + geom_linerange(aes(ymin = Q0.025, ymax = Q0.975),
colour = "gray80", lwd = 0.15) + geom_point(shape = 1, size = 2, col = "orange",
alpha = 1) + geom_hline(yintercept = 0, linetype = "solid", col = "gray80", size = 0.5) +
ylab("Effect Size (log)") + xlab("Location (grid cell)") + theme_classic() +
theme(plot.title = element_text(hjust = 0.5), legend.position = "none", legend.direction = "vertical",
legend.key.width = unit(3, "line"), legend.title = element_blank(), legend.text = element_text(size = 18),
strip.background = element_blank(), title = element_text(face = "bold", size = 22,
hjust = 0.5), strip.text = element_text(face = "bold", size = 22), axis.title.y = element_text(face = "bold",
size = 25), axis.text.x = element_text(face = "bold", size = 10, vjust = 0.5,
hjust = 0.5, angle = 0), axis.text.y = element_text(face = "bold", size = 12),
axis.title.x = element_text(face = "bold", size = 25))
10.0.4 Environmental Effects
Plotting a couple for demonstrative purposes.
names(Model.1$marginals.fixed)[1] "alpha.1" "alpha.2" "E" "Effort" "PH" "OC" "Cult"
[8] "V" "D" ggplot(as.data.frame(Model.1$marginals.fixed[["E"]]), aes(x,y)) +
geom_line(size = 1, linetype = "solid") +
geom_vline(xintercept = 0,
linetype = "solid",
col = "black",
size = 0.5) +
ylab(" ") +
xlab(" ") +
ggtitle("Elevation (E)") +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5),
legend.position = "bottom", #c(0.8, 0.8),
legend.direction = "vertical",
legend.key.width = unit(3,"line"),
legend.title=element_blank(),
legend.text=element_text(size=18),
strip.background = element_blank(),
title = element_text(face="bold", size=22, hjust=0.5),
strip.text = element_text(face="bold", size = 22),
axis.title.y = element_text(face="bold", size = 25),
axis.text.x = element_text(face="bold", size=20, vjust=0.5,
hjust=0.5, angle=0),
axis.text.y = element_text(face="bold", size=20),
axis.title.x = element_text(face="bold", size = 25))
ggplot(as.data.frame(Model.1$marginals.fixed[["V"]]), aes(x,y)) +
geom_line(size = 1, linetype = "solid") +
geom_vline(xintercept = 0,
linetype = "solid",
col = "black",
size = 0.5) +
ylab(" ") +
xlab(" ") +
ggtitle("Vegetation (V)") +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5),
legend.position = "bottom", #c(0.8, 0.8),
legend.direction = "vertical",
legend.key.width = unit(3,"line"),
legend.title=element_blank(),
legend.text=element_text(size=18),
strip.background = element_blank(),
title = element_text(face="bold", size=22, hjust=0.5),
strip.text = element_text(face="bold", size = 22),
axis.title.y = element_text(face="bold", size = 25),
axis.text.x = element_text(face="bold", size=20, vjust=0.5,
hjust=0.5, angle=0),
axis.text.y = element_text(face="bold", size=20),
axis.title.x = element_text(face="bold", size = 25)) 
10.0.5 Fitted Estimates
match.set = data_df %>%
dplyr::select(Year, Region, Fitted)
hex.2020 = hex.2019 = hex.2018 = hex.fort
hex.2020$Year = 2020
hex.2019$Year = 2019
hex.2018$Year = 2018
hex.yrs = rbind(hex.2020, hex.2019, hex.2018)
hex.yrs = left_join(hex.yrs, match.set, by = c("Year", "Region"))
ggplot(hex.yrs, aes(long, lat, group = group, fill = Fitted)) + geom_polygon(col = "gray85") +
scale_fill_viridis(name = "Estimated GH Density", discrete = F, option = "inferno",
direction = -1, na.value = "white") + geom_polygon(data = West.fort, aes(long,
lat, group = group), fill = "transparent", col = "gray65", size = 1) + xlab(" ") +
ylab(" ") + facet_wrap(~Year, ncol = 3) + 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(), legend.position = "bottom", strip.text = element_text(size = 25,
face = "bold"), strip.background = element_blank(), legend.key = element_blank(),
legend.key.size = unit(1, "line"), legend.key.width = unit(3, "line"), legend.text = element_text(size = 12,
face = "bold"), legend.title = element_text(size = 18, face = "bold"), axis.title.x = element_blank(),
axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(),
axis.ticks.x = element_blank(), axis.ticks.y = element_blank(), plot.title = element_blank()) +
guides(fill = guide_colourbar(title.position = "top", title.hjust = 0.5))
10.0.6 ONI SVC
Only three years were used in simulating data, therefore just plotting results to demonstrate some variation within and between years.
oni.SVC = Model.1$summary.random$xi[, c(1:4, 6)]
names(oni.SVC) = c("Region", "Mean", "sd", "Q0.025", "Q0.975")
oni.SVC$Year = rep(seq(2018, 2020, 1), each = dim(hex_grid)[1])
match.set = oni.SVC %>%
dplyr::select(Year, Region, Mean)
hex.2020 = hex.2019 = hex.2018 = hex.fort
hex.2020$Year = 2020
hex.2019$Year = 2019
hex.2018$Year = 2018
hex.yrs = rbind(hex.2020, hex.2019, hex.2018)
hex.yrs = left_join(hex.yrs, match.set, by = c("Year", "Region"))
ggplot(hex.yrs, aes(long, lat, group = group, fill = Mean)) + geom_polygon(col = "gray85") +
scale_fill_viridis(name = "Estimated GH Density", discrete = F, option = "turbo",
direction = -1, alpha = 0.75, na.value = "white") + geom_polygon(data = West.fort,
aes(long, lat, group = group), fill = "transparent", col = "gray65", size = 1) +
xlab(" ") + ylab(" ") + facet_wrap(~Year, ncol = 3) + 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(), legend.position = "bottom", strip.text = element_text(size = 25,
face = "bold"), strip.background = element_blank(), legend.key = element_blank(),
legend.key.size = unit(1, "line"), legend.key.width = unit(3, "line"), legend.text = element_text(size = 12,
face = "bold"), legend.title = element_text(size = 18, face = "bold"), axis.title.x = element_blank(),
axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(),
axis.ticks.x = element_blank(), axis.ticks.y = element_blank(), plot.title = element_blank()) +
guides(fill = guide_colourbar(title.position = "top", title.hjust = 0.5))
# rmdtemplates::line_cite(pkgs) # This creates a single line citing all
# packages
rmdtemplates::list_cite(pkgs) # This creates a 'thightlist' of all packages 11 Used Packages
We used the following packages to create this document:
- Package: INLA by Havard Rue, Lindgren, and Teixeira Krainski (2021)
- Package: spdep by Bivand (2022)
- Package: dagitty by Textor, van der Zander, and Ankan (2021)
- Package: ggdag by Barrett (2021)
- Package: knitr by Xie (2022)
- Package: tidyverse by Wickham (2021)
- Package: dplyr by Wickham et al. (2022)
- Package: rmdformats by Barnier (2021)
- Package: kableExtra by Zhu (2021)
- Package: scales by Wickham and Seidel (2022)
- Package: psych by Revelle (2022)
- Package: rmdtemplates by Calero Valdez (2020)
- Package: viridis by Garnier (2021)
- Package: DT by Xie, Cheng, and Tan (2022)
- Package: ggplot2 by Wickham et al. (2021)