library(tidyverse)
library(mnsentinellakes)
library(mgcv)
library(gratia)
library(mwlaxeref)
library(cowplot)
options(scipen = 999)
zoop <- read.csv("./data/ramsey_county_zooplankton.csv")
# These are calulated from all of white bear ice data (not just post-1980)
early.anomaly = -17.52
late.anomaly = 14
Combine ice data with zooplankton data
Use ice data only from White Bear (DOW = 82016700)
ice_off <- read.csv("./data/ice_off_summarized.csv") %>%
mutate(DNRID = fixlakeid(DOW)) %>%
dplyr::select(DNRID, year, min_ice_off_julian, min_ice_off_date) %>%
dplyr::filter(DNRID == "82016700")
mean_ice <- ice_off %>%
dplyr::filter(year < 1980) %>%
summarise(mean_ice_off = mean(min_ice_off_julian))
# Negative values indicate earlier anomalies
wb.ice = ice_off %>%
add_row(data.frame(year = 2024, min_ice_off_julian = 67, min_ice_off_date = "2024-03-08")) %>%
mutate(ice_off_anomaly = min_ice_off_julian - mean_ice$mean_ice_off) %>%
select(!DNRID) %>%
filter(year > 1980)
temp <- merge(wb.ice, zoop, by = c("year")) %>%
mutate(days_since_ice = JULIANDAY - min_ice_off_julian) %>%
dplyr::filter(!is.na(CYCLO.THOUS.M3)) %>%
mutate(ice_off_anomaly = as.numeric(ice_off_anomaly),
JULIANDAY = as.numeric(JULIANDAY),
year = as.numeric(year),
year_f = as.factor(year),
DNRID = as.factor(DNRID),
ice_off_anomaly_abs = abs(ice_off_anomaly)) %>%
# Subjective but remove early (pre-february) and late (post Decemeber) dates
dplyr::filter(JULIANDAY > 60) %>%
dplyr::filter(JULIANDAY < 330) %>%
rowwise() %>%
mutate('ZOOP.THOUS.M3' = sum(c_across(CYCLO.THOUS.M3:DIAPHAN.THOUS.M3))) %>%
select(year, min_ice_off_julian, min_ice_off_date, ice_off_anomaly,
DNRID, SITE, DATE, TOW, ZOOP.THOUS.M3, everything()) %>% ungroup() %>% as.data.frame()
sampling.events <- temp %>% group_by(DNRID, year) %>% count() #%>% filter(n > 200)
median(sampling.events$n)
## [1] 8
Run all GAMS
# loops through columns 10-20 in the dataframe temp and runs nb GAMMs
# TO DO - Play with offsets and counts
models <- list()
for (i in 10:18) { #10-18
single.model <- gam(temp[,i] ~ s(ice_off_anomaly) +
s(JULIANDAY) +
ti(ice_off_anomaly, JULIANDAY) +
s(DNRID, bs = "re") + s(year_f, bs = "re"),
method = "REML", select = T,
data = temp, family = nb())
print(paste0("Finished ", colnames(temp)[i]))
models[[paste0(colnames(temp)[i])]] <- single.model
}
## [1] "Finished CYCLO.THOUS.M3"
## [1] "Finished CALA.THOUS.M3"
## [1] "Finished NAUP.THOUS.M3"
## [1] "Finished ROTIF.THOUS.M3"
## [1] "Finished DAPH.THOUS.M3"
## [1] "Finished BOSM.THOUS.M3"
## [1] "Finished CHYD.THOUS.M3"
## [1] "Finished CERIOD.THOUS.M3"
## [1] "Finished DIAPHAN.THOUS.M3"
summary(models$CYCLO.THOUS.M3)
##
## Family: Negative Binomial(1.482)
## Link function: log
##
## Formula:
## temp[, i] ~ s(ice_off_anomaly) + s(JULIANDAY) + ti(ice_off_anomaly,
## JULIANDAY) + s(DNRID, bs = "re") + s(year_f, bs = "re")
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.372 0.178 13.33 <0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(ice_off_anomaly) 0.001762 9 0.0 0.965
## s(JULIANDAY) 7.798189 9 2387.2 <0.0000000000000002 ***
## ti(ice_off_anomaly,JULIANDAY) 14.416896 16 668.7 <0.0000000000000002 ***
## s(DNRID) 26.720666 27 10671.9 <0.0000000000000002 ***
## s(year_f) 42.568367 43 17829.9 <0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.287 Deviance explained = 58.9%
## -REML = 19583 Scale est. = 1 n = 5990
summary(models$CALA.THOUS.M3)
##
## Family: Negative Binomial(1.244)
## Link function: log
##
## Formula:
## temp[, i] ~ s(ice_off_anomaly) + s(JULIANDAY) + ti(ice_off_anomaly,
## JULIANDAY) + s(DNRID, bs = "re") + s(year_f, bs = "re")
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.6784 0.1467 11.44 <0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(ice_off_anomaly) 0.005409 9 0.02 0.65
## s(JULIANDAY) 7.173303 9 1970.56 <0.0000000000000002 ***
## ti(ice_off_anomaly,JULIANDAY) 10.558149 16 204.97 <0.0000000000000002 ***
## s(DNRID) 26.077445 27 1622.64 <0.0000000000000002 ***
## s(year_f) 42.424057 43 7486.51 <0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.287 Deviance explained = 42.6%
## -REML = 16278 Scale est. = 1 n = 5990
summary(models$NAUP.THOUS.M3)
##
## Family: Negative Binomial(1.835)
## Link function: log
##
## Formula:
## temp[, i] ~ s(ice_off_anomaly) + s(JULIANDAY) + ti(ice_off_anomaly,
## JULIANDAY) + s(DNRID, bs = "re") + s(year_f, bs = "re")
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.1067 0.1708 18.19 <0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(ice_off_anomaly) 0.001765 9 0.0 0.989
## s(JULIANDAY) 7.502667 9 2756.9 < 0.0000000000000002 ***
## ti(ice_off_anomaly,JULIANDAY) 8.348112 16 75.5 0.0000343 ***
## s(DNRID) 26.744701 27 16289.4 < 0.0000000000000002 ***
## s(year_f) 42.680855 43 28487.7 < 0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.233 Deviance explained = 59.4%
## -REML = 23394 Scale est. = 1 n = 5990
summary(models$ROTIF.THOUS.M3)
##
## Family: Negative Binomial(0.544)
## Link function: log
##
## Formula:
## temp[, i] ~ s(ice_off_anomaly) + s(JULIANDAY) + ti(ice_off_anomaly,
## JULIANDAY) + s(DNRID, bs = "re") + s(year_f, bs = "re")
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.4417 0.2483 9.835 <0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(ice_off_anomaly) 0.001706 9 0.003 0.803
## s(JULIANDAY) 7.469614 9 940.899 < 0.0000000000000002
## ti(ice_off_anomaly,JULIANDAY) 9.328772 16 90.890 0.0000105
## s(DNRID) 26.380036 27 6147.263 < 0.0000000000000002
## s(year_f) 42.606172 43 18818.082 < 0.0000000000000002
##
## s(ice_off_anomaly)
## s(JULIANDAY) ***
## ti(ice_off_anomaly,JULIANDAY) ***
## s(DNRID) ***
## s(year_f) ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.065 Deviance explained = 46.7%
## -REML = 19053 Scale est. = 1 n = 5990
summary(models$DAPH.THOUS.M3)
##
## Family: Negative Binomial(0.832)
## Link function: log
##
## Formula:
## temp[, i] ~ s(ice_off_anomaly) + s(JULIANDAY) + ti(ice_off_anomaly,
## JULIANDAY) + s(DNRID, bs = "re") + s(year_f, bs = "re")
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.4713 0.1572 9.361 <0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(ice_off_anomaly) 0.4554 9 603.3 0.178
## s(JULIANDAY) 7.5484 9 1606.5 <0.0000000000000002 ***
## ti(ice_off_anomaly,JULIANDAY) 12.6630 16 418.9 <0.0000000000000002 ***
## s(DNRID) 26.2236 27 1365.0 <0.0000000000000002 ***
## s(year_f) 41.6616 43 4571.9 <0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.283 Deviance explained = 37.9%
## -REML = 15347 Scale est. = 1 n = 5990
summary(models$BOSM.THOUS.M3)
##
## Family: Negative Binomial(0.437)
## Link function: log
##
## Formula:
## temp[, i] ~ s(ice_off_anomaly) + s(JULIANDAY) + ti(ice_off_anomaly,
## JULIANDAY) + s(DNRID, bs = "re") + s(year_f, bs = "re")
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.0989 0.3403 6.169 0.000000000689 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(ice_off_anomaly) 0.00232 9 0.016 0.188
## s(JULIANDAY) 7.62130 9 1744.395 <0.0000000000000002 ***
## ti(ice_off_anomaly,JULIANDAY) 12.31899 16 355.786 <0.0000000000000002 ***
## s(DNRID) 26.84852 27 10673.965 <0.0000000000000002 ***
## s(year_f) 42.02908 43 14766.037 <0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.15 Deviance explained = 52.7%
## -REML = 17174 Scale est. = 1 n = 5990
summary(models$CHYD.THOUS.M3)
##
## Family: Negative Binomial(0.334)
## Link function: log
##
## Formula:
## temp[, i] ~ s(ice_off_anomaly) + s(JULIANDAY) + ti(ice_off_anomaly,
## JULIANDAY) + s(DNRID, bs = "re") + s(year_f, bs = "re")
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.2556 0.3278 0.78 0.436
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(ice_off_anomaly) 0.002075 9 0.004 0.493
## s(JULIANDAY) 6.820935 9 857.528 <0.0000000000000002
## ti(ice_off_anomaly,JULIANDAY) 6.120625 16 106.526 <0.0000000000000002
## s(DNRID) 26.614806 27 4909.712 <0.0000000000000002
## s(year_f) 41.951927 43 11242.338 <0.0000000000000002
##
## s(ice_off_anomaly)
## s(JULIANDAY) ***
## ti(ice_off_anomaly,JULIANDAY) ***
## s(DNRID) ***
## s(year_f) ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0582 Deviance explained = 48%
## -REML = 9548.8 Scale est. = 1 n = 5990
summary(models$CERIOD.THOUS.M3)
##
## Family: Negative Binomial(0.213)
## Link function: log
##
## Formula:
## temp[, i] ~ s(ice_off_anomaly) + s(JULIANDAY) + ti(ice_off_anomaly,
## JULIANDAY) + s(DNRID, bs = "re") + s(year_f, bs = "re")
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.8868 0.4878 -1.818 0.0691 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(ice_off_anomaly) 0.005018 9 0.001 0.812
## s(JULIANDAY) 6.940875 9 2149.233 <0.0000000000000002 ***
## ti(ice_off_anomaly,JULIANDAY) 7.719588 16 132.333 <0.0000000000000002 ***
## s(DNRID) 26.674936 27 4367.261 <0.0000000000000002 ***
## s(year_f) 41.221625 43 9562.052 <0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = -0.951 Deviance explained = 61.8%
## -REML = 6651.2 Scale est. = 1 n = 5990
appraise(models$CYCLO.THOUS.M3)
appraise(models$CALA.THOUS.M3)
appraise(models$NAUP.THOUS.M3)
appraise(models$ROTIF.THOUS.M3)
appraise(models$DAPH.THOUS.M3)
appraise(models$BOSM.THOUS.M3)
appraise(models$CHYD.THOUS.M3)
appraise(models$CERIOD.THOUS.M3)
k.check(models$CYCLO.THOUS.M3)
## k' edf k-index p-value
## s(ice_off_anomaly) 9 0.001761511 0.9056888 0.0550
## s(JULIANDAY) 9 7.798189230 0.9211881 0.3125
## ti(ice_off_anomaly,JULIANDAY) 16 14.416895885 0.8098815 0.0000
## s(DNRID) 28 26.720666299 NA NA
## s(year_f) 44 42.568366803 NA NA
k.check(models$CALA.THOUS.M3)
## k' edf k-index p-value
## s(ice_off_anomaly) 9 0.005408851 0.9186060 0.9125
## s(JULIANDAY) 9 7.173302533 0.9025275 0.5475
## ti(ice_off_anomaly,JULIANDAY) 16 10.558148938 0.8433064 0.0000
## s(DNRID) 28 26.077445209 NA NA
## s(year_f) 44 42.424057243 NA NA
k.check(models$NAUP.THOUS.M3)
## k' edf k-index p-value
## s(ice_off_anomaly) 9 0.001764526 0.9296690 0.3050
## s(JULIANDAY) 9 7.502666963 0.9289166 0.2475
## ti(ice_off_anomaly,JULIANDAY) 16 8.348112192 0.8387346 0.0000
## s(DNRID) 28 26.744700867 NA NA
## s(year_f) 44 42.680855294 NA NA
k.check(models$ROTIF.THOUS.M3)
## k' edf k-index p-value
## s(ice_off_anomaly) 9 0.001706475 0.8018154 0.425
## s(JULIANDAY) 9 7.469613965 0.7688942 0.000
## ti(ice_off_anomaly,JULIANDAY) 16 9.328772372 0.7267178 0.000
## s(DNRID) 28 26.380035520 NA NA
## s(year_f) 44 42.606172435 NA NA
k.check(models$DAPH.THOUS.M3)
## k' edf k-index p-value
## s(ice_off_anomaly) 9 0.4553822 0.8728431 0.3700
## s(JULIANDAY) 9 7.5484387 0.8706399 0.2975
## ti(ice_off_anomaly,JULIANDAY) 16 12.6629806 0.8288372 0.0000
## s(DNRID) 28 26.2236144 NA NA
## s(year_f) 44 41.6615817 NA NA
k.check(models$BOSM.THOUS.M3)
## k' edf k-index p-value
## s(ice_off_anomaly) 9 0.002320306 0.7887612 0.4725
## s(JULIANDAY) 9 7.621298026 0.7925971 0.5625
## ti(ice_off_anomaly,JULIANDAY) 16 12.318992891 0.7481422 0.0000
## s(DNRID) 28 26.848520122 NA NA
## s(year_f) 44 42.029079770 NA NA
k.check(models$CHYD.THOUS.M3)
## k' edf k-index p-value
## s(ice_off_anomaly) 9 0.00207526 0.7522163 0.1225
## s(JULIANDAY) 9 6.82093486 0.7785335 0.9425
## ti(ice_off_anomaly,JULIANDAY) 16 6.12062519 0.7295041 0.0000
## s(DNRID) 28 26.61480559 NA NA
## s(year_f) 44 41.95192713 NA NA
k.check(models$CERIOD.THOUS.M3)
## k' edf k-index p-value
## s(ice_off_anomaly) 9 0.005017647 0.7613685 0.5925
## s(JULIANDAY) 9 6.940874642 0.7431696 0.0875
## ti(ice_off_anomaly,JULIANDAY) 16 7.719588431 0.7194170 0.0000
## s(DNRID) 28 26.674935555 NA NA
## s(year_f) 44 41.221625156 NA NA
write_rds(models, "./models/zoop.models.rds")
Read-in outputted models
models <- read_rds("./models/zoop.models.rds")
DOY change function
# Ouputs a figure for changes in zoop dynamics across a year
plot.zoop.doy = function(model = NA, var.name = NA, y.label = NA,
early.anomaly = -17.52, late.anomaly = 14, start.date = 80, end.date = 274) {
if (is.na(model)) {
model <- models[[var.name]]
}
# Year and DNRID aren't used but are required for predict.gam
length <- length(seq(start.date, end.date, by = 1))
pred.df <- data.frame(ice_off_anomaly_abs =
c(rep(early.anomaly, length), rep(0, length), rep(late.anomaly, length)),
ice_off_anomaly =
c(rep(early.anomaly, length), rep(0, length), rep(late.anomaly, length)),
JULIANDAY = seq(start.date, end.date, by = 1),
DNRID = "62001200", year_f = "2000")
fit <- predict.gam(model, pred.df,
unconditional = T, # Gives simultaneous
type = "response",
exclude = c('s(DNRID)',
's(year_f)'),
se.fit = T)
pred.df$predictions = fit$fit
pred.df$se = fit$se.fit
y.label <- species.labels %>% filter(variable == var.name) %>%
select(axis.label) %>% as.character()
plot <- pred.df %>%
mutate(ice_off_anomaly = as.factor(ice_off_anomaly),
ice_off_f =
ifelse(ice_off_anomaly == early.anomaly, "Early", "Normal"),
ice_off_f =
ifelse(ice_off_anomaly == late.anomaly, "Late", ice_off_f),
ice_off_f = factor(ice_off_f,
levels = c("Early", "Normal", "Late"))) %>%
ggplot(aes(y = predictions, x = as.Date(JULIANDAY, origin = as.Date("2018-01-01")),
color = ice_off_anomaly)) +
geom_line(aes(y = predictions, x = as.Date(JULIANDAY, origin = as.Date("2018-01-01")),
color = ice_off_f), lwd = 2) +
theme_classic() +
scale_x_date(date_breaks = "months" , date_labels = "%b") +
scale_color_manual(values = c("#FF0000", "#000000", "#0000FF")) +
labs(y = paste0(y.label), x = "Julian Day",
color = "Ice-Off", fill = "Ice-Off") +
theme(axis.text = element_text(size = 14),
axis.title = element_text(size = 18),
legend.title = element_text(size = 16),
legend.text = element_text(size = 16))
return(plot)
}
Bayesian peak timing - FIGURE S2 (and Code for other figures)
toy.df <- data.frame(ice_off_anomaly = c(rep(early.anomaly, 201), rep(0, 201),
rep(late.anomaly, 201)),
JULIANDAY = seq(80, 280, by = 1),
DNRID = "82016700", year_f = "2000") %>%
mutate(.row = row_number()) %>%
mutate(ice_off_f = ifelse(ice_off_anomaly == early.anomaly, "Early", "Normal"),
ice_off_f = ifelse(ice_off_anomaly == late.anomaly, "Late", ice_off_f))
# Models to loop through
iters <- names(models)
timing.cis <- data.frame(group = character(0),
numeric = numeric(0),
timing = character(0),
median = numeric(0),
lower.95 = numeric(0),
upper.95 = numeric(0),
lower.75 = numeric(0),
upper.75 = numeric(0))
pred.df = toy.df
for (i in 1:length(iters)) {
model = models[[iters[i]]]
if (iters[i] %in% c("BOSM.THOUS.M3")) {
pred.df = pred.df %>% filter(JULIANDAY < 182)
}
pred.df <- pred.df %>% mutate(.row = row_number())
# Use posterior dist of the mean to get posterior fitted values
sims <- fitted_samples(model, exclude = c('s(DNRID)', 's(year_f)'),
data = pred.df, n = 1000, seed = 1024) |>
left_join(pred.df |> select(.row, ice_off_anomaly, JULIANDAY, ice_off_f),
by = join_by(.row == .row))
# Get the peak timing when the ice off anomaly is big (ice-off happens late)
early <- sims %>% filter(ice_off_anomaly == early.anomaly) %>%
group_by(.draw) %>% slice_max(.fitted)
early.quant <- quantile(early$JULIANDAY, probs = c(.025,.975))
early.quant.75 <- quantile(early$JULIANDAY, probs = c(.125, .875))
early.median <- median(early$JULIANDAY)
# Get the peak timing when the ice off anomaly is big (ice-off happens late)
late <- sims %>% filter(ice_off_anomaly == late.anomaly) %>%
group_by(.draw) %>% slice_max(.fitted)
late.quant <- quantile(late$JULIANDAY, probs = c(.025,.975))
late.quant.75 <- quantile(late$JULIANDAY, probs = c(.125, .875))
late.median <- median(late$JULIANDAY)
# Get the peak timing when the ice off anomaly is normal (ice-off is average)
average <- sims %>% filter(ice_off_anomaly == 0) %>%
group_by(.draw) %>% slice_max(.fitted)
average.quant <- quantile(average$JULIANDAY, probs = c(.025,.975))
average.quant.75 <- quantile(average$JULIANDAY, probs = c(.125, .875))
average.median <- median(average$JULIANDAY)
# Get the difference between peak timings
diff.early <- early$JULIANDAY - average$JULIANDAY
diff.earlymedian <- median(diff.early)
diff.earlyquant <- quantile(diff.early, probs = c(.025,.975))
diff.earlyquant.75 <- quantile(diff.early, probs = c(.125, .875))
diff.late <- late$JULIANDAY - average$JULIANDAY
diff.latemedian <- median(diff.late)
diff.latequant <- quantile(diff.late, probs = c(.025,.975))
diff.latequant.75 <- quantile(diff.late, probs = c(.125, .875))
# Make new data frame to rbind to big data frame
df <- data.frame(group = iters[[i]],
numeric = i,
timing = c("Early", "Late", "Average", "Early Difference", "Late Difference"),
median = c(early.median, late.median, average.median,
diff.earlymedian, diff.latemedian),
lower.95 = c(early.quant[1], late.quant[1], average.quant[1],
diff.earlyquant[1], diff.latequant[1]),
upper.95 = c(early.quant[2], late.quant[2], average.quant[2],
diff.earlyquant[2], diff.latequant[2]),
lower.75 = c(early.quant.75[1], late.quant.75[1], average.quant.75[1],
diff.earlyquant.75[1], diff.latequant.75[1]),
upper.75 = c(early.quant.75[2], late.quant.75[2], average.quant.75[2],
diff.earlyquant.75[2], diff.latequant.75[2]))
timing.cis <- rbind(timing.cis, df)
pred.df = toy.df
}
row.names(timing.cis) <- NULL
species.labels <- data.frame(variable = names(models),
axis.label = c("Cyclopoid Peak", "Calanoid Peak",
"Nauplii Peak", "Rotifer Peak", "Daphnia Peak",
"Bosmina Peak", "Chydoridae Peak", "Ceriodaphnia Peak",
"Diaphanosoma Peak"))
timing.cis.plot <- timing.cis %>%
mutate(jit = ifelse(timing == "Early", -0.15, 0.15),
jit = ifelse(timing == "Average", 0, jit),
jit = ifelse(timing == "Early Difference", -0.15, jit)) %>%
merge(y = species.labels, by.x = "group", by.y = "variable") %>%
mutate(axis.label = factor(axis.label,
levels = c(
"Cyclopoid", "Calanoid",
"Nauplii", "Rotifer", "Daphnia", "Bosmina",
"Chydoridae", "Ceriodaphnia", "Diaphanosoma")))
timing.cis.plot$timing <- factor(timing.cis.plot$timing,
levels = c("Early", "Average", "Late",
"Early Difference", "Late Difference"))
timing.cis.plot %>%
filter(!grepl("Difference", timing)) %>%
ggplot() +
geom_point(aes(x = as.Date(median, origin = as.Date("2018-01-01")), y = -1*(numeric + jit),
color = timing), size = 2) +
geom_segment(aes(y = -1*(numeric + jit), yend = -1*(numeric + jit),
x = as.Date(lower.95, origin = as.Date("2018-01-01")),
xend = as.Date(upper.95, origin = as.Date("2018-01-01")), color = timing)) +
geom_segment(aes(y = -1*(numeric + jit), yend = -1*(numeric + jit),
x = as.Date(lower.75, origin = as.Date("2018-01-01")),
xend = as.Date(upper.75, origin = as.Date("2018-01-01")), color = timing),
linewidth = 1.5) +
#facet_grid(group ~., switch = "y") +
theme_classic() +
labs(y = "",
x = "Day Peak Density Occurs", color = "Ice-off Anomaly") +
theme(panel.grid.major.x = element_line(),
panel.grid.minor.x = element_line(),
strip.text.y.left = element_text(angle = 0, size = 12, face = 2),
plot.title = element_text(angle = 0, size = 14, face = 2),
strip.placement = "outside",
strip.background = element_blank(),
axis.text.y = element_text(size = 14),
axis.text.x = element_text(size = 14),
axis.title.x = element_text(size = 16),
legend.title = element_text(size = 14),
legend.text = element_text(size = 12)
) +
scale_color_manual(values = c("#FF0000", "#000000", "#0000FF")) +
scale_y_continuous(breaks = -1:-9, labels = levels(timing.cis.plot$axis.label)) +
scale_x_date(date_labels = "%b")
#scale_x_continuous(breaks = seq(100, 275, by = 25))
# Figure S2
ggsave("./figures/FigureS2.jpeg", height = 10, width = 10)
FIGURE 2A - Bayesian peak timing - SUBSET - Calendar DOY
#read in more data
timing.cis.wae <- read.csv("./data/wae_spawning_modeled.csv")
timing.cis.phyto <- read.csv("./data/phyto_modeled_peak.csv")
# y-axis plot labels
subset.labels = c("Chrysophyte","Dinoflagellate",
"Diatom",
"Cyclopoid", "Calanoid",
"Daphnia", "Walleye \n Peak Spawning")
# groups to keep while filtering
groups_to_keep <- c("CHRYSOS", "DINOS","DIATOM",
"CALA.THOUS.M3", "CYCLO.THOUS.M3", "DAPH.THOUS.M3",
"peak")
# x-axis plot labels
labels = c("April 1st", "May 1st",
"June 1st", "July 1st")
# bind data together
timing.cis.plot.all <- rbind(timing.cis.plot, timing.cis.wae, timing.cis.phyto) %>%
filter(group != "start")
plankton_subset_wae <- timing.cis.plot.all %>%
filter(!grepl("Difference", timing)) %>%
filter(group %in% groups_to_keep) %>%
mutate(jit = ifelse(timing == "Early", -0.15, 0.15),
jit = ifelse(timing == "Average", 0, jit),
jit = ifelse(timing == "Early Difference", -0.15, jit)) %>%
mutate(numeric = ifelse(group == "CHRYSOS", 1, numeric),
numeric = ifelse(group == "DINOS", 2, numeric),
numeric = ifelse(group == "DIATOMS", 3, numeric),
numeric = ifelse(group == "CYCLO.THOUS.M3", 4, numeric),
numeric = ifelse(group == "CALA.THOUS.M3", 5, numeric),
numeric = ifelse(group == "DAPH.THOUS.M3", 6, numeric),
numeric = ifelse(group == "peak", 7, numeric)) %>%
ggplot() +
geom_point(aes(x = as.Date(median, origin = as.Date("2018-01-01")), y = -1*(numeric + jit),
color = timing), size = 2) +
geom_segment(aes(y = -1*(numeric + jit), yend = -1*(numeric + jit),
x = as.Date(lower.95, origin = as.Date("2018-01-01")),
xend = as.Date(upper.95, origin = as.Date("2018-01-01")), color = timing)) +
geom_segment(aes(y = -1*(numeric + jit), yend = -1*(numeric + jit),
x = as.Date(lower.75, origin = as.Date("2018-01-01")),
xend = as.Date(upper.75, origin = as.Date("2018-01-01")), color = timing),
linewidth = 1.5) +
#facet_grid(group ~., switch = "y") +
theme_classic() +
labs(y = "",
x = "Day Peak Density Occurs\n ", color = "Ice-off Anomaly") +
theme(panel.grid.major.x = element_line(),
panel.grid.minor.x = element_line(),
strip.text.y.left = element_text(angle = 0, size = 12, face = 2),
plot.title = element_text(angle = 0, size = 14, face = 2),
strip.placement = "outside",
strip.background = element_blank(),
axis.text.y = element_text(size = 14),
axis.text.x = element_text(size = 14),
axis.title.x = element_text(size = 16),
legend.title = element_text(size = 14),
legend.text = element_text(size = 12)
) +
theme(legend.position = c(.8, .1)) +
scale_color_manual(values = c("#FF0000", "#000000", "#0000FF")) +
scale_y_continuous(breaks = -1:-7, labels = subset.labels) +
scale_x_date(breaks = as.Date(c("2018-04-01", "2018-05-01",
"2018-06-01", "2018-07-01")),
label = labels) +
theme(plot.margin = margin(0,0,0,1, "cm"))
plankton_subset_wae
# Figure 2A
ggsave("./figures/Figure2A.jpeg", height = 8, width = 10)
FIGURE 2B - Bayesian peak timing - SUBSET - rate of change
# y-axis plot labels
subset.labels = c("Chrysophyte","Dinoflagellate",
"Diatom",
"Cyclopoid", "Calanoid",
"Daphnia", "Walleye \n Peak Spawning")
# groups to keep while filtering
groups_to_keep <- c("CHRYSOS", "DINOS","DIATOM",
"CALA.THOUS.M3", "CYCLO.THOUS.M3", "DAPH.THOUS.M3",
"peak")
rate_of_change <- timing.cis.plot.all %>%
filter(grepl("Difference", timing)) %>%
filter(group %in% c(groups_to_keep)) %>%
mutate(jit = ifelse(timing == "Early", -0.15, 0.15),
jit = ifelse(timing == "Average", 0, jit),
jit = ifelse(timing == "Early Difference", -0.15, jit)) %>%
mutate(numeric = ifelse(group == "CHRYSOS", 1, numeric),
numeric = ifelse(group == "DINOS", 2, numeric),
numeric = ifelse(group == "DIATOMS", 3, numeric),
numeric = ifelse(group == "CYCLO.THOUS.M3", 4, numeric),
numeric = ifelse(group == "CALA.THOUS.M3", 5, numeric),
numeric = ifelse(group == "DAPH.THOUS.M3", 6, numeric),
numeric = ifelse(group == "peak", 7, numeric)) %>%
mutate(timing = ifelse(timing == "Early Difference",
"Difference between early \nand average ice off",
"Difference between late \nand average ice off")) %>%
mutate(median = ifelse(timing == "Difference between early \nand average ice off",
(median/(early.anomaly))*early.anomaly, (median/late.anomaly)*late.anomaly),
upper.95 = ifelse(timing == "Difference between early \nand average ice off",
(upper.95/(early.anomaly))*early.anomaly, (upper.95/late.anomaly)*late.anomaly),
lower.95 = ifelse(timing == "Difference between early \nand average ice off",
(lower.95/(early.anomaly))*early.anomaly, (lower.95/late.anomaly)*late.anomaly),
upper.75 = ifelse(timing == "Difference between early \nand average ice off",
(upper.75/(early.anomaly))*early.anomaly, (upper.75/late.anomaly)*late.anomaly),
lower.75 = ifelse(timing == "Difference between early \nand average ice off",
(lower.75/(early.anomaly))*early.anomaly, (lower.75/late.anomaly)*late.anomaly)) %>%
ggplot() +
geom_point(aes(x = median, y = -1*(numeric + jit),
color = timing), size = 2) +
geom_segment(aes(y = -1*(numeric + jit), yend = -1*(numeric + jit),
x = lower.95,
xend = upper.95, color = timing)) +
geom_segment(aes(y = -1*(numeric + jit), yend = -1*(numeric + jit),
x = lower.75, xend = upper.75, color = timing), linewidth = 1.5) +
geom_vline(xintercept = 0, color = "black", linewidth = 1) +
geom_vline(xintercept = early.anomaly, color = "grey", linewidth = .75) +
geom_vline(xintercept = late.anomaly, color = "grey", linewidth = .75) +
theme_classic() +
labs(y = "",
x = "Change in Days Event Occurs Per\n Day Change in Ice-Off",
color = "Ice-off Anomaly") +
theme(panel.grid.major.x = element_line(),
panel.grid.minor.x = element_line(),
strip.text.y.left = element_text(angle = 0, size = 12, face = 2),
plot.title = element_text(angle = 0, size = 14, face = 2),
strip.placement = "outside",
strip.background = element_blank(),
axis.text.y = element_text(size = 14),
axis.text.x = element_text(size = 14),
axis.title.x = element_text(size = 16),
legend.title = element_text(size = 14),
legend.text = element_text(size = 12)
) +
scale_color_manual(values = c("#FF0000", "#0000FF")) +
scale_y_continuous(breaks = -1:-7, labels = subset.labels) + theme(legend.position="none")
rate_of_change <- timing.cis.plot.all %>%
filter(grepl("Difference", timing)) %>%
filter(group %in% c(groups_to_keep)) %>%
mutate(jit = ifelse(timing == "Early", -0.15, 0.15),
jit = ifelse(timing == "Average", 0, jit),
jit = ifelse(timing == "Early Difference", -0.15, jit)) %>%
mutate(numeric = ifelse(group == "CHRYSOS", 1, numeric),
numeric = ifelse(group == "DINOS", 2, numeric),
numeric = ifelse(group == "DIATOMS", 3, numeric),
numeric = ifelse(group == "CYCLO.THOUS.M3", 4, numeric),
numeric = ifelse(group == "CALA.THOUS.M3", 5, numeric),
numeric = ifelse(group == "DAPH.THOUS.M3", 6, numeric),
numeric = ifelse(group == "peak", 7, numeric)) %>%
mutate(timing = ifelse(timing == "Early Difference",
"Difference between early \nand average ice off",
"Difference between late \nand average ice off")) %>%
mutate(median = ifelse(timing == "Difference between early \nand average ice off",
median/(-1*early.anomaly), (median/late.anomaly)),
upper.95 = ifelse(timing == "Difference between early \nand average ice off",
upper.95/(-1*early.anomaly), upper.95/late.anomaly),
lower.95 = ifelse(timing == "Difference between early \nand average ice off",
lower.95/(-1*early.anomaly), lower.95/late.anomaly),
upper.75 = ifelse(timing == "Difference between early \nand average ice off",
upper.75/(-1*early.anomaly), upper.75/late.anomaly),
lower.75 = ifelse(timing == "Difference between early \nand average ice off",
lower.75/(-1*early.anomaly), lower.75/late.anomaly)) %>%
ggplot() +
geom_point(aes(x = median, y = -1*(numeric + jit),
color = timing), size = 2) +
geom_segment(aes(y = -1*(numeric + jit), yend = -1*(numeric + jit),
x = lower.95,
xend = upper.95, color = timing)) +
geom_segment(aes(y = -1*(numeric + jit), yend = -1*(numeric + jit),
x = lower.75, xend = upper.75, color = timing), linewidth = 1.5) +
geom_vline(xintercept = 0, color = "black", linewidth = 1) +
geom_vline(xintercept = -1, color = "grey", linewidth = .75) +
geom_vline(xintercept = 1, color = "grey", linewidth = .75) +
theme_classic() +
labs(y = "",
x = "Change in Days Event Occurs Per\n Day Change in Ice-Off",
color = "Ice-off Anomaly") +
theme(panel.grid.major.x = element_line(),
panel.grid.minor.x = element_line(),
strip.text.y.left = element_text(angle = 0, size = 12, face = 2),
plot.title = element_text(angle = 0, size = 14, face = 2),
strip.placement = "outside",
strip.background = element_blank(),
axis.text.y = element_text(size = 14),
axis.text.x = element_text(size = 14),
axis.title.x = element_text(size = 16),
legend.title = element_text(size = 14),
legend.text = element_text(size = 12)
) +
scale_color_manual(values = c("#FF0000", "#0000FF")) +
scale_y_continuous(breaks = -1:-7, labels = subset.labels) + theme(legend.position="none")
rate_of_change
# Figure 2B
ggsave("./figures/Figure2B_alt.jpeg", height = 8, width = 9)
SETUP FIGURE 4 - Recreating Feiner et al. 2022 figure
green_mean_ice_off <- 108.5
feiner.df <- data.frame(ice_off_anomaly = numeric(0),
JULIANDAY = numeric(0),
DNRID = numeric(0), year_f = numeric(0),
name = character(0), count = numeric(0)
)
anomaly = c(early.anomaly, 0, late.anomaly)
for (i in 1:3) {
pred.df <- data.frame(ice_off_anomaly = anomaly[i],
JULIANDAY = seq(80, 280, by = 1),
DNRID = "62001200", year_f = "2000")
fit <- predict.gam(models$CYCLO.THOUS.M3, pred.df,
unconditional = T, # Gives simultaneous
type = "response",
exclude = c('s(DNRID)', 's(year_f)'),
se.fit = T)
pred.df$Cyclopoid = fit$fit
fit.daph <- predict.gam(models$DAPH.THOUS.M3, pred.df,
unconditional = T, # Gives simultaneous
type = "response",
exclude = c('s(DNRID)', 's(year_f)'),
se.fit = T)
pred.df$Daphnia = fit.daph$fit
pred.df.long <- pred.df %>% pivot_longer(cols = Cyclopoid:Daphnia, values_to = "count") %>%
mutate(JULIANDAY = as.Date(JULIANDAY, origin = as.Date("2018-01-01")))
feiner.df = rbind(feiner.df, pred.df.long)
}
peaks <- feiner.df %>% group_by(ice_off_anomaly, name) %>% slice_max(count) %>%
mutate(JULIANDAY = as.Date(JULIANDAY, origin = as.Date("2018-01-01")))
# Make a theme for all plots
personal_theme = theme_classic() +
theme(axis.text.y = element_text(size = 14),
axis.text.x = element_text(size = 14),
axis.title.x = element_text(size = 16),
axis.title.y = element_text(size = 16),
legend.title = element_text(size = 14),
legend.text = element_text(size = 12)) +
theme(plot.title = element_text(size = 16, hjust = 0.5))
walleye_spawn_timing = data.frame(ice_off_anomaly = c(early.anomaly, 0, late.anomaly),
spawn = c(as.Date(106.0, origin = as.Date("2018-01-01")),
as.Date(114.0, origin = as.Date("2018-01-01")),
as.Date(125, origin = as.Date("2018-01-01"))))
# Make function to label plot
add_annotation <- function(group, ice_off) {
subset = peaks %>% filter(name == group) %>% filter(ice_off_anomaly == ice_off)
spawn_date = walleye_spawn_timing %>% filter(ice_off_anomaly == ice_off)
if(group == "Cyclopoid"){
text_location = spawn_date$spawn + (yday(subset$JULIANDAY)+1 - yday(spawn_date$spawn))/2
}
else{
text_location = subset$JULIANDAY - 9
}
# the function will return the last object it creates, ie this list with two objects
list(
annotate("segment", color = "black",
x = subset$JULIANDAY,
xend = spawn_date$spawn,
y = subset$count, yend = subset$count),
annotate("text", color = "black", size = 5,
x = text_location, y = subset$count + 1.5,
label = paste(yday(subset$JULIANDAY) - yday(spawn_date$spawn), "days")
)
)
}
FIGURE 4 - Recreating Feiner et al. 2022 figure
# Plot phenology when ice-off is early
early.anomaly.plot <- feiner.df %>%
filter(ice_off_anomaly == early.anomaly) %>%
ggplot() +
# Plot where ice-off and spawning occur
geom_vline(aes(xintercept = as.Date(green_mean_ice_off + early.anomaly, origin = as.Date("2018-01-01")),
color = "Ice-Off"), linewidth = 1.5) +
geom_vline(aes(xintercept = walleye_spawn_timing[1,2],
color = "Walleye Spawning"), linewidth = 1) +
# Plot where daphnia and cyclopoid peaks occur and label them
geom_segment(data = peaks %>% filter(name == "Cyclopoid" & ice_off_anomaly == early.anomaly),
aes(x = JULIANDAY, y = count,
xend = JULIANDAY, yend = -Inf),
color = "#A52A2A") +
geom_segment(data = peaks %>% filter(name == "Daphnia" & ice_off_anomaly == early.anomaly),
aes(x = JULIANDAY, y = count,
xend = JULIANDAY, yend = -Inf),
color = "#00FF00") +
add_annotation(group = "Cyclopoid", ice_off = early.anomaly) +
add_annotation(group = "Daphnia", ice_off = early.anomaly) +
# Change value with how much Walleye spawning gets earlier
geom_ribbon(aes(x = JULIANDAY, ymax = count, fill = name), ymin = 0, alpha = 0.3) +
coord_cartesian(xlim = as.Date(c("2018-03-31", "2018-10-15")),
ylim = c(0, 40)) +
scale_y_continuous(expand = c(0,0)) +
labs(title = "Early Ice-Off",
y = "Indivduals/L",
x = "Month", fill = "", color = "") +
scale_color_manual(guide = "legend",
breaks = c("Ice-Off", "Walleye Spawning"),
values = c("Ice-Off" = "Grey", "Walleye Spawning" = "Black")) +
scale_fill_manual(guide = "legend", breaks = c("Cyclopoid", "Daphnia"),
values = c("Cyclopoid" = "#A52A2A", "Daphnia" = "green")) +
scale_x_date(breaks = "1 month", date_labels = "%b") +
personal_theme + theme(legend.position = "none")
early.anomaly.plot
# Plot phenology when ice-off is Normal
average.anomaly.plot = feiner.df %>%
filter(ice_off_anomaly == 0) %>%
ggplot() +
geom_vline(aes(xintercept = as.Date(green_mean_ice_off, origin = as.Date("2018-01-01")),
color = "Ice-Off"), linewidth = 1.5) +
geom_vline(aes( xintercept = walleye_spawn_timing[2,2],
color = "Walleye Spawning"),
linewidth = 1) +
# Plot where daphnia and cyclopoid peaks occur and label them
geom_segment(data = peaks %>% filter(name == "Cyclopoid" & ice_off_anomaly == 0),
aes(x = JULIANDAY, y = count,
xend = JULIANDAY, yend = -Inf),
color = "#A52A2A") +
geom_segment(data = peaks %>% filter(name == "Daphnia" & ice_off_anomaly == 0),
aes(x = JULIANDAY, y = count,
xend = JULIANDAY, yend = -Inf),
color = "#00FF00") +
add_annotation(group = "Cyclopoid", ice_off = 0) +
add_annotation(group = "Daphnia", ice_off = 0) +
# Change value with how much Walleye spawning gets earlier
geom_ribbon(aes(x = JULIANDAY, ymax = count, fill = name),
ymin = 0, alpha = 0.3) +
coord_cartesian(xlim = as.Date(c("2018-03-31", "2018-10-15")), ylim = c(0, 40)) +
scale_y_continuous(expand = c(0,0)) +
labs(title = "Average Ice-Off",
y = "Indivduals/L",
x = "Month", fill = "", color = "") +
scale_color_manual(guide = "legend",
breaks = c("Ice-Off", "Walleye Spawning"),
values = c("Ice-Off" = "Grey", "Walleye Spawning" = "Black")) +
scale_fill_manual(guide = "legend", breaks = c("Cyclopoid", "Daphnia"),
values = c("Cyclopoid" = "#A52A2A", "Daphnia" = "green")) +
scale_x_date(breaks = "1 month", date_labels = "%b") +
personal_theme + theme(legend.position = c(.8, .7))
average.anomaly.plot
# Plot phenology when ice-off is late
late.anomaly.plot = feiner.df %>%
filter(ice_off_anomaly == late.anomaly) %>%
ggplot() +
geom_vline(aes(xintercept = as.Date(green_mean_ice_off + late.anomaly, origin = as.Date("2018-01-01")),
color = "Ice-Off"), linewidth = 1.5) +
geom_vline(aes(xintercept = walleye_spawn_timing[3,2],
color = "Walleye Spawning"), linewidth = 1) +
# Plot where daphnia and cyclopoid peaks occur and label them
geom_segment(data = peaks %>% filter(name == "Cyclopoid" & ice_off_anomaly == late.anomaly),
aes(x = JULIANDAY, y = count,
xend = JULIANDAY, yend = -Inf),
color = "#A52A2A") +
geom_segment(data = peaks %>% filter(name == "Daphnia" & ice_off_anomaly == late.anomaly),
aes(x = JULIANDAY, y = count,
xend = JULIANDAY, yend = -Inf),
color = "#00FF00") +
add_annotation(group = "Cyclopoid", ice_off = late.anomaly) +
add_annotation(group = "Daphnia", ice_off = late.anomaly) +
# Change value with how much Walleye spawning gets earlier
geom_ribbon(aes(x = JULIANDAY, ymax = count, fill = name),
ymin = 0, alpha = 0.3) +
coord_cartesian(xlim = as.Date(c("2018-03-31", "2018-10-15")), ylim = c(0, 40)) +
scale_y_continuous(expand = c(0,0)) +
labs(title = "Late Ice-Off",
y = "Indivduals/L",
x = "Month", fill = "", color = "") +
scale_color_manual(guide = "legend",
breaks = c("Ice-Off", "Walleye Spawning"),
values = c("Ice-Off" = "Grey", "Walleye Spawning" = "Black")) +
scale_fill_manual(guide = "legend", breaks = c("Cyclopoid", "Daphnia"),
values = c("Cyclopoid" = "Brown", "Daphnia" = "green")) +
scale_x_date(breaks = "1 month", date_labels = "%b") +
personal_theme + theme(legend.position = "none")
late.anomaly.plot
cowplot::plot_grid(early.anomaly.plot, average.anomaly.plot, late.anomaly.plot,
nrow = 3, ncol = 1, labels = "AUTO")
# Figure 4
ggsave("./figures/Figure4.jpeg", height = 10, width = 10)
Figure S5 - Bayesian peak timing - SUBSET - rate of change
subset.labels = c("Bosmina", "Calanoid", "Ceriodaphnia", "Chrysophyte", "Chydorus",
"Cryptomonad", "Cyanobacteria", "Cyclopoid","Daphnia", "Diaphanosoma",
"Diatom", "Dinoflagellate","Green-algae", "Nauplii", "Rotifera",
"Walleye \n Peak Spawning")
all_rate_of_change <- timing.cis.plot.all %>%
filter(grepl("Difference", timing)) %>%
filter(!group %in% c("ZOOP.THOUS.M3", "ALL.CELLS")) %>%
mutate(jit = ifelse(timing == "Early", -0.15, 0.15),
jit = ifelse(timing == "Average", 0, jit),
jit = ifelse(timing == "Early Difference", -0.15, jit),
numeric = as.numeric(factor(group)),
numeric = ifelse(group == "ROTIF.THOUS.M3", 15, numeric),
numeric = ifelse(group == "peak", 16, numeric)) %>%
mutate(timing = ifelse(timing == "Early Difference",
"Early ice-off",
"Late ice-off")) %>%
mutate(median = ifelse(timing == "Early ice-off",
median/(-1*early.anomaly), median/late.anomaly),
upper.95 = ifelse(timing == "Early ice-off",
upper.95/(-1*early.anomaly), upper.95/late.anomaly),
lower.95 = ifelse(timing == "Early ice-off",
lower.95/(-1*early.anomaly), lower.95/late.anomaly),
upper.75 = ifelse(timing == "Early ice-off",
upper.75/(-1*early.anomaly), upper.75/late.anomaly),
lower.75 = ifelse(timing == "Early ice-off",
lower.75/(-1*early.anomaly), lower.75/late.anomaly)) %>%
ggplot() +
geom_point(aes(x = median, y = -1*(numeric + jit),
color = timing), size = 2) +
geom_segment(aes(y = -1*(numeric + jit), yend = -1*(numeric + jit),
x = lower.95,
xend = upper.95, color = timing)) +
geom_segment(aes(y = -1*(numeric + jit), yend = -1*(numeric + jit),
x = lower.75, xend = upper.75, color = timing), linewidth = 1.5) +
geom_vline(xintercept = 0, color = "black", linewidth = 1) +
geom_vline(xintercept = -1, color = "grey", linewidth = .75) +
geom_vline(xintercept = 1, color = "grey", linewidth = .75) +
theme_classic() +
labs(y = "",
x = "Change in Days Event Occurs Per\n Day Change in Ice-Off",
color = "Ice-off Anomaly") +
theme(panel.grid.major.x = element_line(),
panel.grid.minor.x = element_line(),
strip.text.y.left = element_text(angle = 0, size = 12, face = 2),
plot.title = element_text(angle = 0, size = 14, face = 2),
strip.placement = "outside",
strip.background = element_blank(),
axis.text.y = element_text(size = 14),
axis.text.x = element_text(size = 14),
axis.title.x = element_text(size = 16),
legend.title = element_text(size = 14),
legend.text = element_text(size = 12)
) +
scale_color_manual(values = c("#FF0000", "#0000FF")) +
scale_y_continuous(breaks = -1:-16, labels = subset.labels) + theme(legend.position="none")
all_rate_of_change
# Figure S5
ggsave("./figures/FigureS5.jpeg", height = 8, width = 9)
See if zooplankton bloom size changes with ice-anomaly
Figure S7
species.labels <- data.frame(variable = names(models),
axis.label = c("Cyclopoid", "Calanoid",
"Nauplii", "Rotifer", "Daphnia",
"Bosmina", "Chydoridae", "Ceriodaphnia",
"Diaphanosoma"))
bloom.size <- timing.cis.plot.all %>%
filter(!timing %in% c("Early Difference", "Late Difference")) %>%
# filter out relative changes in timing
filter(!group %in% c("peak")) %>% # filter out peak walleye spawning, not relevant
filter(!group %in% c("CHRYSOS", "CRYPTOS", "CYANOS", "DIATOM", "DINOS", "GREEN")) %>%
mutate(.row = row_number(),
ice_off_anomaly = ifelse(timing == "Early", early.anomaly, 0),
ice_off_anomaly = ifelse(timing == "Late", late.anomaly, ice_off_anomaly))
bloom.size.plot <- data.frame(group = character(0), ice_off_anomaly = numeric(0),
median = numeric(0),
q2.5 = numeric(0), q97.5 = numeric(0),
q25 = numeric(0), q75 = numeric(0))
groups = unique(bloom.size$group)
for (i in 1:length(groups)) {
model = models[[groups[i]]]
group.df <- bloom.size %>% filter(group == groups[i])
pred.df <- data.frame(JULIANDAY =
c(
seq(group.df$median[1] - 3, group.df$median[1] + 3, by = 1 ),
seq(group.df$median[2] - 3, group.df$median[2] + 3, by = 1 ),
seq(group.df$median[3] - 3, group.df$median[3] + 3, by = 1 )),
ice_off_anomaly = c(rep(group.df$ice_off_anomaly[1], 7),
rep(group.df$ice_off_anomaly[2], 7),
rep(group.df$ice_off_anomaly[3], 7)),
DNRID = "82016700", year_f = "2000",
TOW = 6) %>%
mutate(.row = row_number())
sims <- fitted_samples(model, exclude = c('s(DNRID)', 's(year_f)'),
data = pred.df, n = 1000, seed = 1024) |>
left_join(pred.df |> select(.row, ice_off_anomaly, JULIANDAY),
by = join_by(.row == .row))
early.diff <- ((sims[sims$ice_off_anomaly == early.anomaly,]$.fitted -
sims[sims$ice_off_anomaly == 0,]$.fitted) /
sims[sims$ice_off_anomaly == 0,]$.fitted) * 100
late.diff <- ((sims[sims$ice_off_anomaly == late.anomaly,]$.fitted -
sims[sims$ice_off_anomaly == 0,]$.fitted) /
sims[sims$ice_off_anomaly == 0,]$.fitted) * 100
temp_toadd <- data.frame(group = groups[i],
ice_off_anomaly = c("Early", "Late"),
median = c(median(early.diff), median(late.diff)),
q2.5 = c(quantile(early.diff, prob = c(0.025)),
quantile(late.diff, prob = c(0.025))),
q97.5 = c(quantile(early.diff, prob = c(0.975)),
quantile(late.diff, prob = c(0.975))),
q25 = c(quantile(early.diff, prob = c(0.75)),
quantile(late.diff, prob = c(0.75))),
q75 = c(quantile(early.diff, prob = c(0.25)),
quantile(late.diff, prob = c(0.25))))
bloom.size.plot <- rbind(bloom.size.plot, temp_toadd)
}
bloom.size.plot %>%
merge(species.labels, by.y = "variable", by.x = "group") %>%
ggplot() +
geom_point(aes(y = median, x = ice_off_anomaly, color = ice_off_anomaly)) +
geom_segment(aes(y = q2.5, yend = q97.5, x = ice_off_anomaly,
xend = ice_off_anomaly, color = ice_off_anomaly)) +
geom_segment(aes(y = q25, yend = q75, x = ice_off_anomaly,
xend = ice_off_anomaly, color = ice_off_anomaly),
linewidth = 1.5) +
geom_hline(yintercept = 0, linetype = 'dotted') +
facet_wrap(vars(axis.label), scales = "free") +
scale_color_manual(values = c("#FF0000", "#0000FF")) +
labs(y = "% Change in Abundance From Normal Ice-Off",
x = "Ice-Off Timing", color = "Ice-Off Timing") +
theme_classic() +
theme(legend.position = c(.45, .935))
# Figure S7
ggsave("./figures/FigureS7.jpeg", height = 8, width = 9)
GRAVEYARD
NOT USED IN PAPER
DOY change function - BAYESIAN
plot.zoop.doy.bayes = function(model = NA, var.name = NA, type = "plot", slice = 1, y.label = NA,
early.anomaly = -17.52, late.anomaly = -14, start.date = 90, end.date = 280) {
if (is.na(model)) {
model <- models[[var.name]]
}
# Year and DNRID aren't used but are required for predict.gam
length <- length(seq(start.date, end.date, by = 1))
pred.df <- data.frame(ice_off_anomaly_abs =
c(rep(early.anomaly, length), rep(0, length), rep(late.anomaly, length)),
ice_off_anomaly =
c(rep(early.anomaly, length), rep(0, length), rep(late.anomaly, length)),
JULIANDAY = seq(start.date, end.date, by = 1),
days_since_ice = seq(1, length),
towdepth = 10,
DNRID = "62001200", year_f = "2000") %>%
mutate(.row = row_number()) %>%
mutate(
ice_off_f = ifelse(ice_off_anomaly == early.anomaly, "Early", "Normal"),
ice_off_f = ifelse(ice_off_anomaly == late.anomaly, "Late", ice_off_f))
fv <- fitted_values(model, data = pred.df, exclude = c('s(DNRID)', 's(year_f)'))
fs <- fitted_samples(model, data = pred.df, n = 20, seed = 1024,
exclude = c('s(DNRID)', 's(year_f)')) |>
left_join(pred.df |> select(.row, ice_off_anomaly, JULIANDAY),
by = join_by(.row == .row)) %>%
drop_na()
fit <- predict.gam(model, pred.df,
unconditional = T, # Gives simultaneous
type = "response",
exclude = c('s(DNRID)',
's(year_f)'),
se.fit = T)
pred.df$predictions = fit$fit
pred.df$se = fit$se.fit
if (type == "plot") {
y.label <- species.labels %>% filter(variable == var.name) %>%
select(axis.label) %>% as.character()
plot <- ggplot() +
#geom_line(data = fs,
# aes(group = draw,
# x = as.Date(JULIANDAY, origin = as.Date("2018-01-01")),
# y = fitted, colour = ice_off_f), alpha = 0.4) +
geom_line(data = pred.df, aes(y = predictions,
x = as.Date(JULIANDAY, origin = as.Date("2018-01-01")),
color = ice_off_f), lwd = 2) +
geom_ribbon(data = fv,
aes(x = as.Date(JULIANDAY, origin = as.Date("2018-01-01")),
y = .fitted, ymin = lower, ymax = upper,
fill = ice_off_f), alpha = 0.3) +
theme_classic() +
scale_x_date(date_breaks = "months" , date_labels = "%b") +
labs(y = paste0(y.label), x = "Julian Day",
color = "Ice-Off", fill = "Ice-Off") +
theme(axis.text = element_text(size = 14),
axis.title = element_text(size = 18),
legend.title = element_text(size = 16),
legend.text = element_text(size = 16))
return(plot)
} # End plot if statement
if (type == "max.date") {
max.date <- pred.df %>% group_by(ice_off_anomaly) %>% slice_max(predictions, n = slice) %>%
mutate(zoop = var.name)
max.date$predictions <- predict(model, newdata = max.date,
exclude = c("s(DNRID)", "s(year_f)"),
type = "response")
max.date$se <- predict(model, newdata = max.date, exclude = c("s(DNRID)", "s(year_f)"),
type = "response", se.fit = T)$se.fit
return(max.date)
} # End prediction if statement
}
Plot DOY change - CIs
hist(temp$JULIANDAY) # Approx 125-275 for peak sampling, look at 100 - 300
min(temp$JULIANDAY)
# Thoughts to do - Change x-axis to Months
species.labels <- data.frame(variable = names(models),
axis.label = c("Zooplankton Density","Cyclopoid Density", "Calanoid Density",
"Nauplii Density", "Rotifer Density",
"Daphnia Density", "Bosmina Density",
"Chydoridae Density", "Ceriodaphnia Density",
"Diaphanosoma Density"
))
cyclo.doy <- plot.zoop.doy.bayes(var.name = "CYCLO.THOUS.M3") + theme(legend.position="none")
#cyclo.doy
#ggsave("./figures/cyclo_doy.jpeg", height = 10, width = 10)
cala.doy <- plot.zoop.doy.bayes(var.name = "CALA.THOUS.M3") + theme(legend.position="none")
naup.doy <- plot.zoop.doy.bayes(var.name = "NAUP.THOUS.M3") + theme(legend.position="none")
rotif.doy <- plot.zoop.doy.bayes(var.name = "ROTIF.THOUS.M3") + theme(legend.position="none")
daph.doy <- plot.zoop.doy.bayes(var.name = "DAPH.THOUS.M3") + theme(legend.position="none")
bosm.doy <- plot.zoop.doy.bayes(var.name = "BOSM.THOUS.M3") + theme(legend.position="none")
ceriod.doy <- plot.zoop.doy.bayes(var.name = "CERIOD.THOUS.M3") + theme(legend.position="none")
diaphan.doy <- plot.zoop.doy.bayes(var.name = "DIAPHAN.THOUS.M3") +
theme(legend.position = c(0.25, 0.75))
zoop.doy <- plot.zoop.doy.bayes(var.name = "ZOOP.THOUS.M3") + theme(legend.position="none")
zoop <- cowplot::plot_grid(cyclo.doy, cala.doy, naup.doy, rotif.doy, daph.doy,
bosm.doy, ceriod.doy, diaphan.doy, zoop.doy)
title <- ggdraw() +
draw_label(
"Effects of Anomalous Ice-Off on Zooplankton Densities Over the Year",
fontface = 'bold',
size = 23,
x = 0,
hjust = 0
) +
theme(
# add margin on the left of the drawing canvas,
# so title is aligned with left edge of first plot
plot.margin = margin(0, 0, 0, 7)
)
plot_grid(
#title,
zoop,
ncol = 1,
# rel_heights values control vertical title margins
rel_heights = c(0.1, 1)
)
#ggsave("./figures/all_zoop_doy.bayes.jpeg", height = 10, width = 12)
Plot Anomaly results
plot.zoop.anomaly = function(var.name){
model <- models[[var.name]]
sm <- smooth_estimates(model)
sm <- sm %>%
add_constant(coef(model)["(Intercept)"]) %>%
add_confint() %>%
transform_fun(inv_link(model)) %>%
mutate(lower_ci = exp(lower_ci),
upper_ci = exp(upper_ci))
temp_residuals <- temp %>% add_partial_residuals(model) %>%
add_constant(coef(model)["(Intercept)"], column = 26) %>%
transform_fun(inv_link(model), column = 26)
# Creat upper limit for ggplot (1 or 2 really large residuals)
upper.limit <- quantile(temp_residuals$`s(ice_off_anomaly)`, 0.975)
fig <- sm %>% filter(smooth == "s(ice_off_anomaly)" ) %>%
ggplot(aes(x = ice_off_anomaly, y = est)) +
geom_point(aes(x = ice_off_anomaly, y = `s(ice_off_anomaly)`),
alpha = 0.1,
data = temp_residuals) +
geom_line(lwd = 1.5) +
geom_ribbon(aes(ymin = lower_ci , ymax = upper_ci), alpha = 0.3) +
ylim(0, upper.limit) +
theme_classic() +
labs(y = var.name, x = "Ice-Off Anomaly") +
theme(axis.text = element_text(size = 12),
axis.title = element_text(size = 16))
return(fig)
}
cyclo <- plot.zoop.anomaly("CYCLO.THOUS.M3")
cyclo
cala <- plot.zoop.anomaly("CALA.THOUS.M3")
cala
naup <- plot.zoop.anomaly("NAUP.THOUS.M3")
naup
rotif <- plot.zoop.anomaly("ROTIF.THOUS.M3")
rotif
daph <- plot.zoop.anomaly("DAPH.THOUS.M3")
daph
bosm <- plot.zoop.anomaly("BOSM.THOUS.M3")
ceriod <- plot.zoop.anomaly("CERIOD.THOUS.M3")
diaphan <- plot.zoop.anomaly("DIAPHAN.THOUS.M3")
lepto <- plot.zoop.anomaly("LEPTO.THOUS.M3")
cowplot::plot_grid(cyclo, cala, naup, rotif, daph,
bosm, ceriod, diaphan,
labels = "AUTO")
#ggsave("./figures/all_zoop.jpeg", height = 10, width = 15)
Can I do it in one model
long <- temp %>%
pivot_longer(cols = CYCLO.THOUS.M3:DIAPHAN.THOUS.M3,
values_to = "values",
names_to = "taxa") %>%
mutate(taxa = as.factor(taxa)
)
zoop.gam <- gam(values ~
te(JULIANDAY, ice_off_anomaly, m = 2, bs = "tp") +
t2(JULIANDAY, ice_off_anomaly, taxa,
m = 2, bs = c("tp", "tp", "re")) +
s(taxa, DNRID, bs = "re") + s(taxa, year_f, bs = "re"),
method = "fREML", select = T,
data = long, family = tw(link = "log"))
summary(zoop.gam)
appraise(zoop.gam)
draw(zoop.gam, select = 2)
#write_rds(zoop.gam, "./models/multispecies.zoop.model.rds")
#zoop.gam <- readRDS("./models/multispecies.zoop.model.rds")
Alternative plotting
groups <- unique(long$taxa)
n = length(groups)
pred.df <- data.frame(ice_off_anomaly =
c(rep(early.anomaly, 186*n), rep(0, 186*n), rep(late.anomaly, 186*n)),
JULIANDAY = seq(90, 275, by = 1),
taxa = rep(groups, each = 186),
DNRID = "62001300", year_f = "2010")
fit <- predict.gam(zoop.gam, pred.df,
unconditional = T, # Gives simultaneous
type = "response",
exclude = c('s(DNRID)',
's(year_f)'),
se.fit = T)
pred.df$predictions = fit$fit
pred.df$se = fit$se.fit
pred.df %>%
mutate(ice_off_anomaly = as.factor(ice_off_anomaly),
ice_off_descriptor =
ifelse(ice_off_anomaly == early.anomaly, "Early", "Normal"),
ice_off_descriptor =
ifelse(ice_off_anomaly == late.anomaly, "Late", ice_off_descriptor)) %>%
ggplot(aes(y = predictions, x = as.Date(JULIANDAY, origin = as.Date("2018-01-01")),
color = ice_off_anomaly)) +
#geom_point() +
geom_line(aes(y = predictions, x = as.Date(JULIANDAY, origin = as.Date("2018-01-01")),
color = ice_off_descriptor), lwd = 2) +
#geom_ribbon(aes(ymin = predictions - se, ymax = predictions + se,
# fill = ice_off_descriptor, color = ice_off_descriptor),
# lwd = 0.5, alpha = 0.05) +
facet_wrap(vars(taxa), scales = "free") +
theme_classic() +
scale_x_date(date_breaks = "months" , date_labels = "%b") +
labs(y = paste0("Zooplankton density"), x = "Julian Day",
color = "Ice-Off", fill = "Ice-Off") +
theme(axis.text = element_text(size = 14),
axis.title = element_text(size = 18),
legend.title = element_text(size = 16),
legend.text = element_text(size = 16))
ggsave("./multispecies.jpeg", width = 10, height = 8)
