openness_occupancy_modeling.Rmd

---
title: "Code from 'A lidar-based openness index to aid conservation planning for grassland wildlife' (occupancy modeling)"
author: "Mike Allen"
date: "11/10/2021"
output: html_document
---
# Load packages and data
R version 4.1.0 (2021-05-18)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19043)
```{r}
library(raster) # raster_3.4-13
library(sf) # sf_1.0-1 
library(tidyverse) # tidyverse_1.3.1
library(lubridate) # lubridate_1.7.10
library(unmarked) # unmarked_1.1.1
library(ubms) # ubms_1.0.2.9007
library(AICcmodavg) # AICcmodavg_2.3-1
library(ncf) # ncf_1.2-9
select <- dplyr::select # resolves namespace conflicts

# read in Duke Farms field borders shapefile
duke <- 
   read_sf("data/duke_GL_pred_areas.shp") %>%
   st_transform(crs = 32618)

# read in Duke Farms treelines shapefile
trees <- st_read("data/duke_tree_lines.shp") %>%
   mutate(Id = 1:10)

# read in bird and covariate data for occupancy modeling
g.occ <- read.csv("data/grsp.occupancy.openness.csv")

# create sf object of Duke Farms survey points
duke_pts <- g.occ %>% 
   dplyr::select(pt, latitude, longitude) %>%
   st_as_sf(coords = c("longitude", "latitude"),
            crs = 4326) %>%
   st_transform(crs = crs(duke))

# add in UTM coordinates
g.occ <- g.occ %>%
   mutate(y = st_coordinates(duke_pts)[,2],
          x = st_coordinates(duke_pts)[,1])

# Create unmarked data.frame for occupancy analyses
# Note: ordinal date and time are divided by constants to facilitate model convergence
g.umf <- unmarkedFrameOccu(
  y = as.matrix(g.occ[, 5:8]),
  # detection/non-detection measurements
  siteCovs = g.occ[, c(2, 21:34)],
  # site-level covariates
  obsCovs = list(
    # observation-specific covariates
    ord = as.matrix(g.occ[, 13:16]) / 1000,
    # ordinal date
    dtime = as.matrix(g.occ[, 17:20]) / 10,
    # time of day (decimal hours)
    obs = as.matrix(g.occ[, 9:12])
  )
) # observer        
```
# Calculate openness summary stats for Table 2
```{r}
sitecov_sums <-
  rbind.data.frame(
    aggregate(g.occ[, 21:32], by = list(g.occ$site), mean),
    aggregate(g.occ[, 21:32], by = list(g.occ$site), min),
    aggregate(g.occ[, 21:32], by = list(g.occ$site), max),
    t(data.frame(c(
      Group.1 = "All", apply(g.occ[, 21:32], 2, mean)
    ))),
    t(data.frame(c(
      Group.1 = "All", apply(g.occ[, 21:32], 2, min)
    ))),
    t(data.frame(c(
      Group.1 = "All", apply(g.occ[, 21:32], 2, max)
    )))
  ) %>%
  mutate(metric = c("mean", "mean", "min", "min", "max",
                    "max", "mean", "min", "max"))

row.names(sitecov_sums) <- NULL
sitecov_sums
```
# Evaluate all Grasshopper Sparrow detection sub-model structures 
Evaluate all combinations of the 3 covariates based on AICc. Starting values are provided for some models to aid convergence.
```{r}
p._psi. <-
  occu(~ 1 ~ 1,
       data = g.umf)

p.ord_psi. <-
  occu(~ ord ~ 1,
       data = g.umf)

p.ord2_psi. <-
  occu(~ ord + I(ord) ^ 2 ~ 1,
       data = g.umf,
       starts = c(-2,-4,-4, 2))

p.dtime_psi. <-
  occu(~ dtime ~ 1,
       data = g.umf)

p.dtime2_psi. <-
  occu(~ dtime + I(dtime) ^ 2 ~ 1,
       data = g.umf,
       starts = c(0, 0,-4, 2))

p.ord.dtime_psi. <-
  occu(~ ord + dtime ~ 1,
       data = g.umf)

p.ord.dtime2_psi. <-
  occu(~ ord + dtime + I(dtime) ^ 2 ~ 1,
       data = g.umf,
       starts = c(0,-1, 5,-3, 3))

p.ord2.dtime_psi. <-
  occu(~ ord + I(ord) ^ 2 + dtime ~ 1,
       data = g.umf,
       starts = c(0,-1, 5,-3, 3))

p.ord2.dtime2_psi. <-
  occu(
    ~ ord + I(ord) ^ 2 + dtime + I(dtime) ^ 2 ~ 1,
    data = g.umf,
    starts = c(0,-1, 5,-3, 3, 0)
  )

# obs

p.obs_psi. <-
  occu(~ obs ~ 1,
       data = g.umf)

p.ord.obs_psi. <-
  occu(~ ord + obs ~ 1,
       data = g.umf)

p.ord2.obs_psi. <-
  occu(~ obs + ord + I(ord) ^ 2 ~ 1,
       data = g.umf,
       starts = c(-2,-4,-4, 2,-1,-2))

p.dtime.obs_psi. <-
  occu(~ dtime + obs ~ 1,
       data = g.umf)

p.dtime2.obs_psi. <-
  occu(~ obs + dtime + I(dtime) ^ 2 ~ 1,
       data = g.umf,
       starts = c(0,-1,-2, 2,-4, 2))

p.ord.dtime.obs_psi. <-
  occu(~ ord + dtime + obs ~ 1,
       data = g.umf)

p.ord.dtime2.obs_psi. <-
  occu(~ obs + ord + dtime + I(dtime) ^ 2 ~ 1,
       data = g.umf,
       starts = c(0,-1, 1, 1, 4,-3, 2))

p.ord2.dtime.obs_psi. <-
  occu(~ obs + ord + I(ord) ^ 2 + dtime ~ 1,
       data = g.umf,
       starts = c(0,-1, 1, 1, 4,-1, 2))

p.ord2.dtime2.obs_psi. <-
  occu(
    ~ obs + ord + I(ord) ^ 2 + dtime + I(dtime) ^ 2 ~ 1,
    data = g.umf,
    starts = c(0, 1, 1, 4,-1,-2,-4,-2)
  )

  # create AICc model selction table
AICcmodavg::aictab(
  list(
    p._psi. = p._psi.,
    p.ord_psi. = p.ord_psi.,
    p.ord2_psi. = p.ord2_psi.,
    p.dtime_psi. = p.dtime_psi.,
    p.dtime2_psi. = p.dtime2_psi.,
    p.ord.dtime_psi. = p.ord.dtime_psi.,
    p.ord2.dtime_psi. = p.ord2.dtime_psi.,
    p.ord.dtime2_psi. = p.ord.dtime2_psi.,
    p.ord2.dtime2_psi. = p.ord2.dtime2_psi.,
    p.obs_psi. = p.obs_psi.,
    p.ord.obs_psi. = p.ord.obs_psi.,
    p.ord2.obs_psi. = p.ord2.obs_psi.,
    p.dtime.obs_psi. = p.dtime.obs_psi.,
    p.dtime2.obs_psi. = p.dtime2.obs_psi.,
    p.ord.dtime.obs_psi. = p.ord.dtime.obs_psi.,
    p.ord2.dtime.obs_psi. = p.ord2.dtime.obs_psi.,
    p.ord.dtime2.obs_psi. = p.ord.dtime2.obs_psi.,
    p.ord2.dtime2.obs_psi. = p.ord2.dtime2.obs_psi.
  )
)
```
# Evaluate mean vs. max openness models, with vs. without wires
Final model list consisting of the four top detection sub-models combined with the four occupancy-submodel variables of interest:  
mean2018 = mean-angle openness (computed with powerlines present)
max2018 = maximum-angle openness (computed with powerlines present)
No_wires_mean = mean-angle openness (computed with powerlines digitally erased)
No_wires_max = maximum-angle openness (computed with powerlines digitally erased)
```{r}
p.obs_psi. <-
  occu(~ obs ~ site,
       data = g.umf)

p.obs_psi.mean <-
  occu(~ obs ~ site + mean2018,
       data = g.umf)

p.obs_psi.max <-
  occu(~ obs ~ site + max2018,
       data = g.umf)

p.obs_psi.No_wires_mean <-
  occu(~ obs ~ site + No_wires_mean,
       data = g.umf)

p.obs_psi.No_wires_max <-
  occu(~ obs ~ site + No_wires_max,
       data = g.umf)

p.ord.obs_psi. <-
  occu(~ ord + obs ~ site,
       data = g.umf)

p.ord.obs_psi.mean <-
  occu(~ ord + obs ~ site + mean2018,
       data = g.umf)

p.ord.obs_psi.max <-
  occu(~ ord + obs ~ site + max2018,
       data = g.umf)

p.ord.obs_psi.No_wires_mean <-
  occu(~ ord + obs ~ site + No_wires_mean,
       data = g.umf)

p.ord.obs_psi.No_wires_max <-
  occu(~ ord + obs ~ site + No_wires_max,
       data = g.umf)

p.dtime.obs_psi. <-
  occu(~ dtime + obs ~ site,
       data = g.umf)

p.dtime.obs_psi.mean <-
  occu(~ dtime + obs ~ site + mean2018,
       data = g.umf)

p.dtime.obs_psi.max <-
  occu(~ dtime + obs ~ site + max2018,
       data = g.umf)

p.dtime.obs_psi.No_wires_mean <-
  occu(~ dtime + obs ~ site + No_wires_mean,
       data = g.umf)

p.dtime.obs_psi.No_wires_max <-
  occu(~ dtime + obs ~ site + No_wires_max,
       data = g.umf)

p.ord.dtime_psi. <-
  occu(~ ord + dtime + obs ~ site,
       data = g.umf)

p.ord.dtime.obs_psi. <-
  occu(~ ord + dtime + obs ~ 1,
       data = g.umf)

p.ord.dtime.obs_psi.mean <-
  occu(~ ord + dtime + obs ~ site + mean2018,
       data = g.umf)

p.ord.dtime.obs_psi.max <-
  occu(~ ord + dtime + obs ~ site + max2018,
       data = g.umf)

p.ord.dtime.obs_psi.No_wires_mean <-
  occu(~ ord + dtime + obs ~ site + No_wires_mean,
       data = g.umf)

p.ord.dtime.obs_psi.No_wires_max <-
  occu(~ ord + dtime + obs ~ site + No_wires_max,
       data = g.umf)

  # create AICc model selection table
(mod_table <-   AICcmodavg::aictab(
  list(
    p.obs_psi. = p.obs_psi.,
    p.obs_psi.mean = p.obs_psi.mean,
    p.obs_psi.max = p.obs_psi.max,
    p.obs_psi.No_wires_mean = p.obs_psi.No_wires_mean,
    p.obs_psi.No_wires_max = p.obs_psi.No_wires_max,
    
    p.ord.obs_psi. = p.ord.obs_psi.,
    p.ord.obs_psi.mean = p.ord.obs_psi.mean,
    p.ord.obs_psi.max = p.ord.obs_psi.max,
    p.ord.obs_psi.No_wires_mean = p.ord.obs_psi.No_wires_mean,
    p.ord.obs_psi.No_wires_max = p.ord.obs_psi.No_wires_max,
    
    p.dtime.obs_psi. = p.dtime.obs_psi.,
    p.dtime.obs_psi.mean = p.dtime.obs_psi.mean,
    p.dtime.obs_psi.max = p.dtime.obs_psi.max,
    p.dtime.obs_psi.No_wires_mean = p.dtime.obs_psi.No_wires_mean,
    p.dtime.obs_psi.No_wires_max = p.dtime.obs_psi.No_wires_max,
    
    p.ord.dtime.obs_psi. = p.ord.dtime.obs_psi.,
    p.ord.dtime.obs_psi.mean = p.ord.dtime.obs_psi.mean,
    p.ord.dtime.obs_psi.max = p.ord.dtime.obs_psi.max,
    p.ord.dtime.obs_psi.No_wires_mean = p.ord.dtime.obs_psi.No_wires_mean,
    p.ord.dtime.obs_psi.No_wires_max = p.ord.dtime.obs_psi.No_wires_max
  )
))

# create AICc model selection table for only models containing 'max openness'
AICcmodavg::aictab(
  list(
    p.obs_psi.max = p.obs_psi.max,
    p.obs_psi.No_wires_max = p.obs_psi.No_wires_max,
    
    p.ord.obs_psi.max = p.ord.obs_psi.max,
    p.ord.obs_psi.No_wires_max = p.ord.obs_psi.No_wires_max,
    
    p.dtime.obs_psi.max = p.dtime.obs_psi.max,
    p.dtime.obs_psi.No_wires_max = p.dtime.obs_psi.No_wires_max,
    
    p.ord.dtime.obs_psi.max = p.ord.dtime.obs_psi.max,
    p.ord.dtime.obs_psi.No_wires_max = p.ord.dtime.obs_psi.No_wires_max
  )
)    
```
# Calculate & evaluate distance to edge
Calculate distance from survey points to the nearest forest edge. Evaluate a model containing this covariate compared with the model set containing the openness covariates.

