In pbs-assess/gfranges: Groundfish Ranges
Fit spatial temporal models to bottom temperature
library(tidyverse)
library(sdmTMB)
# library(future)
library(dplyr)
library(gfranges)
plot_map_raster <- function(dat, column = "est") {
ggplot(dat, aes_string("X", "Y", fill = column)) +
geom_raster() +
coord_fixed()+
gfplot::theme_pbs() + theme(
axis.title = element_blank(),
axis.ticks = element_blank(),
axis.text = element_blank()
)
}
Model diagnostics function
get_diag <- function(m, response = "temperature_c", variable = "depth_scaled", colour_var = "depth", start_year = 2007, ssids = c(1,3,4,16)) {
predictions <- predict(m)
predictions$residuals <- residuals(m)
predictions <- predictions %>% filter (year >= start_year)
plot_map <- function(dat, column = "est") {
ggplot(dat, aes_string("X", "Y", colour = column)) +
geom_point() +
coord_fixed()+
gfplot::theme_pbs()
}
g <- plot_map(predictions, "est") +
scale_colour_viridis_c() +
ggtitle("Prediction (fixed effects + all random effects)")
print(g)
g <- plot_map(predictions, "est_non_rf") +
ggtitle("Prediction (fixed effects only)") +
scale_colour_viridis_c()
print(g)
g <- plot_map(predictions, "est_rf") +
ggtitle("All random effects only") +
scale_colour_gradient2()
print(g)
g <- plot_map(predictions, "omega_s") +
ggtitle("Spatial random effects only") +
scale_colour_gradient2()
print(g)
g <- plot_map(predictions, "epsilon_st") +
ggtitle("Spatiotemporal random effects only") +
facet_wrap(~year) +
scale_colour_gradient2()
print(g)
g <- ggplot(predictions, aes_string("est", response)) +
geom_point(alpha = 0.2) +
facet_wrap(~year) +
coord_fixed() +
geom_abline()+
gfplot::theme_pbs()
print(g)
g <- ggplot(predictions, aes(est, residuals)) +
geom_point(alpha = 0.2) +
geom_smooth()+
gfplot::theme_pbs()
print(g)
g <- ggplot(predictions, aes_string(variable, "residuals", colour = colour_var)) +
geom_point(alpha = 0.4, size = 1.2) +
geom_smooth(colour="grey", size = 1.2) +
scale_colour_viridis_c(option = "B", direction = -1, begin = 0, end = 0.7,
limits= c(min(predictions[colour_var]), max(predictions[colour_var]))) + #, trans= sqrt
# facet_wrap(~year, scales = "free_x")+
gfplot::theme_pbs() + theme(
axis.text.y = element_text(size = 14),
axis.text.x = element_blank(), axis.ticks.x = element_blank())
print(g)
g <- ggplot(predictions, aes_string(variable, "residuals", colour = colour_var)) +
geom_point(alpha = 0.4, size = 0.7) +
geom_smooth(colour="grey", size = 0.7) +
scale_colour_viridis_c(option = "B", direction = -1, begin = 0, end = 0.7,
limits= c(min(predictions[colour_var]), max(predictions[colour_var]))) + #, trans= sqrt
facet_wrap(~year, scales = "free_x")+
gfplot::theme_pbs() + theme(axis.text.x = element_blank(), axis.ticks.x = element_blank())
print(g)
g <- ggplot(filter(predictions, ssid %in% ssids), aes_string(variable, "residuals", colour = colour_var)) +
geom_point(alpha = 0.4, size = 1) +
geom_smooth(colour="grey", size = 1) +
scale_colour_viridis_c(option = "B", direction = -1, begin = 0, end = 0.7,
limits= c(min(predictions[colour_var]), max(predictions[colour_var]))) + #, trans= sqrt
facet_wrap(~ssid, scales = "free_x")+
gfplot::theme_pbs()
print(g)
g <- ggplot(predictions, aes(X, Y, colour = residuals)) +
geom_point(alpha = 0.2) +
coord_fixed() +
scale_color_gradient2() +
facet_wrap(~year)+
gfplot::theme_pbs()
print(g)
aic <- function(m) {
k <- length(m$model$par)
nll <- m$model$objective
2 * k - 2 * (-nll)
}
print("R^2:")
r2 <- cor(predictions$est, predictions[[response]])^2
print(r2)
print("")
print("AIC:")
print(aic(m))
}
Load data if not using SST
all_sensor <- readRDS(here::here("analysis/tmb-sensor-explore/data/dat-sensor-trawl-processed.rds")) # no SST
all_depth <- all_sensor %>%
# filter(depth_max>10) %>% # maybe we should filter incase missing depth values also indicate a sensor fail?
# filter missing location data and trial year
dplyr::filter(!is.na(latitude), !is.na(longitude)) %>%
# convert lat and lon to UTMs
dplyr::mutate(X = longitude, Y = latitude) %>%
gfplot:::ll2utm(., utm_zone = 9) %>%
# interpolate missing depth values
dplyr::rename(depth = depth_m) %>%
gfplot:::interp_survey_bathymetry()
d_trawl <- all_depth$data %>%
# remove obvious sensor fails, d_trawl %>% filter(depth_max<10), fishing event ids = 481861, 2179089
filter(depth_max>10) %>%
gfplot:::scale_survey_predictors()
d_trawl <- d_trawl[!is.na(d_trawl$depth), ]
d_trawl <- d_trawl[!is.na(d_trawl$temperature_c), ]
d_trawl <- mutate(d_trawl,
DOY = lubridate::yday(date),
DOY_scaled = arm::rescale(lubridate::yday(date)),
DOY_scaled2 = arm::rescale(lubridate::yday(date))^2)
EXLORE TEMPERATURE at DEPTH
ggplot(d_trawl, aes_string("X", "Y", colour = "temperature_c")) +
geom_point() +
facet_wrap(~year) +
coord_fixed() +
scale_color_viridis_c()
ggplot(d_trawl, aes(DOY, temperature_c, colour = -depth, alpha=0.5)) + geom_point() + facet_wrap(~year)
Current temperature model
spde <- sdmTMB::make_mesh(d_trawl, c("X", "Y"), n_knots = 800, type = "kmeans")
plot(spde)
Original model code
# spde <- make_spde(d_trawl$X, d_trawl$Y, n_knots = 800)
# m_temp_rw <- sdmTMB(
# temperature_c ~ 0 + as.factor(year), #+ DOY_scaled + as.factor(ssid)
# time_varying = ~ 0 + depth_scaled + depth_scaled2,
# data = d_trawl,
# ar1_fields = TRUE, # changed AR1 to true
# include_spatial = TRUE,
# time = "year", spde = spde,
# family = gaussian(link = "identity"),
# nlminb_loops = 2,
# newton_steps = 1,
# silent = FALSE
# )
# saveRDS(m_temp_rw, file = (here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-800kn.rds")))
New CV model code
m_temp_cv <- sdmTMB_cv(
temperature_c ~ 0 + as.factor(year), #+ DOY_scaled + as.factor(ssid)
time_varying = ~ 0 + depth_scaled + depth_scaled2,
data = d_trawl,
ar1_fields = TRUE, # changed AR1 to true
include_spatial = TRUE,
time = "year", spde = spde,
family = gaussian(link = "identity"),
k_folds = 5,
# nlminb_loops = 2,
# newton_steps = 1,
silent = FALSE
)
saveRDS(m_temp_cv, file = (here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-800kn-cv.rds")))
m_temp_cv <- readRDS((here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-800kn-cv.rds")))
m_temp_cv$fold_loglik
m_temp_cv$sum_loglik
m_temp_cv$max_gradients
Check predictions against raw data
# m_temp_rw <- readRDS(here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-500kn-2.rds"))
m_temp_rw <- readRDS(here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-800kn.rds"))
predictions <- predict(m_temp_rw)
predictions$residuals <- residuals(m_temp_rw)
cor(predictions$est, predictions$temperature_c)^2
sum(abs(predictions$residuals))
# 2908.879-2789.322
# 119.557/2908.879 # 4% reduction in abs errors
# test against ROMS
events_roms <- readRDS("data/roms_temp_by_events2.rds") %>% select(fishing_event_id, roms, year)
predictions <- left_join(predictions, events_roms)
cor(predictions$roms, predictions$est, use = "pairwise.complete.obs")
cor(predictions$roms, predictions$temperature_c, use = "pairwise.complete.obs")
ggplot(predictions, aes(roms, est, colour = -depth)) + geom_point(alpha=0.3) +
facet_wrap(~year) + coord_cartesian(xlim = c(2,13), ylim = c(2,13)) +
geom_abline(intercept = 0, slope = 1) +
gfplot::theme_pbs()
ggplot(predictions, aes(roms, temperature_c, colour = DOY)) + geom_point(alpha=0.3) +
facet_wrap(~year) + coord_cartesian(xlim = c(2,13), ylim = c(2,13)) +
geom_abline(intercept = 0, slope = 1) +
scale_colour_viridis_c() +
gfplot::theme_pbs()
ggplot(predictions, aes(DOY, (roms-temperature_c), colour = -depth)) + geom_point(alpha=0.3) +
facet_wrap(~year) + #coord_cartesian(xlim = c(2,13), ylim = c(2,13)) +
geom_abline(intercept = 0, slope = 0) +
scale_colour_viridis_c() +
gfplot::theme_pbs()
ggplot(predictions, aes(DOY, (est-roms), colour = -depth)) + geom_point(alpha=0.3) +
facet_wrap(~year) + #coord_cartesian(xlim = c(2,13), ylim = c(2,13)) +
geom_abline(intercept = 0, slope = 0) +
scale_colour_viridis_c() +
gfplot::theme_pbs()
p_nd <- readRDS((here::here("analysis/VOCC/data/predicted-DO-2020-06-20-more2016-roms.rds")))
cor(p_nd$roms_temp, p_nd$temp, use = "pairwise.complete.obs")
ggplot(p_nd, aes(roms_temp, temp, colour = -depth)) + geom_point(alpha=0.04) +
facet_wrap(~year) + coord_cartesian(xlim = c(2,13), ylim = c(2,13)) +
geom_abline(intercept = 0, slope = 1) +
gfplot::theme_pbs()
m_temp_rw <- readRDS((here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-500kn-2.rds")))
get_diag(m_temp_rw, response = "temperature_c", variable = "depth_scaled")
WOULD SST HELP?
Load and prep SST data
all_sensor <- readRDS(here::here("analysis/tmb-sensor-explore/data/dat-sensor-trawl-meanSST-all.rds"))
d_trawl2 <-all_sensor
# all_sensor <- readRDS(here::here("analysis/tmb-sensor-explore/data/dat-sensor-trawl-processed.rds")) # no SST
# glimpse(all_sensor)
#
# all_depth <- all_sensor %>%
# # filter(depth_max>10) %>% # maybe we should filter incase missing depth values also indicate a sensor fail?
# # filter missing location data and trial year
# dplyr::filter(!is.na(latitude), !is.na(longitude)) %>%
#
# # convert lat and lon to UTMs
# dplyr::mutate(X = longitude, Y = latitude) %>%
# gfplot:::ll2utm(., utm_zone = 9) %>%
#
# # interpolate missing depth values
# dplyr::rename(depth = depth_m) %>%
# gfplot:::interp_survey_bathymetry()
# d_trawl <- all_depth$data %>%
# # remove obvious sensor fails, d_trawl %>% filter(depth_max<10), fishing event ids = 481861, 2179089
# filter(depth_max>10) %>%
# gfplot:::scale_survey_predictors()
olddata <- readRDS((here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-400kn.rds")))
d_trawl <- olddata$data
d_trawl <- d_trawl[!is.na(d_trawl$depth), ]
d_trawl <- d_trawl[!is.na(d_trawl$temperature_c), ]
d_trawl2 <- d_trawl[!is.na(d_trawl$SST), ]
d_trawl2 <- d_trawl2[!is.na(d_trawl2$meanSST), ]
d_trawl2 <- mutate(d_trawl2,
DOY = lubridate::yday(date),
DOY_scaled = arm::rescale(DOY),
DOY_scaled2 = DOY_scaled^2,
depth_scaled = arm::rescale(log(depth)),
SST_scaled = arm::rescale(meanSST),
depth_scaled3 = depth_scaled^3,
shallow = ifelse(depth>35,0,1),
deep = ifelse(depth<35,0,1)
)
# d_trawl2[is.na(d_trawl2$shallow), ]
# d_trawl2[is.na(d_trawl2$SST_scaled), ]
Explore SST on survey day
ggplot(d_trawl2, aes_string("X", "Y", colour = "SST")) +
geom_point() +
facet_wrap(~year) +
coord_fixed() +
scale_color_viridis_c()
Is bottom temp related to SST in shallower samples?
ggplot(dplyr::filter(d_trawl2, depth <35), aes(SST,temperature_c, colour = -depth)) + geom_point()
ggplot(dplyr::filter(d_trawl2, depth <35), aes(SST,temperature_c, colour = DOY)) + geom_point()
Explore mean SST
ggplot(d_trawl2, aes_string("X", "Y", colour = "meanSST")) +
geom_point() +
facet_wrap(~year) +
coord_fixed() +
scale_color_viridis_c()
ggplot(d_trawl2, aes_string("depth_m", "meanSST", colour = "meanSST")) +
geom_point() +
facet_wrap(~ssid, scales = "free") +
# coord_fixed() +
scale_color_viridis_c()
ggplot(dplyr::filter(d_trawl2, shallow==1), aes(SST_scaled,temperature_c, colour = -depth)) + geom_point()
Does a fixed interaction with SST_scaled do any better?
Rerun main model with SST dataset for comparisons
spde <- make_mesh(d_trawl2, c("X", "Y"), n_knots = 800)
m_temp_new <- sdmTMB(
temperature_c ~ 0 + as.factor(year),
time_varying = ~ 0 + depth_scaled + depth_scaled2,
data = d_trawl2, spde = spde,
ar1_fields = TRUE,
include_spatial = TRUE,
time = "year",
family = gaussian(link = "identity"),
nlminb_loops = 2,
newton_steps = 1,
silent = FALSE
)
saveRDS(m_temp_new, file = (here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-800kn-sstdata.rds")))
m_temp_new
AIC(m_temp_new)
spde <- make_mesh(d_trawl2, c("X", "Y"), n_knots = 800)
m_temp_sst <- sdmTMB(
temperature_c ~ 0 + as.factor(year) + SST_scaled * shallow,
data = d_trawl2, spde = spde,
time_varying = ~ 0 + depth_scaled + depth_scaled2,
ar1_fields = TRUE,
include_spatial = TRUE,
time = "year",
family = gaussian(link = "identity"),
nlminb_loops = 2,
newton_steps = 1,
silent = TRUE
)
saveRDS(m_temp_sst, file = (here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-800kn-SST.rds")))
m_temp_sst
AIC(m_temp_sst)
m_temp_new <- readRDS((here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-800kn-sstdata.rds")))
m_temp_sst <- readRDS((here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-800kn-SST.rds")))
get_diag(m_temp_new, response = "temperature_c")
get_diag(m_temp_sst, response = "temperature_c")
try with fewer knots
spde3 <- sdmTMB::make_mesh(d_trawl2, c("X", "Y"), n_knots = 400, type = "kmeans")
m_temp1aic <- sdmTMB(
temperature_c ~ 0 + as.factor(year),
time_varying = ~ 0 + depth_scaled + depth_scaled2,
data = d_trawl2, spde = spde3,
ar1_fields = TRUE,
include_spatial = TRUE,
time = "year",
family = gaussian(link = "identity"),
nlminb_loops = 2,
newton_steps = 1,
silent = FALSE
)
saveRDS(m_temp1aic, file = (here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-400kn-new.rds")))
m_temp1aic <- readRDS((here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-400kn-new.rds")))
# m_temp1aic<- sdmTMB:::update_model(m_temp1aic)
AIC(m_temp1aic)
Run CV on main model
# spde2 <- make_spde(d_trawl2$X, d_trawl2$Y, n_knots = 800)
spde2 <- sdmTMB::make_mesh(d_trawl2, c("X", "Y"), n_knots = 200, type = "kmeans")
plot(spde2)
m_temp1 <- sdmTMB_cv(
temperature_c ~ 0 + as.factor(year),
time_varying = ~ 0 + depth_scaled + depth_scaled2,
data = d_trawl2, spde = spde3,
ar1_fields = TRUE,
include_spatial = TRUE,
time = "year",
family = gaussian(link = "identity"),
k_folds = 20,
# nlminb_loops = 2,
# newton_steps = 1,
silent = FALSE
)
# saveRDS(m_temp1, file = (here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-400kn-cv20.rds")))
saveRDS(m_temp1, file = (here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-150kn-cv20.rds")))
# m_temp1 <- readRDS((here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-300kn-cv.rds")))
m_temp1$fold_loglik
m_temp1$sum_loglik
m_temp1$max_gradients
m_temp1$pdHess
Check if DOY helps model
spde3 <- sdmTMB::make_mesh(d_trawl2, c("X", "Y"), n_knots = 400, type = "kmeans")
m_temp2aic <- sdmTMB(
temperature_c ~ 0 + as.factor(year) + DOY_scaled + DOY_scaled2,
data = d_trawl2, spde = spde3,
time_varying = ~ 0 + depth_scaled + depth_scaled2,
ar1_fields = TRUE,
include_spatial = TRUE,
time = "year",
family = gaussian(link = "identity"),
nlminb_loops = 2,
newton_steps = 1,
# k_folds = 20,
silent = TRUE
)
saveRDS(m_temp2aic, file = (here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-DOY-400kn.rds")))
m_temp2aic <- readRDS((here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-DOY-400kn.rds")))
m_temp2aic
AIC(m_temp2aic)
# m_temp2cv <- sdmTMB_cv(
# temperature_c ~ 0 + as.factor(year) + DOY_scaled + DOY_scaled2,
# data = d_trawl2, spde = spde2,
# time_varying = ~ 0 + depth_scaled + depth_scaled2,
# ar1_fields = TRUE,
# include_spatial = TRUE,
# time = "year",
# family = gaussian(link = "identity"),
# # nlminb_loops = 2,
# # newton_steps = 1,
# k_folds = 20,
# silent = TRUE
# )
#
# saveRDS(m_temp2, file = (here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-DOY-cv20.rds")))
# m_temp2cv <- readRDS((here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-DOY-cv20.rds")))
# m_temp2cv
# m_temp2cv$fold_loglik
# m_temp2cv$sum_loglik
# m_temp2cv$max_gradients
Check if SST helps model
m_temp3aic <- sdmTMB(
temperature_c ~ 0 + as.factor(year) + SST_scaled * shallow,
data = d_trawl2, spde = spde3,
time_varying = ~ 0 + depth_scaled + depth_scaled2,
ar1_fields = TRUE,
include_spatial = TRUE,
time = "year",
family = gaussian(link = "identity"),
nlminb_loops = 2,
newton_steps = 1,
silent = TRUE
)
saveRDS(m_temp3aic, file = (here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-SST-400kn.rds")))
# m_temp3aic <- readRDS(here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-SST-400kn.rds"))
m_temp3aic <- readRDS((here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-SST-400kn.rds")))
AIC(m_temp1aic, m_temp2aic, m_temp3aic)
m_temp3aic
Check if SST and DOY help model
m_temp4aic <- sdmTMB(
temperature_c ~ 0 + as.factor(year) + SST_scaled * DOY_scaled * shallow,
data = d_trawl2, spde = spde3,
time_varying = ~ 0 + depth_scaled + depth_scaled2,
ar1_fields = TRUE,
include_spatial = TRUE,
time = "year",
family = gaussian(link = "identity"),
nlminb_loops = 2,
newton_steps = 1,
silent = TRUE
)
saveRDS(m_temp4aic, file = (here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-SST-DOY-400kn.rds")))
m_temp4aic <- readRDS(here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-SST-DOY-400kn.rds"))
# m_temp4aic
AIC(m_temp1aic, m_temp2aic, m_temp3aic, m_temp4aic)
m_temp3 <- sdmTMB_cv(
temperature_c ~ 0 + as.factor(year) + SST_scaled * shallow,
data = d_trawl2, spde = spde2,
time_varying = ~ 0 + depth_scaled + depth_scaled2,
ar1_fields = TRUE,
include_spatial = TRUE,
time = "year",
family = gaussian(link = "identity"),
# nlminb_loops = 2,
# newton_steps = 1,
k_folds = 20,
silent = TRUE
)
# 500 knot early version
# saveRDS(m_temp3, file = (here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-SST.rds")))
# new 800 knot version for cross validation
saveRDS(m_temp3, file = (here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-SST-cv20.rds")))
# m_temp3 <- readRDS(here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-SST-cv20.rds"))
# m_temp3
m_temp3$fold_loglik
m_temp3$sum_loglik
m_temp3$max_gradients
Check if SST and DOY help model using cv
m_temp4 <- sdmTMB_cv(
temperature_c ~ 0 + as.factor(year) + SST_scaled * shallow + SST_scaled * DOY_scaled * shallow,
data = d_trawl2, spde = spde2,
time_varying = ~ 0 + depth_scaled + depth_scaled2,
ar1_fields = TRUE,
include_spatial = TRUE,
time = "year",
family = gaussian(link = "identity"),
# nlminb_loops = 2,
# newton_steps = 1,
k_folds = 20,
silent = TRUE
)
# 500 knot early version
# saveRDS(m_temp3, file = (here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-SST.rds")))
# new 800 knot version for cross validation
saveRDS(m_temp4, file = (here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-SST-DOY-cv20.rds")))
# m_temp4 <- readRDS(here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-SST-DOY-cv20.rds"))
# m_temp4
m_temp4$fold_loglik
m_temp4$sum_loglik
m_temp4$max_gradients
m_temp4$pdHess
m_temp_sst <- readRDS(here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-SST.rds"))
m_temp_sst <- sdmTMB:::update_model(m_temp_sst)
predictions3 <- predict(m_temp3)
predictions3$residuals <- residuals(m_temp3)
cor(predictions3$est, predictions3$temperature_c)^2 #
get_diag(m_temp_new, response = "temperature_c", variable = "depth_scaled", colour_var = "depth", start_year = 2002)
get_diag(m_temp_new, response = "temperature_c", variable = "DOY", colour_var = "depth", start_year = 2002, ssids = c(1,3,4, 16))
Plotting predictions from saved temperature model
# m_temp_rw <- readRDS((here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-500kn.rds")))
nd_all <- readRDS(here::here("analysis/VOCC/data/nd_all_synoptic.rds")) %>% filter(year<2019)
unique(nd_all$year)
glimpse(nd_all)
p_nd <- predict(m_temp_rw, newdata = nd_all)
# saveRDS(p_nd, file = "../VOCC/data/predicted_temp_allyears_800kn.rds")
p_nd <- readRDS("../VOCC/data/predicted_temp_allyears_800kn.rds")
plot_map_raster(p_nd, "est") +
scale_fill_viridis_c(trans = sqrt) +
facet_wrap(~year) +
ggtitle("Prediction (fixed effects + all random effects)")
plot_facet_map(p_nd, "est", viridis_option = "D",
raster_limits = c(min(p_nd$est), quantile(p_nd$est, 0.99)),
transform_col = sqrt
)
p_nd1 <- p_nd %>% filter(ssid!=16) %>% filter(ssid!=4)
breaks <- scales::trans_breaks(no_trans[["transform"]], no_trans[["inverse"]],n = 6)
labels <- function(x) format(x, digits = 2, scientific = FALSE)
pred_1 <- plot_map_raster(p_nd1, "est") +
scale_fill_viridis_c(labels = labels,
option = "C",
limits = c(4, 11),
na.value = "red") +
facet_wrap(~year) +
labs(fill = "°C") +
ggtitle("Prediction (fixed effects + all random effects)")
pred_1
Time-varying depth
# pd <- expand.grid(
# depth_scaled = seq(min(m_temp_rw$data$depth_scaled),
# max(m_temp_rw$data$depth_scaled), length.out = 50),
# year = unique(m_temp_rw$data$year)
# )
# pd$depth_scaled2 <- pd$depth_scaled^2
#
# # p <- predict(m_temp_rw, newdata = pd, se_fit = TRUE, re_form = NA, xy_cols = c("X", "Y"))
# p <- predict(m_temp_new, newdata = pd, se_fit = TRUE, re_form = NA)
#
# p$Depth <- exp(p$depth_scaled*0.661 + 5.07)
# p$Temperature <- p$est
# #
# p$lowCI <- (p$est - 1.96 * p$est_se)
# p$highCI <- (p$est + 1.96 * p$est_se)
# #
# saveRDS(p, here::here("analysis/tmb-sensor-explore/data/depth-tv-m-temp-rw-800-kn.rds"))
p <- readRDS(here::here("analysis/tmb-sensor-explore/data/depth-tv-m-temp-rw-800-kn.rds"))
ggplot(filter(p, Depth >15 & year >2007), aes(Depth, Temperature,
ymin = lowCI,
ymax = highCI,
group = as.factor(year), fill=year, colour = year)) +
geom_line(size = 0.5, alpha =0.85) +
# geom_smooth(method= "loess", size = 0.4, alpha =0.85, se =F ) +
# geom_smooth(span = 0.5, size = 0.4, alpha =0.85, se =F ) +
scale_fill_viridis_c(option = "C") +
scale_colour_viridis_c(option = "C") +
# ylab("DO") +
geom_ribbon(alpha = 0.2, colour = NA) +
# scale_y_continuous(limits = c(0,4.5)) +
scale_x_reverse(limits = c(600,20)) +
ylim( 2, 15) +
coord_flip() +
# scale_x_continuous(limits = c(18,210)) +
gfplot::theme_pbs() +
theme(legend.position = c(.8,.4))
ggsave(filename = here::here("analysis/tmb-sensor-explore/figs/depth-tv-m-temp-rw-800-kn.png"), width = 3, height = 2.5)
# ggplot(filter(p, Depth >15 & year == 2018), aes(Depth, Temperature,
# ymin = (Temperature - 1.96 * est_se),
# ymax = (Temperature + 1.96 * est_se))) +
# geom_line(size = 1.5, alpha =0.85) +
# # scale_color_viridis_c(option = "C") +
# geom_ribbon(alpha=0.2) +
# scale_x_reverse(limits = c(210,20)) +
# ylim(5, 15) +
# coord_flip() +
# gfplot::theme_pbs() +
# theme(legend.position = c(.8,.4))
p_nd4 <- p_nd %>% filter(ssid==4)
pred_4 <- plot_map_raster(p_nd4, "est") +
scale_fill_viridis_c(labels = labels,
option = "C",
limits = c(5, 11),
na.value = "black")+
facet_wrap(~year) +
theme(legend.position = c(.85,.12)) +
labs(fill = "°C")
pred_4
plot_map_raster(p_nd, "est_non_rf") +
ggtitle("Prediction (fixed effects only)") +
facet_wrap(~year) +
scale_fill_viridis_c()
plot_map_raster(p_nd, "est_rf") +
ggtitle("Prediction (all random effects only)") +
facet_wrap(~year) +
scale_fill_viridis_c()
(omega_s_all <- plot_map_raster(p_nd, "omega_s") +
ggtitle("Spatial random effects only") +
scale_fill_gradient2( low = "Steel Blue 4",
high = "Red 3", limits = c(-1.5,1.5),
na.value = "red") + labs(fill = "") + theme(legend.position = c(0.2, 0.3))
)
ggsave(filename = here::here("analysis/tmb-sensor-explore/figs/omegas-temp-rw-800-kn.png"), width = 3, height = 3)
omega_s_all
plot_map_raster(p_nd, "epsilon_st") +
ggtitle("Spatiotemporal random effects only") +
facet_wrap(~year) + scale_fill_gradient2( low = "Steel Blue 4",
high = "Red 3", limits = c(-1.5,1.5),
na.value = "red")
p_nd <- readRDS("../VOCC/data/predicted_temp_allyears_800kn.rds")
p_nd1 <- p_nd %>% filter(ssid!=16) %>% filter(ssid!=4) %>% filter(year>2007)
(st_1n3 <- plot_map_raster(p_nd1, "epsilon_st") +
ggtitle("Spatiotemporal random effects") +
facet_grid(~year) +
scale_fill_gradient2( low = "Steel Blue 4",
high = "Red 3", limits = c(-1.5,1.5), #trans = sqrt,
na.value = "red") +
labs(fill = "epsilon") +
theme(
legend.position = "none"))
ggsave(filename = here::here("analysis/tmb-sensor-explore/figs/temp-epsilon-1n3.png"), width =6, height = 2)
# png(
# file = "figs/temp-predictions-1n3.png",
# res = 400,
# units = "in",
# width = 8,
# height = 8
# )
# st_1n3
# dev.off()
p_nd <- readRDS("../VOCC/data/predicted_temp_allyears_800kn.rds")
p_nd4 <- p_nd %>% filter(ssid==4) %>% filter(year>2007)#%>% filter(ssid!=4)
(st_4 <- plot_map_raster(p_nd4, "epsilon_st") +
#ggtitle(" ") +
facet_grid(~year) +
scale_fill_gradient2( low = "Steel Blue 4",
high = "Red 3", limits = c(-1.5,1.5), #trans = sqrt,
na.value = "red") +
labs(fill = " ") + theme( legend.position = "none") )
ggsave(filename = here::here("analysis/tmb-sensor-explore/figs/temp-epsilon-4.png"), width = 6, height = 2)
st_4
p_nd <- readRDS("../VOCC/data/predicted_temp_allyears_800kn.rds")
p_nd16 <- p_nd %>% filter(ssid==16) %>% filter(year>2007)#%>% filter(ssid!=4)
(st_16 <- plot_map_raster(p_nd16, "epsilon_st") +
#ggtitle(" ") +
facet_grid(~year) +
scale_fill_gradient2( low = "Steel Blue 4",
high = "Red 3", limits = c(-1.5,1.5),
na.value = "red") +
xlim(20,360) + ylim(5850,6050)+
# xlim(65,345) + ylim(5800,6100)+
theme(panel.border = element_blank(), legend.position = "none") )
ggsave(filename = here::here("analysis/tmb-sensor-explore/figs/temp-epsilon-16.png"), width = 6, height = 1)
st_16
png(
file = "figs/temp-predictions.png",
res = 600,
units = "in",
width = 11,
height = 10
)
gridExtra::grid.arrange(
grobs = list(pred_1, st_1n3, omega_s_all, pred_4, st_4, st_16 ),
widths = c(2, 2, 1.75),
heights = c(2, 1.35),
layout_matrix = rbind(c(1, 2, 3),
c(4, 5, 6)),
top = grid::textGrob("Bottom temperature estimates from sdmTMB model")
)
# gridExtra::grid.arrange( nrow = 2, heights = c(2,1.5),
# top = grid::textGrob("Dissolved oxygen estimates from sdmTMB model"))
dev.off()
BOOSTED REGRESSION MODEL
library(gbm)
m <- gbm( temperature_c ~ depth_scaled + X + Y, data = d_trawl, n.trees = 2000, interaction.depth = 3, shrinkage = 0.02)
m2 <- gbm( temperature_c ~ SST_scaled + depth_scaled + X + Y, data = d_trawl, n.trees = 2000, interaction.depth = 3, shrinkage = 0.02)
plot(m,i.var=2)
plot(m,i.var=3)
plot(m,i.var=4)
plot(m,i.var=5)
plot(m2,i.var=1)
plot(m2,i.var=2)
plot(m,i.var=3)
plot(m,i.var=4)
plot(m2,i.var=1:2)
plot(m2,i.var=c(1,3))
plot(m2,i.var=c(1,4))
d_trawl$r <- predict(m, n.trees = 2000) - d_trawl1$temperature_c
ggplot(d_trawl, aes(depth_scaled, r)) +
geom_point(alpha=0.4) +
ylim(-8,8) +
geom_smooth()
CENTER OF GRAVITY
p_cog <- predict(m_temp_rw, newdata = nd_all, return_tmb_object = TRUE)
temp_cog <- get_cog(p_cog)
temp_cog %>%
reshape2::dcast(year ~ coord, value.var = "est") %>%
ggplot(aes(X, Y, colour = year)) + geom_path(arrow = arrow()) +
scale_color_viridis_c()
temp_cog %>%
ggplot(aes(year, est, ymin = lwr, ymax = upr)) +
geom_line() +
geom_ribbon(alpha = 0.2) +
facet_wrap(~coord, scales = "free_y")
pbs-assess/gfranges documentation built on Dec. 13, 2021, 4:50 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Fit spatial temporal models to bottom temperature
library(tidyverse) library(sdmTMB) # library(future) library(dplyr) library(gfranges) plot_map_raster <- function(dat, column = "est") { ggplot(dat, aes_string("X", "Y", fill = column)) + geom_raster() + coord_fixed()+ gfplot::theme_pbs() + theme( axis.title = element_blank(), axis.ticks = element_blank(), axis.text = element_blank() ) }
Model diagnostics function
get_diag <- function(m, response = "temperature_c", variable = "depth_scaled", colour_var = "depth", start_year = 2007, ssids = c(1,3,4,16)) { predictions <- predict(m) predictions$residuals <- residuals(m) predictions <- predictions %>% filter (year >= start_year) plot_map <- function(dat, column = "est") { ggplot(dat, aes_string("X", "Y", colour = column)) + geom_point() + coord_fixed()+ gfplot::theme_pbs() } g <- plot_map(predictions, "est") + scale_colour_viridis_c() + ggtitle("Prediction (fixed effects + all random effects)") print(g) g <- plot_map(predictions, "est_non_rf") + ggtitle("Prediction (fixed effects only)") + scale_colour_viridis_c() print(g) g <- plot_map(predictions, "est_rf") + ggtitle("All random effects only") + scale_colour_gradient2() print(g) g <- plot_map(predictions, "omega_s") + ggtitle("Spatial random effects only") + scale_colour_gradient2() print(g) g <- plot_map(predictions, "epsilon_st") + ggtitle("Spatiotemporal random effects only") + facet_wrap(~year) + scale_colour_gradient2() print(g) g <- ggplot(predictions, aes_string("est", response)) + geom_point(alpha = 0.2) + facet_wrap(~year) + coord_fixed() + geom_abline()+ gfplot::theme_pbs() print(g) g <- ggplot(predictions, aes(est, residuals)) + geom_point(alpha = 0.2) + geom_smooth()+ gfplot::theme_pbs() print(g) g <- ggplot(predictions, aes_string(variable, "residuals", colour = colour_var)) + geom_point(alpha = 0.4, size = 1.2) + geom_smooth(colour="grey", size = 1.2) + scale_colour_viridis_c(option = "B", direction = -1, begin = 0, end = 0.7, limits= c(min(predictions[colour_var]), max(predictions[colour_var]))) + #, trans= sqrt # facet_wrap(~year, scales = "free_x")+ gfplot::theme_pbs() + theme( axis.text.y = element_text(size = 14), axis.text.x = element_blank(), axis.ticks.x = element_blank()) print(g) g <- ggplot(predictions, aes_string(variable, "residuals", colour = colour_var)) + geom_point(alpha = 0.4, size = 0.7) + geom_smooth(colour="grey", size = 0.7) + scale_colour_viridis_c(option = "B", direction = -1, begin = 0, end = 0.7, limits= c(min(predictions[colour_var]), max(predictions[colour_var]))) + #, trans= sqrt facet_wrap(~year, scales = "free_x")+ gfplot::theme_pbs() + theme(axis.text.x = element_blank(), axis.ticks.x = element_blank()) print(g) g <- ggplot(filter(predictions, ssid %in% ssids), aes_string(variable, "residuals", colour = colour_var)) + geom_point(alpha = 0.4, size = 1) + geom_smooth(colour="grey", size = 1) + scale_colour_viridis_c(option = "B", direction = -1, begin = 0, end = 0.7, limits= c(min(predictions[colour_var]), max(predictions[colour_var]))) + #, trans= sqrt facet_wrap(~ssid, scales = "free_x")+ gfplot::theme_pbs() print(g) g <- ggplot(predictions, aes(X, Y, colour = residuals)) + geom_point(alpha = 0.2) + coord_fixed() + scale_color_gradient2() + facet_wrap(~year)+ gfplot::theme_pbs() print(g) aic <- function(m) { k <- length(m$model$par) nll <- m$model$objective 2 * k - 2 * (-nll) } print("R^2:") r2 <- cor(predictions$est, predictions[[response]])^2 print(r2) print("") print("AIC:") print(aic(m)) }
Load data if not using SST
all_sensor <- readRDS(here::here("analysis/tmb-sensor-explore/data/dat-sensor-trawl-processed.rds")) # no SST all_depth <- all_sensor %>% # filter(depth_max>10) %>% # maybe we should filter incase missing depth values also indicate a sensor fail? # filter missing location data and trial year dplyr::filter(!is.na(latitude), !is.na(longitude)) %>% # convert lat and lon to UTMs dplyr::mutate(X = longitude, Y = latitude) %>% gfplot:::ll2utm(., utm_zone = 9) %>% # interpolate missing depth values dplyr::rename(depth = depth_m) %>% gfplot:::interp_survey_bathymetry() d_trawl <- all_depth$data %>% # remove obvious sensor fails, d_trawl %>% filter(depth_max<10), fishing event ids = 481861, 2179089 filter(depth_max>10) %>% gfplot:::scale_survey_predictors() d_trawl <- d_trawl[!is.na(d_trawl$depth), ] d_trawl <- d_trawl[!is.na(d_trawl$temperature_c), ] d_trawl <- mutate(d_trawl, DOY = lubridate::yday(date), DOY_scaled = arm::rescale(lubridate::yday(date)), DOY_scaled2 = arm::rescale(lubridate::yday(date))^2)
EXLORE TEMPERATURE at DEPTH
ggplot(d_trawl, aes_string("X", "Y", colour = "temperature_c")) + geom_point() + facet_wrap(~year) + coord_fixed() + scale_color_viridis_c()
ggplot(d_trawl, aes(DOY, temperature_c, colour = -depth, alpha=0.5)) + geom_point() + facet_wrap(~year)
Current temperature model
spde <- sdmTMB::make_mesh(d_trawl, c("X", "Y"), n_knots = 800, type = "kmeans") plot(spde)
Original model code
# spde <- make_spde(d_trawl$X, d_trawl$Y, n_knots = 800) # m_temp_rw <- sdmTMB( # temperature_c ~ 0 + as.factor(year), #+ DOY_scaled + as.factor(ssid) # time_varying = ~ 0 + depth_scaled + depth_scaled2, # data = d_trawl, # ar1_fields = TRUE, # changed AR1 to true # include_spatial = TRUE, # time = "year", spde = spde, # family = gaussian(link = "identity"), # nlminb_loops = 2, # newton_steps = 1, # silent = FALSE # ) # saveRDS(m_temp_rw, file = (here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-800kn.rds")))
New CV model code
m_temp_cv <- sdmTMB_cv( temperature_c ~ 0 + as.factor(year), #+ DOY_scaled + as.factor(ssid) time_varying = ~ 0 + depth_scaled + depth_scaled2, data = d_trawl, ar1_fields = TRUE, # changed AR1 to true include_spatial = TRUE, time = "year", spde = spde, family = gaussian(link = "identity"), k_folds = 5, # nlminb_loops = 2, # newton_steps = 1, silent = FALSE ) saveRDS(m_temp_cv, file = (here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-800kn-cv.rds"))) m_temp_cv <- readRDS((here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-800kn-cv.rds"))) m_temp_cv$fold_loglik m_temp_cv$sum_loglik m_temp_cv$max_gradients
Check predictions against raw data
# m_temp_rw <- readRDS(here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-500kn-2.rds")) m_temp_rw <- readRDS(here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-800kn.rds")) predictions <- predict(m_temp_rw) predictions$residuals <- residuals(m_temp_rw) cor(predictions$est, predictions$temperature_c)^2 sum(abs(predictions$residuals)) # 2908.879-2789.322 # 119.557/2908.879 # 4% reduction in abs errors
# test against ROMS events_roms <- readRDS("data/roms_temp_by_events2.rds") %>% select(fishing_event_id, roms, year) predictions <- left_join(predictions, events_roms) cor(predictions$roms, predictions$est, use = "pairwise.complete.obs") cor(predictions$roms, predictions$temperature_c, use = "pairwise.complete.obs") ggplot(predictions, aes(roms, est, colour = -depth)) + geom_point(alpha=0.3) + facet_wrap(~year) + coord_cartesian(xlim = c(2,13), ylim = c(2,13)) + geom_abline(intercept = 0, slope = 1) + gfplot::theme_pbs() ggplot(predictions, aes(roms, temperature_c, colour = DOY)) + geom_point(alpha=0.3) + facet_wrap(~year) + coord_cartesian(xlim = c(2,13), ylim = c(2,13)) + geom_abline(intercept = 0, slope = 1) + scale_colour_viridis_c() + gfplot::theme_pbs() ggplot(predictions, aes(DOY, (roms-temperature_c), colour = -depth)) + geom_point(alpha=0.3) + facet_wrap(~year) + #coord_cartesian(xlim = c(2,13), ylim = c(2,13)) + geom_abline(intercept = 0, slope = 0) + scale_colour_viridis_c() + gfplot::theme_pbs() ggplot(predictions, aes(DOY, (est-roms), colour = -depth)) + geom_point(alpha=0.3) + facet_wrap(~year) + #coord_cartesian(xlim = c(2,13), ylim = c(2,13)) + geom_abline(intercept = 0, slope = 0) + scale_colour_viridis_c() + gfplot::theme_pbs() p_nd <- readRDS((here::here("analysis/VOCC/data/predicted-DO-2020-06-20-more2016-roms.rds"))) cor(p_nd$roms_temp, p_nd$temp, use = "pairwise.complete.obs") ggplot(p_nd, aes(roms_temp, temp, colour = -depth)) + geom_point(alpha=0.04) + facet_wrap(~year) + coord_cartesian(xlim = c(2,13), ylim = c(2,13)) + geom_abline(intercept = 0, slope = 1) + gfplot::theme_pbs()
m_temp_rw <- readRDS((here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-500kn-2.rds"))) get_diag(m_temp_rw, response = "temperature_c", variable = "depth_scaled")
WOULD SST HELP?
Load and prep SST data
all_sensor <- readRDS(here::here("analysis/tmb-sensor-explore/data/dat-sensor-trawl-meanSST-all.rds")) d_trawl2 <-all_sensor # all_sensor <- readRDS(here::here("analysis/tmb-sensor-explore/data/dat-sensor-trawl-processed.rds")) # no SST # glimpse(all_sensor) # # all_depth <- all_sensor %>% # # filter(depth_max>10) %>% # maybe we should filter incase missing depth values also indicate a sensor fail? # # filter missing location data and trial year # dplyr::filter(!is.na(latitude), !is.na(longitude)) %>% # # # convert lat and lon to UTMs # dplyr::mutate(X = longitude, Y = latitude) %>% # gfplot:::ll2utm(., utm_zone = 9) %>% # # # interpolate missing depth values # dplyr::rename(depth = depth_m) %>% # gfplot:::interp_survey_bathymetry() # d_trawl <- all_depth$data %>% # # remove obvious sensor fails, d_trawl %>% filter(depth_max<10), fishing event ids = 481861, 2179089 # filter(depth_max>10) %>% # gfplot:::scale_survey_predictors() olddata <- readRDS((here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-400kn.rds"))) d_trawl <- olddata$data d_trawl <- d_trawl[!is.na(d_trawl$depth), ] d_trawl <- d_trawl[!is.na(d_trawl$temperature_c), ] d_trawl2 <- d_trawl[!is.na(d_trawl$SST), ] d_trawl2 <- d_trawl2[!is.na(d_trawl2$meanSST), ] d_trawl2 <- mutate(d_trawl2, DOY = lubridate::yday(date), DOY_scaled = arm::rescale(DOY), DOY_scaled2 = DOY_scaled^2, depth_scaled = arm::rescale(log(depth)), SST_scaled = arm::rescale(meanSST), depth_scaled3 = depth_scaled^3, shallow = ifelse(depth>35,0,1), deep = ifelse(depth<35,0,1) ) # d_trawl2[is.na(d_trawl2$shallow), ] # d_trawl2[is.na(d_trawl2$SST_scaled), ]
Explore SST on survey day
ggplot(d_trawl2, aes_string("X", "Y", colour = "SST")) + geom_point() + facet_wrap(~year) + coord_fixed() + scale_color_viridis_c()
Is bottom temp related to SST in shallower samples?
ggplot(dplyr::filter(d_trawl2, depth <35), aes(SST,temperature_c, colour = -depth)) + geom_point() ggplot(dplyr::filter(d_trawl2, depth <35), aes(SST,temperature_c, colour = DOY)) + geom_point()
Explore mean SST
ggplot(d_trawl2, aes_string("X", "Y", colour = "meanSST")) + geom_point() + facet_wrap(~year) + coord_fixed() + scale_color_viridis_c()
ggplot(d_trawl2, aes_string("depth_m", "meanSST", colour = "meanSST")) + geom_point() + facet_wrap(~ssid, scales = "free") + # coord_fixed() + scale_color_viridis_c()
ggplot(dplyr::filter(d_trawl2, shallow==1), aes(SST_scaled,temperature_c, colour = -depth)) + geom_point()
Does a fixed interaction with SST_scaled do any better? Rerun main model with SST dataset for comparisons
spde <- make_mesh(d_trawl2, c("X", "Y"), n_knots = 800) m_temp_new <- sdmTMB( temperature_c ~ 0 + as.factor(year), time_varying = ~ 0 + depth_scaled + depth_scaled2, data = d_trawl2, spde = spde, ar1_fields = TRUE, include_spatial = TRUE, time = "year", family = gaussian(link = "identity"), nlminb_loops = 2, newton_steps = 1, silent = FALSE ) saveRDS(m_temp_new, file = (here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-800kn-sstdata.rds"))) m_temp_new AIC(m_temp_new)
spde <- make_mesh(d_trawl2, c("X", "Y"), n_knots = 800) m_temp_sst <- sdmTMB( temperature_c ~ 0 + as.factor(year) + SST_scaled * shallow, data = d_trawl2, spde = spde, time_varying = ~ 0 + depth_scaled + depth_scaled2, ar1_fields = TRUE, include_spatial = TRUE, time = "year", family = gaussian(link = "identity"), nlminb_loops = 2, newton_steps = 1, silent = TRUE ) saveRDS(m_temp_sst, file = (here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-800kn-SST.rds"))) m_temp_sst AIC(m_temp_sst)
m_temp_new <- readRDS((here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-800kn-sstdata.rds"))) m_temp_sst <- readRDS((here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-800kn-SST.rds"))) get_diag(m_temp_new, response = "temperature_c") get_diag(m_temp_sst, response = "temperature_c")
try with fewer knots
spde3 <- sdmTMB::make_mesh(d_trawl2, c("X", "Y"), n_knots = 400, type = "kmeans") m_temp1aic <- sdmTMB( temperature_c ~ 0 + as.factor(year), time_varying = ~ 0 + depth_scaled + depth_scaled2, data = d_trawl2, spde = spde3, ar1_fields = TRUE, include_spatial = TRUE, time = "year", family = gaussian(link = "identity"), nlminb_loops = 2, newton_steps = 1, silent = FALSE ) saveRDS(m_temp1aic, file = (here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-400kn-new.rds"))) m_temp1aic <- readRDS((here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-400kn-new.rds"))) # m_temp1aic<- sdmTMB:::update_model(m_temp1aic) AIC(m_temp1aic)
Run CV on main model
# spde2 <- make_spde(d_trawl2$X, d_trawl2$Y, n_knots = 800) spde2 <- sdmTMB::make_mesh(d_trawl2, c("X", "Y"), n_knots = 200, type = "kmeans") plot(spde2)
m_temp1 <- sdmTMB_cv( temperature_c ~ 0 + as.factor(year), time_varying = ~ 0 + depth_scaled + depth_scaled2, data = d_trawl2, spde = spde3, ar1_fields = TRUE, include_spatial = TRUE, time = "year", family = gaussian(link = "identity"), k_folds = 20, # nlminb_loops = 2, # newton_steps = 1, silent = FALSE ) # saveRDS(m_temp1, file = (here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-400kn-cv20.rds"))) saveRDS(m_temp1, file = (here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-150kn-cv20.rds"))) # m_temp1 <- readRDS((here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-300kn-cv.rds"))) m_temp1$fold_loglik m_temp1$sum_loglik m_temp1$max_gradients m_temp1$pdHess
Check if DOY helps model
spde3 <- sdmTMB::make_mesh(d_trawl2, c("X", "Y"), n_knots = 400, type = "kmeans") m_temp2aic <- sdmTMB( temperature_c ~ 0 + as.factor(year) + DOY_scaled + DOY_scaled2, data = d_trawl2, spde = spde3, time_varying = ~ 0 + depth_scaled + depth_scaled2, ar1_fields = TRUE, include_spatial = TRUE, time = "year", family = gaussian(link = "identity"), nlminb_loops = 2, newton_steps = 1, # k_folds = 20, silent = TRUE ) saveRDS(m_temp2aic, file = (here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-DOY-400kn.rds"))) m_temp2aic <- readRDS((here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-DOY-400kn.rds"))) m_temp2aic AIC(m_temp2aic) # m_temp2cv <- sdmTMB_cv( # temperature_c ~ 0 + as.factor(year) + DOY_scaled + DOY_scaled2, # data = d_trawl2, spde = spde2, # time_varying = ~ 0 + depth_scaled + depth_scaled2, # ar1_fields = TRUE, # include_spatial = TRUE, # time = "year", # family = gaussian(link = "identity"), # # nlminb_loops = 2, # # newton_steps = 1, # k_folds = 20, # silent = TRUE # ) # # saveRDS(m_temp2, file = (here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-DOY-cv20.rds"))) # m_temp2cv <- readRDS((here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-DOY-cv20.rds"))) # m_temp2cv # m_temp2cv$fold_loglik # m_temp2cv$sum_loglik # m_temp2cv$max_gradients
Check if SST helps model
m_temp3aic <- sdmTMB( temperature_c ~ 0 + as.factor(year) + SST_scaled * shallow, data = d_trawl2, spde = spde3, time_varying = ~ 0 + depth_scaled + depth_scaled2, ar1_fields = TRUE, include_spatial = TRUE, time = "year", family = gaussian(link = "identity"), nlminb_loops = 2, newton_steps = 1, silent = TRUE ) saveRDS(m_temp3aic, file = (here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-SST-400kn.rds"))) # m_temp3aic <- readRDS(here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-SST-400kn.rds")) m_temp3aic <- readRDS((here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-SST-400kn.rds"))) AIC(m_temp1aic, m_temp2aic, m_temp3aic) m_temp3aic
Check if SST and DOY help model
m_temp4aic <- sdmTMB( temperature_c ~ 0 + as.factor(year) + SST_scaled * DOY_scaled * shallow, data = d_trawl2, spde = spde3, time_varying = ~ 0 + depth_scaled + depth_scaled2, ar1_fields = TRUE, include_spatial = TRUE, time = "year", family = gaussian(link = "identity"), nlminb_loops = 2, newton_steps = 1, silent = TRUE ) saveRDS(m_temp4aic, file = (here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-SST-DOY-400kn.rds"))) m_temp4aic <- readRDS(here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-SST-DOY-400kn.rds")) # m_temp4aic AIC(m_temp1aic, m_temp2aic, m_temp3aic, m_temp4aic)
m_temp3 <- sdmTMB_cv( temperature_c ~ 0 + as.factor(year) + SST_scaled * shallow, data = d_trawl2, spde = spde2, time_varying = ~ 0 + depth_scaled + depth_scaled2, ar1_fields = TRUE, include_spatial = TRUE, time = "year", family = gaussian(link = "identity"), # nlminb_loops = 2, # newton_steps = 1, k_folds = 20, silent = TRUE ) # 500 knot early version # saveRDS(m_temp3, file = (here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-SST.rds"))) # new 800 knot version for cross validation saveRDS(m_temp3, file = (here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-SST-cv20.rds"))) # m_temp3 <- readRDS(here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-SST-cv20.rds")) # m_temp3 m_temp3$fold_loglik m_temp3$sum_loglik m_temp3$max_gradients
Check if SST and DOY help model using cv
m_temp4 <- sdmTMB_cv( temperature_c ~ 0 + as.factor(year) + SST_scaled * shallow + SST_scaled * DOY_scaled * shallow, data = d_trawl2, spde = spde2, time_varying = ~ 0 + depth_scaled + depth_scaled2, ar1_fields = TRUE, include_spatial = TRUE, time = "year", family = gaussian(link = "identity"), # nlminb_loops = 2, # newton_steps = 1, k_folds = 20, silent = TRUE ) # 500 knot early version # saveRDS(m_temp3, file = (here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-SST.rds"))) # new 800 knot version for cross validation saveRDS(m_temp4, file = (here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-SST-DOY-cv20.rds"))) # m_temp4 <- readRDS(here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-SST-DOY-cv20.rds")) # m_temp4 m_temp4$fold_loglik m_temp4$sum_loglik m_temp4$max_gradients m_temp4$pdHess
m_temp_sst <- readRDS(here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-SST.rds")) m_temp_sst <- sdmTMB:::update_model(m_temp_sst)
predictions3 <- predict(m_temp3) predictions3$residuals <- residuals(m_temp3) cor(predictions3$est, predictions3$temperature_c)^2 #
get_diag(m_temp_new, response = "temperature_c", variable = "depth_scaled", colour_var = "depth", start_year = 2002) get_diag(m_temp_new, response = "temperature_c", variable = "DOY", colour_var = "depth", start_year = 2002, ssids = c(1,3,4, 16))
Plotting predictions from saved temperature model
# m_temp_rw <- readRDS((here::here("analysis/tmb-sensor-explore/models/model-temp-all-years-500kn.rds"))) nd_all <- readRDS(here::here("analysis/VOCC/data/nd_all_synoptic.rds")) %>% filter(year<2019) unique(nd_all$year) glimpse(nd_all) p_nd <- predict(m_temp_rw, newdata = nd_all) # saveRDS(p_nd, file = "../VOCC/data/predicted_temp_allyears_800kn.rds")
p_nd <- readRDS("../VOCC/data/predicted_temp_allyears_800kn.rds") plot_map_raster(p_nd, "est") + scale_fill_viridis_c(trans = sqrt) + facet_wrap(~year) + ggtitle("Prediction (fixed effects + all random effects)")
plot_facet_map(p_nd, "est", viridis_option = "D", raster_limits = c(min(p_nd$est), quantile(p_nd$est, 0.99)), transform_col = sqrt )
p_nd1 <- p_nd %>% filter(ssid!=16) %>% filter(ssid!=4) breaks <- scales::trans_breaks(no_trans[["transform"]], no_trans[["inverse"]],n = 6) labels <- function(x) format(x, digits = 2, scientific = FALSE) pred_1 <- plot_map_raster(p_nd1, "est") + scale_fill_viridis_c(labels = labels, option = "C", limits = c(4, 11), na.value = "red") + facet_wrap(~year) + labs(fill = "°C") + ggtitle("Prediction (fixed effects + all random effects)") pred_1
Time-varying depth
# pd <- expand.grid( # depth_scaled = seq(min(m_temp_rw$data$depth_scaled), # max(m_temp_rw$data$depth_scaled), length.out = 50), # year = unique(m_temp_rw$data$year) # ) # pd$depth_scaled2 <- pd$depth_scaled^2 # # # p <- predict(m_temp_rw, newdata = pd, se_fit = TRUE, re_form = NA, xy_cols = c("X", "Y")) # p <- predict(m_temp_new, newdata = pd, se_fit = TRUE, re_form = NA) # # p$Depth <- exp(p$depth_scaled*0.661 + 5.07) # p$Temperature <- p$est # # # p$lowCI <- (p$est - 1.96 * p$est_se) # p$highCI <- (p$est + 1.96 * p$est_se) # # # saveRDS(p, here::here("analysis/tmb-sensor-explore/data/depth-tv-m-temp-rw-800-kn.rds")) p <- readRDS(here::here("analysis/tmb-sensor-explore/data/depth-tv-m-temp-rw-800-kn.rds")) ggplot(filter(p, Depth >15 & year >2007), aes(Depth, Temperature, ymin = lowCI, ymax = highCI, group = as.factor(year), fill=year, colour = year)) + geom_line(size = 0.5, alpha =0.85) + # geom_smooth(method= "loess", size = 0.4, alpha =0.85, se =F ) + # geom_smooth(span = 0.5, size = 0.4, alpha =0.85, se =F ) + scale_fill_viridis_c(option = "C") + scale_colour_viridis_c(option = "C") + # ylab("DO") + geom_ribbon(alpha = 0.2, colour = NA) + # scale_y_continuous(limits = c(0,4.5)) + scale_x_reverse(limits = c(600,20)) + ylim( 2, 15) + coord_flip() + # scale_x_continuous(limits = c(18,210)) + gfplot::theme_pbs() + theme(legend.position = c(.8,.4)) ggsave(filename = here::here("analysis/tmb-sensor-explore/figs/depth-tv-m-temp-rw-800-kn.png"), width = 3, height = 2.5) # ggplot(filter(p, Depth >15 & year == 2018), aes(Depth, Temperature, # ymin = (Temperature - 1.96 * est_se), # ymax = (Temperature + 1.96 * est_se))) + # geom_line(size = 1.5, alpha =0.85) + # # scale_color_viridis_c(option = "C") + # geom_ribbon(alpha=0.2) + # scale_x_reverse(limits = c(210,20)) + # ylim(5, 15) + # coord_flip() + # gfplot::theme_pbs() + # theme(legend.position = c(.8,.4))
p_nd4 <- p_nd %>% filter(ssid==4) pred_4 <- plot_map_raster(p_nd4, "est") + scale_fill_viridis_c(labels = labels, option = "C", limits = c(5, 11), na.value = "black")+ facet_wrap(~year) + theme(legend.position = c(.85,.12)) + labs(fill = "°C") pred_4
plot_map_raster(p_nd, "est_non_rf") + ggtitle("Prediction (fixed effects only)") + facet_wrap(~year) + scale_fill_viridis_c()
plot_map_raster(p_nd, "est_rf") + ggtitle("Prediction (all random effects only)") + facet_wrap(~year) + scale_fill_viridis_c()
(omega_s_all <- plot_map_raster(p_nd, "omega_s") + ggtitle("Spatial random effects only") + scale_fill_gradient2( low = "Steel Blue 4", high = "Red 3", limits = c(-1.5,1.5), na.value = "red") + labs(fill = "") + theme(legend.position = c(0.2, 0.3)) ) ggsave(filename = here::here("analysis/tmb-sensor-explore/figs/omegas-temp-rw-800-kn.png"), width = 3, height = 3) omega_s_all
plot_map_raster(p_nd, "epsilon_st") + ggtitle("Spatiotemporal random effects only") + facet_wrap(~year) + scale_fill_gradient2( low = "Steel Blue 4", high = "Red 3", limits = c(-1.5,1.5), na.value = "red")
p_nd <- readRDS("../VOCC/data/predicted_temp_allyears_800kn.rds") p_nd1 <- p_nd %>% filter(ssid!=16) %>% filter(ssid!=4) %>% filter(year>2007) (st_1n3 <- plot_map_raster(p_nd1, "epsilon_st") + ggtitle("Spatiotemporal random effects") + facet_grid(~year) + scale_fill_gradient2( low = "Steel Blue 4", high = "Red 3", limits = c(-1.5,1.5), #trans = sqrt, na.value = "red") + labs(fill = "epsilon") + theme( legend.position = "none")) ggsave(filename = here::here("analysis/tmb-sensor-explore/figs/temp-epsilon-1n3.png"), width =6, height = 2) # png( # file = "figs/temp-predictions-1n3.png", # res = 400, # units = "in", # width = 8, # height = 8 # ) # st_1n3 # dev.off()
p_nd <- readRDS("../VOCC/data/predicted_temp_allyears_800kn.rds") p_nd4 <- p_nd %>% filter(ssid==4) %>% filter(year>2007)#%>% filter(ssid!=4) (st_4 <- plot_map_raster(p_nd4, "epsilon_st") + #ggtitle(" ") + facet_grid(~year) + scale_fill_gradient2( low = "Steel Blue 4", high = "Red 3", limits = c(-1.5,1.5), #trans = sqrt, na.value = "red") + labs(fill = " ") + theme( legend.position = "none") ) ggsave(filename = here::here("analysis/tmb-sensor-explore/figs/temp-epsilon-4.png"), width = 6, height = 2) st_4
p_nd <- readRDS("../VOCC/data/predicted_temp_allyears_800kn.rds") p_nd16 <- p_nd %>% filter(ssid==16) %>% filter(year>2007)#%>% filter(ssid!=4) (st_16 <- plot_map_raster(p_nd16, "epsilon_st") + #ggtitle(" ") + facet_grid(~year) + scale_fill_gradient2( low = "Steel Blue 4", high = "Red 3", limits = c(-1.5,1.5), na.value = "red") + xlim(20,360) + ylim(5850,6050)+ # xlim(65,345) + ylim(5800,6100)+ theme(panel.border = element_blank(), legend.position = "none") ) ggsave(filename = here::here("analysis/tmb-sensor-explore/figs/temp-epsilon-16.png"), width = 6, height = 1) st_16
png( file = "figs/temp-predictions.png", res = 600, units = "in", width = 11, height = 10 ) gridExtra::grid.arrange( grobs = list(pred_1, st_1n3, omega_s_all, pred_4, st_4, st_16 ), widths = c(2, 2, 1.75), heights = c(2, 1.35), layout_matrix = rbind(c(1, 2, 3), c(4, 5, 6)), top = grid::textGrob("Bottom temperature estimates from sdmTMB model") ) # gridExtra::grid.arrange( nrow = 2, heights = c(2,1.5), # top = grid::textGrob("Dissolved oxygen estimates from sdmTMB model")) dev.off()
BOOSTED REGRESSION MODEL
library(gbm) m <- gbm( temperature_c ~ depth_scaled + X + Y, data = d_trawl, n.trees = 2000, interaction.depth = 3, shrinkage = 0.02) m2 <- gbm( temperature_c ~ SST_scaled + depth_scaled + X + Y, data = d_trawl, n.trees = 2000, interaction.depth = 3, shrinkage = 0.02) plot(m,i.var=2) plot(m,i.var=3) plot(m,i.var=4) plot(m,i.var=5) plot(m2,i.var=1) plot(m2,i.var=2) plot(m,i.var=3) plot(m,i.var=4) plot(m2,i.var=1:2) plot(m2,i.var=c(1,3)) plot(m2,i.var=c(1,4)) d_trawl$r <- predict(m, n.trees = 2000) - d_trawl1$temperature_c ggplot(d_trawl, aes(depth_scaled, r)) + geom_point(alpha=0.4) + ylim(-8,8) + geom_smooth()
CENTER OF GRAVITY
p_cog <- predict(m_temp_rw, newdata = nd_all, return_tmb_object = TRUE) temp_cog <- get_cog(p_cog) temp_cog %>% reshape2::dcast(year ~ coord, value.var = "est") %>% ggplot(aes(X, Y, colour = year)) + geom_path(arrow = arrow()) + scale_color_viridis_c()
temp_cog %>% ggplot(aes(year, est, ymin = lwr, ymax = upr)) + geom_line() + geom_ribbon(alpha = 0.2) + facet_wrap(~coord, scales = "free_y")
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.