inst/doc/spotr-intro.R

In spotr: Estimate Spatial Population Indices from Ecological Abundance Data

## ----setup, message=FALSE-----------------------------------------------------
library(spotr)
library(mgcv)
library(brms)
library(sf)
library(ggplot2)
library(dplyr)

## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  fig.width = 5
)

## -----------------------------------------------------------------------------
data(cuckoo)
head(cuckoo)

## ----fig.height = 7-----------------------------------------------------------
cuckoo = st_as_sf(cuckoo, coords = c("lon", "lat"), crs = st_crs(4326), remove = FALSE)
data(swe_map)

## ----fig.height = 7-----------------------------------------------------------
ggplot(swe_map) +  geom_sf() + geom_sf(aes(col = count), data = subset(cuckoo, yr == 2010)) 

## -----------------------------------------------------------------------------
gam_fit = gam(count ~ te(yr, lat, lon, d = c(1,2), k = c(8,30)), data = cuckoo, family = nb())

## ----fig.height = 7-----------------------------------------------------------
spat_grid = st_make_grid(cuckoo, cellsize = .5, what = "centers") |> st_sf() |>
  st_intersection(swe_map)
# Add columns with longitude and latitude
spat_grid[, c("lon", "lat")] = st_coordinates(spat_grid)
ggplot(swe_map) + geom_sf() + geom_sf(data = spat_grid)

## -----------------------------------------------------------------------------
pred_points = merge(data.frame(yr = 2000:2020), spat_grid, by = NULL)

## -----------------------------------------------------------------------------
swe_ind = index(gam_fit, timevar = "yr", newdata = pred_points)
head(swe_ind)

## -----------------------------------------------------------------------------
swe_ind |> 
  ggplot(aes(x = yr, y = index)) + 
  geom_linerange(aes(ymin = ci_2.5, ymax = ci_97.5), col = "gray50") + 
  geom_linerange(aes(ymin = ci_10, ymax = ci_90), col = "red3", lwd =1) +
  geom_point() + ylim(.8,1.2)

## -----------------------------------------------------------------------------
index(gam_fit, timevar = "yr", newdata = pred_points, baseline = 2000:2010) |> 
  ggplot(aes(x = yr, y = index)) +   
  geom_linerange(aes(ymin = ci_2.5, ymax = ci_97.5), col = "gray50") + 
  geom_linerange(aes(ymin = ci_10, ymax = ci_90), col = "red3", lwd =1) +
  geom_point() + ylim(.8,1.2)

## -----------------------------------------------------------------------------
index(gam_fit, timevar = "yr", newdata = pred_points, type = "delta") |> 
  ggplot(aes(x = yr, y = index)) +   
  geom_linerange(aes(ymin = ci_2.5, ymax = ci_97.5), col = "gray50") + 
  geom_linerange(aes(ymin = ci_10, ymax = ci_90), col = "red3", lwd =1) +
  geom_point() 

## ----fig.height = 10, fig.width=7---------------------------------------------
index(gam_fit, timevar = "yr", byvar = "county", newdata = pred_points, baseline = 2000:2010) |> 
  ggplot(aes(x = yr, y = index, ymin = ci_2.5, ymax = ci_97.5)) +   
  geom_linerange(aes(ymin = ci_2.5, ymax = ci_97.5), col = "gray50") + 
  geom_linerange(aes(ymin = ci_10, ymax = ci_90), col = "red3", lwd = 1) +
  geom_point(size=.5) + 
  scale_x_continuous(breaks = c(2005,2015)) + facet_wrap(vars(county))

## ----fig.height = 10, fig.width=7---------------------------------------------
index(gam_fit, timevar = "yr", byvar = "county", newdata = pred_points, baseline = 2000:2010, type = "global") |> 
  ggplot(aes(x = yr, y = index, ymin = ci_2.5, ymax = ci_97.5)) +   
  geom_linerange(aes(ymin = ci_2.5, ymax = ci_97.5), col = "gray50") + 
  geom_linerange(aes(ymin = ci_10, ymax = ci_90), col = "red3", lwd = 1) +
  geom_point(size=.5) + 
  scale_x_continuous(breaks = c(2005,2015)) + facet_wrap(vars(county))

## ----fig.height = 10, fig.width=7---------------------------------------------
# Count number of predictions points per county 
pred_points = pred_points |> group_by(county, yr) |> mutate(ppc = n())
w = 1/pred_points$ppc

# Denominator weights
bw = 1/(nrow(spat_grid) * length(2000:2010))
  
index(gam_fit, timevar = "yr", byvar = "county", newdata = pred_points, weights = w, bweights = bw, baseline = 2000:2010, type = "global")  |> 
  ggplot(aes(x = yr, y = index, ymin = ci_2.5, ymax = ci_97.5)) +   
  geom_linerange(aes(ymin = ci_2.5, ymax = ci_97.5), col = "gray50") + 
  geom_linerange(aes(ymin = ci_10, ymax = ci_90), col = "red3", lwd = 1) +
  geom_point(size=.5) + 
  scale_x_continuous(breaks = c(2005,2015)) + facet_wrap(vars(county))