Distance code adapted from a gis.stackexchange.com post (see References). 
```{r}
# get index value for nearest treeline segment to each survey point
nearest = st_nearest_feature(duke_pts, trees)

# calculate distance to nearest treeline for each survey point
dist = st_distance(duke_pts, trees[nearest, ], by_element = TRUE)

# attach index values to sf object of points
pt_tree_dist <-
  cbind(duke_pts, st_drop_geometry(trees)[nearest, ]) %>%
  rename(treelineID = 2)

# attach distance values to sf object of points
pt_tree_dist$dist <- dist

# make distance-to-nearest-treeline object into a data.frame
pt_tree_df <- pt_tree_dist %>%
  st_drop_geometry() %>%
  as.data.frame() %>%
  mutate(dist = as.numeric(dist)) %>%
  dplyr::select(-treelineID)

# add distance to edge to site covariates
site.cov_dist <- g.occ[, c(1, 2, 21:32)] %>%
  left_join(pt_tree_df, by = "pt") %>%
  mutate(max2018 = max2018,
         dist = dist)

# Create unmarked data.frame for occupancy analyses
# Note: ordinal date and time are divided by constants to facilitate model convergence
g.umf_dist <- unmarkedFrameOccu(
  y = as.matrix(g.occ[, 5:8]),
  # detection/non-detection observations
  siteCovs = site.cov_dist,
  # site-level covariates
  obsCovs = list(
    # observer-specific covariates
    ord = as.matrix(g.occ[, 13:16]) / 1000,
    # ordinal date
    dtime = as.matrix(g.occ[, 17:20]) / 10,
    # time of day (decimal hours)
    obs = as.matrix(g.occ[, 9:12])
  )
) # observer    

# Run top openness model and distance to edge model
p.obs_psi.max <-
  occu(~ obs ~ site + max2018,
       data = g.umf_dist)

p.obs_psi.dist <-
  occu(~ obs ~ site + dist,
       data = g.umf_dist)

# coefficients
 coef(p.obs_psi.dist)
 confint(p.obs_psi.dist, type = "state")
 
 # Model selection table including the 'distance to edge' model
 AICcmodavg::aictab(
   list(
     p.obs_psi. = p.obs_psi.,
     p.obs_psi.mean = p.obs_psi.mean,
     p.obs_psi.max = p.obs_psi.max,
     p.obs_psi.No_wires_mean = p.obs_psi.No_wires_mean,
     p.obs_psi.No_wires_max = p.obs_psi.No_wires_max,
     
     p.ord.obs_psi. = p.ord.obs_psi.,
     p.ord.obs_psi.mean = p.ord.obs_psi.mean,
     p.ord.obs_psi.max = p.ord.obs_psi.max,
     p.ord.obs_psi.No_wires_mean = p.ord.obs_psi.No_wires_mean,
     p.ord.obs_psi.No_wires_max = p.ord.obs_psi.No_wires_max,
     
     p.dtime.obs_psi. = p.dtime.obs_psi.,
     p.dtime.obs_psi.mean = p.dtime.obs_psi.mean,
     p.dtime.obs_psi.max = p.dtime.obs_psi.max,
     p.dtime.obs_psi.No_wires_mean = p.dtime.obs_psi.No_wires_mean,
     p.dtime.obs_psi.No_wires_max = p.dtime.obs_psi.No_wires_max,
     
     p.ord.dtime.obs_psi. = p.ord.dtime.obs_psi.,
     p.ord.dtime.obs_psi.mean = p.ord.dtime.obs_psi.mean,
     p.ord.dtime.obs_psi.max = p.ord.dtime.obs_psi.max,
     p.ord.dtime.obs_psi.No_wires_mean = p.ord.dtime.obs_psi.No_wires_mean,
     p.ord.dtime.obs_psi.No_wires_max = p.ord.dtime.obs_psi.No_wires_max,
     p.obs_psi.dist = p.obs_psi.dist
   )
 )  

```
# Re-run top model as spatial occupancy model in package ubms 
Code for this restricted spatial regression (RSR) occupancy model was adapted from code in kenkellner.com (see References).
```{r}
# view how threshold value affects the spatial neighborhood definition
with(g.occ[,c("site", "max2018", "x", "y")], 
     RSR(x, y, threshold=150, plot_site=215))

# Define the model
form_RSR150 <- ~obs ~site + I(max2018/100) + RSR(x, y, threshold=150)

# Fit the model in STAN
# NOTE: takes several minutes
options(mc.cores=3)
fit_ubms_RSR150 <- stan_occu(form_RSR150, g.umf, 
                             chains=3, iter = 15000, cores = 3)

# Examine model output
fit_ubms_RSR150

# save as an rds file to save time
# saveRDS(fit_ubms_RSR150, "output/fit_ubms_RSR150.rds")
fit_ubms_RSR150 <- readRDS("output/fit_ubms_RSR150.rds")
```
# Detection probability
Print model-estimated detection probability by observer,
```{r}
predict(fit_ubms_RSR150,
        submodel = "det",
        newdata = data.frame(obs = c(
          "Mike Allen", "Charles Barreca", "Thom Almendinger"
        )))
```
# Create Fig. 2 - Grasshopper sparrow: Graph relationship with max 
```{r}
# Format model predicted occupancy for plotting
psi_plot_data <- predict(
  fit_ubms_RSR150,
  newdata = expand.grid(
    site = c("K", "S"),
    max2018 = seq(12.9, 66.6, length.out = 100)
  ),
  "state",
  re.form = NA
) %>%
  rename(q2.5 = 3, q97.5 = 4) %>%
  bind_cols(expand.grid(
    site = c("K", "S"),
    max2018 = seq(12.9, 66.6, length.out = 100)
  ))

# Format raw detection/non-detection data for rug plot
raw_data <- g.occ %>%
  select(1, 2, 5:8, max2018) %>%
  mutate(detect = apply(g.umf@y, 1, max) * 100)

ggplot(psi_plot_data) +
  geom_ribbon(aes(
    x = max2018,
    ymin = q2.5 * 100,
    ymax = q97.5 * 100,
    fill = site
  ),
  alpha = 0.5) +
  geom_line(aes(x = max2018, y = Predicted * 100, group = site), size = 1.5) +
  scale_fill_manual(values = c("firebrick", "steelblue")) +
  theme_bw() +
  theme(axis.text = element_text(size = 14),
        axis.title = element_text(size = 14)) +
  labs(x = "Openness index (90 - max degrees to horizon)",
       y = "Grasshopper Sparrow Occupancy (%)", fill = "Site") +
  geom_rug(
    aes(x = max2018, y = detect, color = site),
    data = filter(raw_data, detect == 0),
    sides = "b"
  ) +
  geom_rug(
    aes(x = max2018, y = detect, color = site),
    data = filter(raw_data, detect == 100),
    sides = "t"
  ) +
  scale_color_manual(values = c("firebrick", "steelblue")) +
  guides(color = "none") +
  scale_x_continuous(breaks = seq(20, 60, by = 10))
#ggsave("figures/Fig_2_relationship_max_openness_RSR150.jpg", height = 5, width = 7, dpi = 600)
```
# Change in predicted occupancy in various openness scenarios
Predicted change in patch-level occupancy based on UBMS model. Predictions and credible intervals and are generated by resampling from the parameter posterior distributions.
```{r}
# load in the pre-fit model to save time (or run it above)
fit_ubms_RSR150 <- readRDS("output/fit_ubms_RSR150.rds")

# get posteriors for fixed parameters of top model
topmod <- cbind(
  extract(fit_ubms_RSR150, "beta_state[(Intercept)]") %>%
    do.call(c, .),
  extract(fit_ubms_RSR150, "beta_state[siteS]") %>%
    do.call(c, .),
  extract(fit_ubms_RSR150,
          "beta_state[I(max2018/100)]") %>%
    do.call(c, .)
)
row.names(topmod) <- NULL

# predicted occupancy (resampling posteriors of the paraemter estimates)
# number of samples to draw
ni = 7500 

# create a randomly-ordered index
random_index <- sample(1:ni, replace = T)

# draw ni samples of each parameter estimate from the posterior
int <- topmod[random_index,1]
beta.site <- topmod[random_index,2]
beta.open <- topmod[random_index,3]

# subset occupancy data by field
g.occ.S <- g.occ %>% filter(site == "S")
g.occ.K <- g.occ %>% filter(site == "K")

## Calculate % change in number of occupied patches from 'no action' scenario

# Skeet - no wires
no_action_S <- 
  lapply(1:ni,
         function(x){
           sum(
             plogis(
             int[x] + beta.site[x] + 
             (g.occ.S$max2018/100)*beta.open[x]
           )
           )
         }
  ) %>%
  do.call(c, .)

no_wires_S <- 
  lapply(1:ni,
         function(x){
           sum(
             plogis(
             int[x] + beta.site[x] + 
             (g.occ.S$No_wires_max/100)*beta.open[x]
           )
           )
         }
  ) %>%
  do.call(c, .)

no_wires <- 100*(no_wires_S - no_action_S) / no_action_S

# Kaufman - no action
no_action_K <- 
  lapply(1:ni,
         function(x){
           sum(
             plogis(
             int[x] + 
             (g.occ.K$max2018/100)*beta.open[x]
           )
           )
         }
  ) %>%
  do.call(c, .)

# no N tree line in Kaufman
no_N_K <- 
  lapply(1:ni,
         function(x){
           sum(
             plogis(
             int[x] + 
             (g.occ.K$No_N_max/100)*beta.open[x]
           )
           )
         }
  ) %>%
  do.call(c, .)

no_N <- 100*(no_N_K - no_action_K) / no_action_K

# no SW tree line in Kaufman
no_SW_K <- 
  lapply(1:ni,
         function(x){
           sum(
             plogis(
             int[x] + 
             (g.occ.K$No_SW_max/100)*beta.open[x]
           )
           )
         }
  ) %>%
  do.call(c, .)

no_SW <- 100*(no_SW_K - no_action_K) / no_action_K

# no SE tree line in Kaufman
no_SE_K <- 
  lapply(1:ni,
         function(x){
           sum(
             plogis(
             int[x] + 
             (g.occ.K$No_SE_max/100)*beta.open[x]
           )
           )
         }
  ) %>%
  do.call(c, .)

no_SE <- 100*(no_SE_K - no_action_K) / no_action_K

# no SE tree line in Kaufman
no_trees_K <- 
  lapply(1:ni,
         function(x){
           sum(
             plogis(
             int[x] + 
             (g.occ.K$No_trees_max/100)*beta.open[x]
           )
           )
         }
  ) %>%
  do.call(c, .)

no_trees <- 100*(no_trees_K - no_action_K) / no_action_K

# Collect all predicted occupancy values into one data frame
scenario_preds <- cbind(no_SW, no_SE, no_N, 
                        no_trees, no_wires) %>%
  apply(., 2, 
        function(x)quantile(x,c(0.025, .1, 0.5, .9, 0.975))) %>%
  t() %>%
  as.data.frame() %>%
  rename(q2.5 = 1, q10 = 2, med = 3, q90 = 4, q97.5 = 5) %>%
  mutate(field = c(rep("K", 4), "S"),
   action = c("Remove SW tree line", "Remove SE tree line", 
              "Remove N tree line", "Remove all tree lines", 
              "Remove power lines")
   ) %>%
   mutate(action = fct_reorder(.f = action, .x = c(1, 2, 3, 4, 5), 
                               .desc = T),
          type = c("No SW", "No SE", "No N", 
                   "No trees", "No wires"))
```
# Create Fig. 3 - predicted change in occupancy in different management scenarios
```{r}
ggplot(scenario_preds) +
  geom_errorbar(
    aes(
      xmin = q2.5,
      xmax = q97.5,
      y = action,
      color = field
    ),
    width = 0,
    size = 1.25
  ) +
  geom_errorbar(aes(
    xmin = q10,
    xmax = q90,
    y = action,
    color = field
  ),
  width = 0,
  size = 2.25) +
  geom_point(aes(x = med, y = action, color = field),
             size = 4.75,
             pch = 16) +
  scale_color_manual(values = c("firebrick", "steelblue")) +
  theme_bw() +
  theme(axis.text = element_text(size = 14),
        axis.title = element_text(size = 14)) +
  labs(y = "", x = "% change in Grasshopper Sparrow occupancy",
       color = "Site") +
  xlim(0, 20.5)

# ggsave(
#   "figures/Fig_3_management_scenarios_RSR150.jpg",
#   height = 6,
#   width = 7,
#   dpi = 600
# )
```
# Appendix A. Testing for spatial autocorrelation
Bjornstad & Bjornstad (2016)
```{r}
# Define the models: non-RSR model and RSR models with 100m threshold
# RSR model with 150 m threshold was already fit above

form_noRSR <- ~obs ~site + I(max2018/100)
form_RSR100 <- ~obs ~site + I(max2018/100) + RSR(x, y, threshold=100)

# Fit the models in STAN

# no RSR model
options(mc.cores=3)
fit_ubms_noRSR <- stan_occu(form_noRSR, g.umf, chains=3, 
                            iter = 15000, cores = 3) 

# examin model output
fit_ubms_noRSR

# load from rds file to save time
# saveRDS(fit_ubms_noRSR, "output/fit_ubms_noRSR.rds")
fit_ubms_noRSR <- readRDS("output/fit_ubms_noRSR.rds")

# RSR 100 model
fit_ubms_RSR100 <- stan_occu(form_RSR100, g.umf, chains=3, 
                             iter = 15000, cores = 3)

# examine model output
fit_ubms_RSR100

# load from rds file to save time
# saveRDS(fit_ubms_RSR100, "output/fit_ubms_RSR100.rds")
fit_ubms_RSR100 <- readRDS("output/fit_ubms_RSR100.rds")

# load RSR 150 m model from rds file (or run it above)
# RSR 150 model
fit_ubms_RSR150 <- readRDS("output/fit_ubms_RSR150.rds")

# extract residuals
resids_noRSR <- residuals(fit_ubms_noRSR, submodel = "state") %>%
  apply(., 2, median) 

resids_RSR100 <- residuals(fit_ubms_RSR100, submodel = "state") %>%
  apply(., 2, median) 

resids_RSR150 <- residuals(fit_ubms_RSR150, submodel = "state") %>%
  apply(., 2, median) 

# combine residuals from the 3 models into one data frame
rsd <- cbind.data.frame(resids_noRSR,
                        resids_RSR100,
                        resids_RSR150) %>%
  mutate(det = apply(g.umf@y, 1, max),
         x = g.occ$x,
         y = g.occ$y) 

# perform spatial autocorrelation test (no RSR)
sac_model_noRSR <- 
  ncf::correlog(x = rsd$x, 
                y = rsd$y, 
                z = rsd$resids_noRSR, 
                increment = 100, 
                  resamp = 3000, latlon = F)

# load output from rds file to save time
# saveRDS(sac_model_noRSR, "output/sac_model_noRSR.rsd")
sac_model_noRSR <- readRDS("output/sac_model_noRSR.rsd")

# perform spatial autocorrelation test (RSR100)
sac_model_RSR100 <- 
  ncf::correlog(x = rsd$x, 
                y = rsd$y, 
                z = rsd$resids_RSR100, 
                increment = 100, 
                  resamp = 3000, latlon = F)

# load from rds file to save time
# saveRDS(sac_model_RSR100, "output/sac_model_RSR100.rsd")
sac_model_RSR100 <- readRDS("output/sac_model_RSR100.rsd")

# perform spatial autocorrelation test (RSR150)
sac_model_RSR150 <- 
  ncf::correlog(x = rsd$x, 
                y = rsd$y, 
                z = rsd$resids_RSR150, 
                increment = 100, 
                  resamp = 3000, latlon = F)

# load output from rds file to save time
# saveRDS(sac_model_RSR150, "output/sac_model_RSR150.rsd")
sac_model_RSR150 <- readRDS("output/sac_model_RSR150.rsd")

# format autocorrelation data for plotting
plot_data_noRSR <- data.frame(
   MoransI = sac_model_noRSR$correlation, 
   lag = sac_model_noRSR$mean.of.class,
   n = sac_model_noRSR$n)

plot_data_RSR100 <- data.frame(
   MoransI = sac_model_RSR100$correlation, 
   lag = sac_model_RSR100$mean.of.class,
   n = sac_model_RSR100$n)

plot_data_RSR150 <- data.frame(
   MoransI = sac_model_RSR150$correlation, 
   lag = sac_model_RSR150$mean.of.class,
   n = sac_model_RSR150$n)

# make a function to plot autocorrelation results
plot_sac <- function(plotdata, sac_anno){
plotdata %>%
   filter(lag < 1600) %>%
ggplot() +
   geom_abline(aes(slope = 0, intercept = 0), linetype = 2, color = "red") +
   geom_line(aes(x = lag, y = MoransI)) +
   geom_point(aes(x = lag, y = MoransI), size = 3) +
   geom_text(aes(x = lag, y = MoransI-0.025, label = as.character(n)), size = 2) +
    annotate(geom = "text", x = 1600, y = 0.3, 
             label = sac_anno, hjust = 1) +
   ylim(-0.31, 0.31) +
    xlim(0,1600) +
   labs(y = "Moran's I", x = "Lag distance (m)") +
   theme_bw() +
   theme(text = element_text(size = 15))
}

# save each plot
(sac_plot_noRSR <- plot_sac(plot_data_noRSR, "Non-spatial model"))
(sac_plot_RSR100 <- plot_sac(plot_data_RSR100, 
                             "RSR model (100 m threshold)"))
(sac_plot_RSR150 <- plot_sac(plot_data_RSR150, 
                             "RSR model (150 m threshold)"))

# create final stacked 3-panel plot
library(patchwork) # patchwork_1.1.0.9000
(sac_plots <- sac_plot_noRSR / sac_plot_RSR100 / sac_plot_RSR150)

ggsave("figures/FigA1 - residual spatial autocorr - RSR.emf",
       plot = sac_plots,
       dpi = 600,
       height = 9, width = 6)

```

# References

Code to calculate distance from survey points to nearest treeline adapted from: https://gis.stackexchange.com/questions/349955/getting-a-new-column-with-distance-to-the-nearest-feature-in-r

Code to run spatial (RSR) occupancy model and extract residuals adapted from:
https://kenkellner.com/blog/ubms-spatial.html
https://kenkellner.com/ubms/reference/residuals-ubmsFit-method.html

Spatial subsampling code ('pruning') modified from here:
https://www.jla-data.net/eng/creating-and-pruning-random-points-and-polygons/