standard prediction grid
pred_grid <- tibble(biodoy = 1:365)
salt_features <- Bottom %>%
filter(!is.na(bottom_salt), !is.na(biodoy)) %>%
group_by(bioyear, xi, eta) %>%
nest() %>%
mutate(
# fit GAM
model = purrr::map(
data,
~ mgcv::gam(
bottom_salt ~ s(biodoy, bs = "cc", k = 10),
data = .x,
method = "REML",
knots = list(biodoy = c(1, 366))
)
),
# predict onto full seasonal cycle
pred = purrr::map(
model,
~ predict(.x, newdata = pred_grid)
),
# mean salinity (observed)
mean_salt = purrr::map_dbl(data, ~ mean(.x$bottom_salt, na.rm = TRUE)),
# seasonal range (from GAM)
max_salt = purrr::map_dbl(pred, max, na.rm = TRUE),
min_salt = purrr::map_dbl(pred, min, na.rm = TRUE),
range_salt = max_salt - min_salt,
# timing of extremes
doy_max = purrr::map_int(pred, which.max),
doy_min = purrr::map_int(pred, which.min),
# threshold metrics (choose thresholds based on biology)
days_below_32 = purrr::map_int(pred, ~ sum(.x < 32, na.rm = TRUE)),
days_32_to_34 = purrr::map_int(pred, ~ sum(.x >= 32 & .x <= 34, na.rm = TRUE)),
days_above_34 = purrr::map_int(pred, ~ sum(.x > 34, na.rm = TRUE)),
# variability
sd_salt = purrr::map_dbl(pred, sd, na.rm = TRUE),
# model complexity
edf = purrr::map_dbl(model, ~ summary(.x)$s.table[1, "edf"]),
dev_expl = purrr::map_dbl(model, ~ summary(.x)$dev.expl)
) %>%
select(
bioyear, xi, eta,
mean_salt,
max_salt,
min_salt,
range_salt,
doy_max,
doy_min,
days_below_32,
days_32_to_34,
days_above_34,
sd_salt,
edf,
dev_expl
) %>%
ungroup()