Fit models to size-at-age data with temperature as covariate
Max Lindmark
11 October, 2022
Read data
load(file = "data/druksiai2.RData")
druks_df4 <- druks_df4 %>%
mutate(lifetime_avg_temp_ct = lifetime_avg_temp - mean(lifetime_avg_temp),
birth_temp_ct = birth_temp - mean(birth_temp))
#> mutate: new variable 'lifetime_avg_temp_ct' (double) with 375 unique values and 0% NA
#> new variable 'birth_temp_ct' (double) with 43 unique values and 0% NA
ggplot(druks_df4, aes(birth_temp, lifetime_avg_temp_ct)) +
geom_point() +
facet_wrap(~age)
ggplot(druks_df4, aes(lifetime_avg_temp)) +
geom_histogram() +
facet_wrap(~species, scales = "free")
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggplot(druks_df4, aes(birth_temp)) +
geom_histogram() +
facet_wrap(~species, scales = "free")
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggplot(druks_df4, aes(lifetime_avg_temp)) +
geom_histogram()
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# Filter species
bream <- druks_df4 %>% filter(species == "bream")
#> filter: removed 6,392 rows (76%), 1,999 rows remaining
perch <- druks_df4 %>% filter(species == "perch")
#> filter: removed 7,212 rows (86%), 1,179 rows remaining
cisco <- druks_df4 %>% filter(species == "cisco")
#> filter: removed 7,592 rows (90%), 799 rows remaining
pike <- druks_df4 %>% filter(species == "pike")
#> filter: removed 7,851 rows (94%), 540 rows remaining
roach <- druks_df4 %>% filter(species == "roach")
#> filter: removed 4,517 rows (54%), 3,874 rows remainingAnalysis
In this script, we fit Bayesian mixed models of length as a function of age, age squared, lifetime average temperature and the interaction between lifetime average temperature and age. The model can be written as follows:
\[ \begin{aligned} \operatorname{l}_{i} &\sim N \left(\mu, \sigma \right) \\ \mu &=\alpha_{j[i],k[i],l[i]} + \beta_{1}(\operatorname{age\_ct}) + \beta_{2}(\operatorname{lifetime\_avg\_temp}) + \beta_{3}(\operatorname{age\_sq})\ + \\ &\quad \beta_{4}(\operatorname{age\_ct} \times \operatorname{lifetime\_avg\_temp}) \\ \alpha_{j} &\sim N \left(\mu_{\alpha_{j}}, \sigma_{\alpha_{j}} \right) \text{, for year j = 1,} \dots \text{,J} \\ \alpha_{k} &\sim N \left(\mu_{\alpha_{k}}, \sigma_{\alpha_{k}} \right) \text{, for month k = 1,} \dots \text{,K} \\ \alpha_{l} &\sim N \left(\mu_{\alpha_{l}}, \sigma_{\alpha_{l}} \right) \text{, for cohort l = 1,} \dots \text{,L} \\ \end{aligned} \]
Since we fit it in a Bayesian framework using brms, we
also need to put priors on parameters. For main coefficients, \(\beta_{1:4}\), we use flat priors, and for
\(\sigma_{\alpha_{j, k, l}}\) and for
\(\sigma\), we use a \(student\)-\(t\) prior: \(student\)-\(t(3,
0, 5.2)\). These priors are uninformative defaults from
brms for our data. After we have fitted the models, we plot
parameter estimates and check model fit and chain convergence.
The figure below shows the parameter combinations we can get, and I’ve highlighted the combination of parameters that gives TSR-like temperature effects.
pal <- rev(brewer.pal(n = 6, name = "Paired")[c(6, 2)])
# This is just to get ballpark estimates, we round and change them anyway
summary(lm(l ~ age_ct*lifetime_avg_temp_ct + age_sq, data = perch))
#>
#> Call:
#> lm(formula = l ~ age_ct * lifetime_avg_temp_ct + age_sq, data = perch)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -10.0770 -2.2516 -0.5126 1.4240 17.2359
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 17.949701 0.151810 118.238 <2e-16 ***
#> age_ct 2.279536 0.043815 52.026 <2e-16 ***
#> lifetime_avg_temp_ct -1.752813 0.095107 -18.430 <2e-16 ***
#> age_sq -0.020622 0.009681 -2.130 0.0334 *
#> age_ct:lifetime_avg_temp_ct 0.065220 0.032843 1.986 0.0473 *
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 3.898 on 1174 degrees of freedom
#> Multiple R-squared: 0.77, Adjusted R-squared: 0.7692
#> F-statistic: 982.5 on 4 and 1174 DF, p-value: < 2.2e-16
plot_df <- data.frame(expand.grid(age_ct = seq(from = min(perch$age_ct),
to = max(perch$age_ct),
length.out = length(unique(perch$age_ct))),
lifetime_avg_temp_ct = c(min(round(perch$lifetime_avg_temp_ct)),
max(round(perch$lifetime_avg_temp_ct))))) %>%
mutate(age_sq = age_ct*age_ct,
age = age_ct + mean(perch$age),
temp = ifelse(lifetime_avg_temp_ct == min(lifetime_avg_temp_ct), "Cold", "Warm"), # for plotting
l_nn = 18 + 2*age_ct - 3*lifetime_avg_temp_ct - 0.01*age_sq - 0.4*age_ct*lifetime_avg_temp_ct,
l_pn = 18 + 2*age_ct + 3*lifetime_avg_temp_ct - 0.01*age_sq - 0.4*age_ct*lifetime_avg_temp_ct,
l_pp = 18 + 2*age_ct + 3*lifetime_avg_temp_ct - 0.01*age_sq + 0.4*age_ct*lifetime_avg_temp_ct,
l_np = 18 + 2*age_ct - 3*lifetime_avg_temp_ct - 0.01*age_sq + 0.4*age_ct*lifetime_avg_temp_ct) %>%
pivot_longer(6:9) %>%
mutate(temp_effect = ifelse(name %in% c("l_pn", "l_pp"),
"italic(\u03B2[temp] > 0)",
"italic(\u03B2[temp] < 0)"),
interaction_effect = ifelse(name %in% c("l_pn", "l_nn"),
"italic(\u03B2[temp%*%age] < 0)",
"italic(\u03B2[temp%*%age] > 0)"),
alpha_group = ifelse(name %in% c("l_pn"), "N", "Y"),
name = factor(name))
#> mutate: new variable 'age_sq' (double) with 17 unique values and 0% NA
#> new variable 'age' (double) with 17 unique values and 0% NA
#> new variable 'temp' (character) with 2 unique values and 0% NA
#> new variable 'l_nn' (double) with 34 unique values and 0% NA
#> new variable 'l_pn' (double) with 34 unique values and 0% NA
#> new variable 'l_pp' (double) with 34 unique values and 0% NA
#> new variable 'l_np' (double) with 34 unique values and 0% NA
#> pivot_longer: reorganized (l_nn, l_pn, l_pp, l_np) into (name, value) [was 34x9, now 136x7]
#> mutate: converted 'name' from character to factor (0 new NA)
#> new variable 'temp_effect' (character) with 2 unique values and 0% NA
#> new variable 'interaction_effect' (character) with 2 unique values and 0% NA
#> new variable 'alpha_group' (character) with 2 unique values and 0% NA
ggplot(plot_df, aes(age, value, color = temp, alpha = alpha_group)) +
geom_line(size = 1.5) +
scale_color_manual(values = pal) +
scale_alpha_manual(values = c(1, 0.3)) +
facet_grid(temp_effect ~ interaction_effect, labeller = label_parsed) +
guides(alpha = "none") +
labs(x = "Age [yrs]", y = "Length [cm]", color = "Temperature") +
theme(legend.position = "bottom") +
NULL
Fig. 1: Possible combinations of temperature (\(beta_3\)) and interaction (\(beta_4\)) coefficients (rows and columns, respectively). Combination leading to TSR-like growth is highlighted (bottom left plot).
ggsave("figures/Fig2_example_TSR_coefs.png", width = 6.5, height = 6.5, dpi = 600)In other words, \(\beta_2\) needs to be positive, and \(\beta_4\) needs to be negative for a growth according to TSR to occur.
Fit temperature models
#For now, fit toy model to extract default priors...
m0 <- brm(l ~ age_ct*lifetime_avg_temp_ct + age_sq + (1|year) + (1|cohort) + (1|month),
data = roach,
chains = 1,
iter = 1,
family = gaussian())
#> Compiling Stan program...
#> Start sampling
#>
#> SAMPLING FOR MODEL '5dfecd883da02436a96d6efd680552c9' NOW (CHAIN 1).
#> Chain 1:
#> Chain 1: Gradient evaluation took 0.000405 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 4.05 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1:
#> Chain 1:
#> Chain 1: WARNING: No variance estimation is
#> Chain 1: performed for num_warmup < 20
#> Chain 1:
#> Chain 1: Iteration: 1 / 1 [100%] (Sampling)
#> Chain 1:
#> Chain 1: Elapsed Time: 2e-06 seconds (Warm-up)
#> Chain 1: 0.000562 seconds (Sampling)
#> Chain 1: 0.000564 seconds (Total)
#> Chain 1:
#> Warning: There were 1 divergent transitions after warmup. See
#> http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
#> to find out why this is a problem and how to eliminate them.
#> Warning: Examine the pairs() plot to diagnose sampling problems
prior <- prior_summary(m0)m0b <- brm(l ~ age_ct*birth_temp_ct + age_sq + (1|year) + (1|cohort) + (1|month),
data = roach,
chains = 1,
iter = 1,
family = gaussian())
#> Compiling Stan program...
#> recompiling to avoid crashing R session
#> Start sampling
#>
#> SAMPLING FOR MODEL '5dfecd883da02436a96d6efd680552c9' NOW (CHAIN 1).
#> Chain 1:
#> Chain 1: Gradient evaluation took 0.00029 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 2.9 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1:
#> Chain 1:
#> Chain 1: WARNING: No variance estimation is
#> Chain 1: performed for num_warmup < 20
#> Chain 1:
#> Chain 1: Iteration: 1 / 1 [100%] (Sampling)
#> Chain 1:
#> Chain 1: Elapsed Time: 2e-06 seconds (Warm-up)
#> Chain 1: 0.000744 seconds (Sampling)
#> Chain 1: 0.000746 seconds (Total)
#> Chain 1:
#> Warning: There were 1 divergent transitions after warmup. See
#> http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
#> to find out why this is a problem and how to eliminate them.
#> Warning: Examine the pairs() plot to diagnose sampling problems
priorb <- prior_summary(m0b)roach_temp <- brm(l ~ age_ct*lifetime_avg_temp_ct + age_sq + (1|year) + (1|cohort) + (1|month),
data = roach,
chains = 4,
iter = 4000,
prior = prior,
family = gaussian())
print(roach_temp)
#> Family: gaussian
#> Links: mu = identity; sigma = identity
#> Formula: l ~ age_ct * lifetime_avg_temp_ct + age_sq + (1 | year) + (1 | cohort) + (1 | month)
#> Data: roach (Number of observations: 3874)
#> Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
#> total post-warmup draws = 8000
#>
#> Group-Level Effects:
#> ~cohort (Number of levels: 52)
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sd(Intercept) 0.67 0.10 0.49 0.89 1.00 1745 3065
#>
#> ~month (Number of levels: 10)
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sd(Intercept) 1.03 0.32 0.60 1.82 1.00 2378 3715
#>
#> ~year (Number of levels: 25)
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sd(Intercept) 2.56 0.39 1.92 3.46 1.00 1541 2810
#>
#> Population-Level Effects:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
#> Intercept 17.38 0.63 16.14 18.64 1.00 1055
#> age_ct 1.85 0.02 1.82 1.89 1.00 3207
#> lifetime_avg_temp_ct 1.22 0.15 0.92 1.52 1.00 2383
#> age_sq -0.03 0.00 -0.03 -0.02 1.00 6455
#> age_ct:lifetime_avg_temp_ct 0.01 0.02 -0.03 0.05 1.00 2852
#> Tail_ESS
#> Intercept 1805
#> age_ct 4692
#> lifetime_avg_temp_ct 3789
#> age_sq 6065
#> age_ct:lifetime_avg_temp_ct 4609
#>
#> Family Specific Parameters:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sigma 1.28 0.01 1.25 1.31 1.00 9249 5959
#>
#> Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
#> and Tail_ESS are effective sample size measures, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).roach_tempb <- brm(l ~ age_ct*birth_temp_ct + age_sq + (1|year) + (1|cohort) + (1|month),
data = roach,
chains = 4,
iter = 4000,
prior = priorb,
family = gaussian())
print(roach_tempb)
#> Family: gaussian
#> Links: mu = identity; sigma = identity
#> Formula: l ~ age_ct * birth_temp_ct + age_sq + (1 | year) + (1 | cohort) + (1 | month)
#> Data: roach (Number of observations: 3874)
#> Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
#> total post-warmup draws = 8000
#>
#> Group-Level Effects:
#> ~cohort (Number of levels: 52)
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sd(Intercept) 0.93 0.13 0.71 1.21 1.00 2150 4271
#>
#> ~month (Number of levels: 10)
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sd(Intercept) 0.94 0.29 0.54 1.63 1.00 3969 5292
#>
#> ~year (Number of levels: 25)
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sd(Intercept) 2.29 0.36 1.71 3.14 1.00 2421 3807
#>
#> Population-Level Effects:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
#> Intercept 17.38 0.57 16.24 18.49 1.00 2334
#> age_ct 1.85 0.02 1.81 1.89 1.00 4910
#> birth_temp_ct 0.39 0.10 0.19 0.60 1.00 3595
#> age_sq -0.03 0.00 -0.03 -0.02 1.00 7715
#> age_ct:birth_temp_ct -0.03 0.01 -0.04 -0.01 1.00 8803
#> Tail_ESS
#> Intercept 3295
#> age_ct 5656
#> birth_temp_ct 4174
#> age_sq 6697
#> age_ct:birth_temp_ct 6411
#>
#> Family Specific Parameters:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sigma 1.28 0.02 1.25 1.31 1.00 10943 5226
#>
#> Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
#> and Tail_ESS are effective sample size measures, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).pike_temp <- brm(l ~ age_ct*lifetime_avg_temp_ct + age_sq + (1|year) + (1|cohort) + (1|month),
data = pike,
chains = 4,
control = list(adapt_delta = 0.99, max_treedepth = 11), # need bc diveregent transitions
iter = 4000,
prior = prior,
family = gaussian())
print(pike_temp)
#> Family: gaussian
#> Links: mu = identity; sigma = identity
#> Formula: l ~ age_ct * lifetime_avg_temp_ct + age_sq + (1 | year) + (1 | cohort) + (1 | month)
#> Data: pike (Number of observations: 540)
#> Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
#> total post-warmup draws = 8000
#>
#> Group-Level Effects:
#> ~cohort (Number of levels: 37)
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sd(Intercept) 0.36 0.28 0.01 1.04 1.00 3808 3816
#>
#> ~month (Number of levels: 8)
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sd(Intercept) 0.42 0.37 0.01 1.38 1.00 4455 4410
#>
#> ~year (Number of levels: 20)
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sd(Intercept) 3.34 0.68 2.26 4.96 1.00 2867 4367
#>
#> Population-Level Effects:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
#> Intercept 47.65 0.94 45.80 49.48 1.00 3714
#> age_ct 7.31 0.19 6.94 7.68 1.00 9133
#> lifetime_avg_temp_ct -0.10 0.58 -1.22 1.06 1.00 4092
#> age_sq 0.31 0.07 0.17 0.44 1.00 11502
#> age_ct:lifetime_avg_temp_ct -0.03 0.12 -0.28 0.21 1.00 10575
#> Tail_ESS
#> Intercept 5143
#> age_ct 6312
#> lifetime_avg_temp_ct 4782
#> age_sq 6648
#> age_ct:lifetime_avg_temp_ct 6209
#>
#> Family Specific Parameters:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sigma 4.62 0.15 4.34 4.93 1.00 14667 5711
#>
#> Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
#> and Tail_ESS are effective sample size measures, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).pike_tempb <- brm(l ~ age_ct*birth_temp_ct + age_sq + (1|year) + (1|cohort) + (1|month),
data = pike,
chains = 4,
control = list(adapt_delta = 0.99, max_treedepth = 11), # need bc diveregent transitions
iter = 4000,
prior = priorb,
family = gaussian())
print(pike_tempb)
#> Family: gaussian
#> Links: mu = identity; sigma = identity
#> Formula: l ~ age_ct * birth_temp_ct + age_sq + (1 | year) + (1 | cohort) + (1 | month)
#> Data: pike (Number of observations: 540)
#> Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
#> total post-warmup draws = 8000
#>
#> Group-Level Effects:
#> ~cohort (Number of levels: 37)
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sd(Intercept) 0.35 0.28 0.01 1.03 1.00 3591 4000
#>
#> ~month (Number of levels: 8)
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sd(Intercept) 0.42 0.38 0.02 1.40 1.00 3910 4073
#>
#> ~year (Number of levels: 20)
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sd(Intercept) 3.29 0.63 2.28 4.72 1.00 3122 4775
#>
#> Population-Level Effects:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
#> Intercept 47.68 0.91 45.88 49.42 1.00 2527
#> age_ct 7.30 0.17 6.97 7.63 1.00 9979
#> birth_temp_ct 0.14 0.38 -0.60 0.88 1.00 5200
#> age_sq 0.31 0.07 0.17 0.44 1.00 10123
#> age_ct:birth_temp_ct -0.02 0.12 -0.25 0.23 1.00 10339
#> Tail_ESS
#> Intercept 4156
#> age_ct 6540
#> birth_temp_ct 5676
#> age_sq 5850
#> age_ct:birth_temp_ct 6164
#>
#> Family Specific Parameters:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sigma 4.63 0.14 4.35 4.92 1.00 12109 6133
#>
#> Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
#> and Tail_ESS are effective sample size measures, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).perch_temp <- brm(l ~ age_ct*lifetime_avg_temp_ct + age_sq + (1|year) + (1|cohort) + (1|month),
data = perch,
chains = 4,
iter = 4000,
prior = prior,
family = gaussian())
print(perch_temp)
#> Family: gaussian
#> Links: mu = identity; sigma = identity
#> Formula: l ~ age_ct * lifetime_avg_temp_ct + age_sq + (1 | year) + (1 | cohort) + (1 | month)
#> Data: perch (Number of observations: 1179)
#> Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
#> total post-warmup draws = 8000
#>
#> Group-Level Effects:
#> ~cohort (Number of levels: 49)
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sd(Intercept) 1.38 0.20 1.04 1.83 1.00 2339 3821
#>
#> ~month (Number of levels: 10)
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sd(Intercept) 2.21 0.70 1.28 3.99 1.00 3765 4980
#>
#> ~year (Number of levels: 19)
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sd(Intercept) 3.75 0.72 2.63 5.42 1.00 2904 4619
#>
#> Population-Level Effects:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
#> Intercept 18.02 1.15 15.74 20.25 1.00 3156
#> age_ct 2.25 0.05 2.15 2.36 1.00 4646
#> lifetime_avg_temp_ct -0.81 0.44 -1.66 0.06 1.00 3075
#> age_sq 0.00 0.01 -0.01 0.01 1.00 8602
#> age_ct:lifetime_avg_temp_ct -0.02 0.04 -0.10 0.07 1.00 4249
#> Tail_ESS
#> Intercept 3916
#> age_ct 5287
#> lifetime_avg_temp_ct 4115
#> age_sq 6666
#> age_ct:lifetime_avg_temp_ct 5267
#>
#> Family Specific Parameters:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sigma 1.70 0.04 1.63 1.77 1.00 9335 5274
#>
#> Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
#> and Tail_ESS are effective sample size measures, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).perch_tempb <- brm(l ~ age_ct*birth_temp_ct + age_sq + (1|year) + (1|cohort) + (1|month),
data = perch,
chains = 4,
iter = 4000,
prior = priorb,
family = gaussian())
print(perch_tempb)
#> Family: gaussian
#> Links: mu = identity; sigma = identity
#> Formula: l ~ age_ct * birth_temp_ct + age_sq + (1 | year) + (1 | cohort) + (1 | month)
#> Data: perch (Number of observations: 1179)
#> Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
#> total post-warmup draws = 8000
#>
#> Group-Level Effects:
#> ~cohort (Number of levels: 49)
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sd(Intercept) 1.31 0.19 0.99 1.73 1.00 2418 3983
#>
#> ~month (Number of levels: 10)
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sd(Intercept) 2.24 0.73 1.28 4.07 1.00 4047 4939
#>
#> ~year (Number of levels: 19)
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sd(Intercept) 3.99 0.75 2.79 5.68 1.00 3124 3852
#>
#> Population-Level Effects:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
#> Intercept 18.19 1.24 15.75 20.61 1.00 2394
#> age_ct 2.23 0.05 2.14 2.33 1.00 5421
#> birth_temp_ct -0.30 0.19 -0.68 0.06 1.00 4514
#> age_sq -0.00 0.01 -0.02 0.01 1.00 8316
#> age_ct:birth_temp_ct -0.00 0.03 -0.06 0.06 1.00 6019
#> Tail_ESS
#> Intercept 3379
#> age_ct 5172
#> birth_temp_ct 4903
#> age_sq 5615
#> age_ct:birth_temp_ct 5137
#>
#> Family Specific Parameters:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sigma 1.70 0.04 1.63 1.77 1.00 9227 5383
#>
#> Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
#> and Tail_ESS are effective sample size measures, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).bream_temp <- brm(l ~ age_ct*lifetime_avg_temp_ct + age_sq + (1|year) + (1|cohort) + (1|month),
data = bream,
chains = 4,
control = list(adapt_delta = 0.95), # need bc diveregent transitions
iter = 4000,
prior = prior,
family = gaussian())
print(bream_temp)
#> Family: gaussian
#> Links: mu = identity; sigma = identity
#> Formula: l ~ age_ct * lifetime_avg_temp_ct + age_sq + (1 | year) + (1 | cohort) + (1 | month)
#> Data: bream (Number of observations: 1999)
#> Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
#> total post-warmup draws = 8000
#>
#> Group-Level Effects:
#> ~cohort (Number of levels: 51)
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sd(Intercept) 0.94 0.14 0.70 1.25 1.00 2264 3582
#>
#> ~month (Number of levels: 9)
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sd(Intercept) 0.78 0.27 0.42 1.46 1.00 2530 4125
#>
#> ~year (Number of levels: 23)
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sd(Intercept) 1.95 0.32 1.43 2.68 1.00 2190 3874
#>
#> Population-Level Effects:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
#> Intercept 26.97 0.53 25.92 28.03 1.00 1261
#> age_ct 3.31 0.03 3.25 3.37 1.00 2851
#> lifetime_avg_temp_ct -0.07 0.23 -0.51 0.38 1.00 1713
#> age_sq -0.09 0.01 -0.10 -0.08 1.00 6332
#> age_ct:lifetime_avg_temp_ct 0.14 0.03 0.07 0.21 1.00 2056
#> Tail_ESS
#> Intercept 2428
#> age_ct 4123
#> lifetime_avg_temp_ct 3272
#> age_sq 6435
#> age_ct:lifetime_avg_temp_ct 4012
#>
#> Family Specific Parameters:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sigma 2.04 0.03 1.97 2.10 1.00 10663 5836
#>
#> Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
#> and Tail_ESS are effective sample size measures, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).bream_tempb <- brm(l ~ age_ct*birth_temp_ct + age_sq + (1|year) + (1|cohort) + (1|month),
data = bream,
chains = 4,
control = list(adapt_delta = 0.95), # need bc diveregent transitions
iter = 4000,
prior = priorb,
family = gaussian())
print(bream_tempb)
#> Family: gaussian
#> Links: mu = identity; sigma = identity
#> Formula: l ~ age_ct * birth_temp_ct + age_sq + (1 | year) + (1 | cohort) + (1 | month)
#> Data: bream (Number of observations: 1999)
#> Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
#> total post-warmup draws = 8000
#>
#> Group-Level Effects:
#> ~cohort (Number of levels: 51)
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sd(Intercept) 1.13 0.15 0.86 1.46 1.00 2188 3703
#>
#> ~month (Number of levels: 9)
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sd(Intercept) 0.77 0.27 0.41 1.43 1.00 2937 4448
#>
#> ~year (Number of levels: 23)
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sd(Intercept) 2.08 0.34 1.54 2.88 1.00 1929 3965
#>
#> Population-Level Effects:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
#> Intercept 26.94 0.57 25.81 28.04 1.00 1460
#> age_ct 3.28 0.03 3.22 3.35 1.00 2746
#> birth_temp_ct -0.33 0.13 -0.59 -0.07 1.00 2220
#> age_sq -0.09 0.01 -0.10 -0.08 1.00 5692
#> age_ct:birth_temp_ct 0.02 0.02 -0.01 0.05 1.00 4510
#> Tail_ESS
#> Intercept 2736
#> age_ct 4501
#> birth_temp_ct 3594
#> age_sq 5826
#> age_ct:birth_temp_ct 5702
#>
#> Family Specific Parameters:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sigma 2.04 0.03 1.97 2.11 1.00 10701 5460
#>
#> Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
#> and Tail_ESS are effective sample size measures, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).cisco_temp <- brm(l ~ age_ct*lifetime_avg_temp_ct + age_sq + (1|year) + (1|cohort) + (1|month),
data = cisco,
chains = 4,
control = list(adapt_delta = 0.99, max_treedepth = 11), # need bc diveregent transitions
iter = 4000,
prior = prior,
family = gaussian())
print(cisco_temp)
#> Family: gaussian
#> Links: mu = identity; sigma = identity
#> Formula: l ~ age_ct * lifetime_avg_temp_ct + age_sq + (1 | year) + (1 | cohort) + (1 | month)
#> Data: cisco (Number of observations: 799)
#> Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
#> total post-warmup draws = 8000
#>
#> Group-Level Effects:
#> ~cohort (Number of levels: 28)
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sd(Intercept) 0.75 0.15 0.51 1.09 1.00 3446 5234
#>
#> ~month (Number of levels: 6)
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sd(Intercept) 1.48 0.66 0.72 3.19 1.00 5497 5177
#>
#> ~year (Number of levels: 19)
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sd(Intercept) 1.28 0.27 0.86 1.92 1.00 4516 5080
#>
#> Population-Level Effects:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
#> Intercept 16.04 0.73 14.62 17.50 1.00 4672
#> age_ct 2.58 0.10 2.40 2.78 1.00 7713
#> lifetime_avg_temp_ct 0.06 0.21 -0.35 0.48 1.00 5061
#> age_sq -0.46 0.05 -0.55 -0.37 1.00 9237
#> age_ct:lifetime_avg_temp_ct -0.15 0.06 -0.27 -0.03 1.00 7559
#> Tail_ESS
#> Intercept 5333
#> age_ct 6265
#> lifetime_avg_temp_ct 5470
#> age_sq 6594
#> age_ct:lifetime_avg_temp_ct 5947
#>
#> Family Specific Parameters:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sigma 0.74 0.02 0.71 0.78 1.00 11816 5855
#>
#> Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
#> and Tail_ESS are effective sample size measures, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).cisco_tempb <- brm(l ~ age_ct*birth_temp_ct + age_sq + (1|year) + (1|cohort) + (1|month),
data = cisco,
chains = 4,
control = list(adapt_delta = 0.99, max_treedepth = 11), # need bc diveregent transitions
iter = 4000,
prior = priorb,
family = gaussian())
print(cisco_tempb)
#> Family: gaussian
#> Links: mu = identity; sigma = identity
#> Formula: l ~ age_ct * birth_temp_ct + age_sq + (1 | year) + (1 | cohort) + (1 | month)
#> Data: cisco (Number of observations: 799)
#> Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
#> total post-warmup draws = 8000
#>
#> Group-Level Effects:
#> ~cohort (Number of levels: 28)
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sd(Intercept) 0.75 0.15 0.51 1.09 1.00 3574 5408
#>
#> ~month (Number of levels: 6)
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sd(Intercept) 1.47 0.63 0.73 3.10 1.00 5715 5766
#>
#> ~year (Number of levels: 19)
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sd(Intercept) 1.25 0.25 0.86 1.84 1.00 4181 5189
#>
#> Population-Level Effects:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
#> Intercept 15.98 0.71 14.59 17.38 1.00 5417
#> age_ct 2.60 0.10 2.40 2.80 1.00 6515
#> birth_temp_ct 0.09 0.21 -0.32 0.50 1.00 5932
#> age_sq -0.44 0.05 -0.53 -0.35 1.00 9042
#> age_ct:birth_temp_ct -0.13 0.06 -0.25 -0.01 1.00 7415
#> Tail_ESS
#> Intercept 5202
#> age_ct 5690
#> birth_temp_ct 5205
#> age_sq 6218
#> age_ct:birth_temp_ct 5727
#>
#> Family Specific Parameters:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sigma 0.74 0.02 0.71 0.78 1.00 11475 5681
#>
#> Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
#> and Tail_ESS are effective sample size measures, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).Plot coefficients
# Extract temperature coefficients
roach_post_temp_temp <- roach_temp %>%
gather_draws(b_lifetime_avg_temp_ct) %>%
mutate(species = "roach")
pike_post_temp_temp <- pike_temp %>%
gather_draws(b_lifetime_avg_temp_ct) %>%
mutate(species = "pike")
perch_post_temp_temp <- perch_temp %>%
gather_draws(b_lifetime_avg_temp_ct) %>%
mutate(species = "perch")
bream_post_temp_temp <- bream_temp %>%
gather_draws(b_lifetime_avg_temp_ct) %>%
mutate(species = "bream")
cisco_post_temp_temp <- cisco_temp %>%
gather_draws(b_lifetime_avg_temp_ct) %>%
mutate(species = "cisco")
post_temp_temp <- bind_rows(perch_post_temp_temp, roach_post_temp_temp,
bream_post_temp_temp, cisco_post_temp_temp,
pike_post_temp_temp) %>%
mutate(parameter = "italic(\u03B2[temp])")
# post_temp_temp %>%
# ggplot(aes(x = .value, y = species, fill = species, color = species)) +
# #ggplot(aes(x = .value, y = species)) +
# geom_vline(xintercept = 0, linetype = 2, color = "grey20") +
# stat_halfeye(alpha = 0.5, size = 5, .width = c(0.7)) +
# #stat_dots(quantiles = 100, color = "grey10") +
# scale_fill_brewer(palette = "Dark2") +
# scale_color_brewer(palette = "Dark2") +
# guides(fill = "none", color = "none") +
# labs(x = expression(paste(italic(beta[2]), " (temperature coefficient)")), fill = "", y = "density") +
# theme(legend.position = c(0.9, 0.9),
# legend.key.size = unit(0.2, "cm"),
# legend.background = element_blank())
# Extract temperature coefficients from BIRTH temperature model
roach_post_temp_tempb <- roach_tempb %>%
gather_draws(b_birth_temp_ct) %>%
mutate(species = "roach")
pike_post_temp_tempb <- pike_tempb %>%
gather_draws(b_birth_temp_ct) %>%
mutate(species = "pike")
perch_post_temp_tempb <- perch_tempb %>%
gather_draws(b_birth_temp_ct) %>%
mutate(species = "perch")
bream_post_temp_tempb <- bream_tempb %>%
gather_draws(b_birth_temp_ct) %>%
mutate(species = "bream")
cisco_post_temp_tempb <- cisco_tempb %>%
gather_draws(b_birth_temp_ct) %>%
mutate(species = "cisco")
post_temp_tempb <- bind_rows(perch_post_temp_tempb, roach_post_temp_tempb,
bream_post_temp_tempb, cisco_post_temp_tempb,
pike_post_temp_tempb) %>%
mutate(parameter = "italic(\u03B2[temp])")# Extract interaction coefficients from temperature model
perch_post_temp_inter <- perch_temp %>%
gather_draws(`b_age_ct:lifetime_avg_temp_ct`) %>%
mutate(species = "perch")
roach_post_temp_inter <- roach_temp %>%
gather_draws(`b_age_ct:lifetime_avg_temp_ct`) %>%
mutate(species = "roach")
bream_post_temp_inter <- bream_temp %>%
gather_draws(`b_age_ct:lifetime_avg_temp_ct`) %>%
mutate(species = "bream")
cisco_post_temp_inter <- cisco_temp %>%
gather_draws(`b_age_ct:lifetime_avg_temp_ct`) %>%
mutate(species = "cisco")
pike_post_temp_inter <- pike_temp %>%
gather_draws(`b_age_ct:lifetime_avg_temp_ct`) %>%
mutate(species = "pike")
post_temp_inter <- bind_rows(perch_post_temp_inter, roach_post_temp_inter,
bream_post_temp_inter, cisco_post_temp_inter,
pike_post_temp_inter) %>%
mutate(parameter = "italic(\u03B2[temp%*%age])")
# Extract interaction coefficients from birth temperature model
perch_post_temp_interb <- perch_tempb %>%
gather_draws(`b_age_ct:birth_temp_ct`) %>%
mutate(species = "perch")
roach_post_temp_interb <- roach_tempb %>%
gather_draws(`b_age_ct:birth_temp_ct`) %>%
mutate(species = "roach")
bream_post_temp_interb <- bream_tempb %>%
gather_draws(`b_age_ct:birth_temp_ct`) %>%
mutate(species = "bream")
cisco_post_temp_interb <- cisco_tempb %>%
gather_draws(`b_age_ct:birth_temp_ct`) %>%
mutate(species = "cisco")
pike_post_temp_interb <- pike_tempb %>%
gather_draws(`b_age_ct:birth_temp_ct`) %>%
mutate(species = "pike")
post_temp_interb <- bind_rows(perch_post_temp_interb, roach_post_temp_interb,
bream_post_temp_interb, cisco_post_temp_interb,
pike_post_temp_interb) %>%
mutate(parameter = "italic(\u03B2[temp%*%age])")post_temp_temp <- bind_rows(mutate(post_temp_temp, temperature = "avg. lifetime"),
mutate(post_temp_tempb, temperature = "birth-year"))
#> mutate (grouped): new variable 'temperature' (character) with one unique value and 0% NA
#> mutate (grouped): new variable 'temperature' (character) with one unique value and 0% NA
post_temp_inter <- bind_rows(mutate(post_temp_inter, temperature = "avg. lifetime"),
mutate(post_temp_interb, temperature = "birth-year"))
#> mutate (grouped): new variable 'temperature' (character) with one unique value and 0% NA
#> mutate (grouped): new variable 'temperature' (character) with one unique value and 0% NAtemp_coefs <- bind_rows(post_temp_temp, post_temp_inter)
ggplot(temp_coefs, aes(x = .value, y = species, color = temperature)) +
geom_rect(data = data.frame(parameter = "italic(\u03B2[temp])"),
aes(ymin = -Inf, ymax = Inf, xmin = 0, xmax = Inf, fill = "TSR-region"),
inherit.aes = FALSE, color = NA) +
geom_rect(data = data.frame(parameter = "italic(\u03B2[temp%*%age])"),
aes(ymin = -Inf, ymax = Inf, xmin = -Inf, xmax = 0, fill = "TSR-region"),
inherit.aes = FALSE, color = NA) +
stat_pointinterval(position = position_dodge(width = 0.2), .width = c(0.75, 0.95)) +
facet_wrap(~parameter, scales = "free_x", labeller = label_parsed) +
scale_color_brewer(palette = "Dark2", name = "Temperature-variable") +
scale_fill_manual(values = "white") +
labs(y = "Species", fill = "", x = "Value") +
theme(legend.position = "bottom",
panel.background = element_rect(fill = 'grey85', colour = 'grey85')) +
guides(fill = guide_legend(override.aes = list(color = "grey50"))) +
NULL
Fig. 2: Estimates of temperature coefficient (beta_temp) and interaction coefficient (beta_temp%%age*), indicated by color. Lines indicate 90% credible intervals
ggsave("figures/Fig3_temp_coeffs.png", width = 6.5, height = 6.5, dpi = 600)
ggplot(temp_coefs, aes(x = .value, y = species, color = temperature)) +
geom_rect(data = data.frame(parameter = "italic(\u03B2[temp])"),
aes(ymin = -Inf, ymax = Inf, xmin = 0, xmax = Inf, fill = "TSR-region"),
inherit.aes = FALSE, color = NA) +
geom_rect(data = data.frame(parameter = "italic(\u03B2[temp%*%age])"),
aes(ymin = -Inf, ymax = Inf, xmin = -Inf, xmax = 0, fill = "TSR-region"),
inherit.aes = FALSE, color = NA) +
stat_pointinterval(position = position_dodge(width = 0.2), .width = c(0.75, 0.95)) +
facet_wrap(~parameter, scales = "free_x", labeller = label_parsed) +
scale_color_brewer(palette = "Dark2", name = "Temperature-variable") +
scale_fill_manual(values = "white") +
labs(x = "Species", fill = "", y = "Value") +
theme(legend.position = "bottom",
panel.background = element_rect(fill = 'grey85', colour = 'grey85')) +
guides(fill = guide_legend(override.aes = list(color = "grey50"))) +
NULL
Fig. 2: Estimates of temperature coefficient (beta_temp) and interaction coefficient (beta_temp%%age*), indicated by color. Lines indicate 90% credible intervals
Plot predictions and data
# I want to predict at 16 and 19 C, because that's the range for most species
roach_pred_temp <- roach %>%
data_grid(age_ct = seq_range(age_ct, by = 0.5),
lifetime_avg_temp_ct = c(16 - mean(roach$lifetime_avg_temp), 19 - mean(roach$lifetime_avg_temp))) %>%
mutate(age_sq = age_ct*age_ct,
age = age_ct + mean(roach$age),
lifetime_avg_temp = lifetime_avg_temp_ct + mean(roach$lifetime_avg_temp)) %>%
add_predicted_draws(roach_temp, re_formula = NA) %>%
mutate(species = "roach")
#> mutate: new variable 'age_sq' (double) with 31 unique values and 0% NA
#> new variable 'age' (double) with 31 unique values and 0% NA
#> new variable 'lifetime_avg_temp' (double) with 2 unique values and 0% NA
#> mutate (grouped): new variable 'species' (character) with one unique value and 0% NA
pike_pred_temp <- pike %>%
data_grid(age_ct = seq_range(age_ct, by = 0.5),
lifetime_avg_temp_ct = c(16 - mean(pike$lifetime_avg_temp), 19 - mean(pike$lifetime_avg_temp))) %>%
mutate(age_sq = age_ct*age_ct,
age = age_ct + mean(pike$age),
lifetime_avg_temp = lifetime_avg_temp_ct + mean(pike$lifetime_avg_temp)) %>%
add_predicted_draws(pike_temp, re_formula = NA) %>%
mutate(species = "pike")
#> mutate: new variable 'age_sq' (double) with 19 unique values and 0% NA
#> new variable 'age' (double) with 19 unique values and 0% NA
#> new variable 'lifetime_avg_temp' (double) with 2 unique values and 0% NA
#> mutate (grouped): new variable 'species' (character) with one unique value and 0% NA
perch_pred_temp <- perch %>%
data_grid(age_ct = seq_range(age_ct, by = 0.5),
lifetime_avg_temp_ct = c(16 - mean(perch$lifetime_avg_temp), 19 - mean(perch$lifetime_avg_temp))) %>%
mutate(age_sq = age_ct*age_ct,
age = age_ct + mean(perch$age),
lifetime_avg_temp = lifetime_avg_temp_ct + mean(perch$lifetime_avg_temp)) %>%
add_predicted_draws(perch_temp, re_formula = NA) %>%
mutate(species = "perch")
#> mutate: new variable 'age_sq' (double) with 33 unique values and 0% NA
#> new variable 'age' (double) with 33 unique values and 0% NA
#> new variable 'lifetime_avg_temp' (double) with 2 unique values and 0% NA
#> mutate (grouped): new variable 'species' (character) with one unique value and 0% NA
bream_pred_temp <- bream %>%
data_grid(age_ct = seq_range(age_ct, by = 0.5),
lifetime_avg_temp_ct = c(16 - mean(bream$lifetime_avg_temp), 19 - mean(bream$lifetime_avg_temp))) %>%
mutate(age_sq = age_ct*age_ct,
age = age_ct + mean(bream$age),
lifetime_avg_temp = lifetime_avg_temp_ct + mean(bream$lifetime_avg_temp)) %>%
add_predicted_draws(bream_temp, re_formula = NA) %>%
mutate(species = "bream")
#> mutate: new variable 'age_sq' (double) with 33 unique values and 0% NA
#> new variable 'age' (double) with 33 unique values and 0% NA
#> new variable 'lifetime_avg_temp' (double) with 2 unique values and 0% NA
#> mutate (grouped): new variable 'species' (character) with one unique value and 0% NA
cisco_pred_temp <- cisco %>%
data_grid(age_ct = seq_range(age_ct, by = 0.5),
lifetime_avg_temp_ct = c(16 - mean(cisco$lifetime_avg_temp), 19 - mean(cisco$lifetime_avg_temp))) %>%
mutate(age_sq = age_ct*age_ct,
age = age_ct + mean(cisco$age),
lifetime_avg_temp = lifetime_avg_temp_ct + mean(cisco$lifetime_avg_temp)) %>%
add_predicted_draws(cisco_temp, re_formula = NA) %>%
mutate(species = "cisco")
#> mutate: new variable 'age_sq' (double) with 11 unique values and 0% NA
#> new variable 'age' (double) with 11 unique values and 0% NA
#> new variable 'lifetime_avg_temp' (double) with 2 unique values and 0% NA
#> mutate (grouped): new variable 'species' (character) with one unique value and 0% NA
all_pred_temp <- bind_rows(bream_pred_temp, cisco_pred_temp, perch_pred_temp,
pike_pred_temp, roach_pred_temp)
# Sample size facets
f_labels <- data.frame(species = c("bream", "cisco", "perch", "pike", "roach"),
label = c(paste("n =", nrow(bream)),
paste("n =", nrow(cisco)),
paste("n =", nrow(perch)),
paste("n =", nrow(pike)),
paste("n =", nrow(roach))))
stage_pal <- RColorBrewer::brewer.pal(n = 9, name = "Set1")[c(9, 2, 1)]
plot(1:2, col = stage_pal[2], type = "l", lwd = 50)
Fig. 2: Predicted length at age for two temperatures (16°C in blue and 19°C in red). Bands indicate 90% credible intervals.
plot(1:2, col = stage_pal[3], type = "l", lwd = 50)
Fig. 2: Predicted length at age for two temperatures (16°C in blue and 19°C in red). Bands indicate 90% credible intervals.
ggplot(all_pred_temp, aes(x = age, y = l, color = factor(lifetime_avg_temp),
fill = factor(lifetime_avg_temp))) +
geom_jitter(data = filter(druks_df4, lifetime_avg_temp_ct < 0), aes(age, l),
inherit.aes = FALSE, width = 0.2, height = 0,
alpha = 0.6, size = 1.2, shape = 21, fill = stage_pal[2], color = "white") +
geom_jitter(data = filter(druks_df4, lifetime_avg_temp_ct > 0), aes(age, l),
inherit.aes = FALSE, width = 0.2, height = 0,
alpha = 0.6, size = 1.2, shape = 21, fill = stage_pal[3], color = "white") +
facet_wrap(~species, scales = "free", ncol = 2) +
stat_lineribbon(aes(y = .prediction), .width = c(0.9), alpha = 0.2, size = 0.8) +
stat_lineribbon(aes(y = .prediction), .width = c(0), alpha = 0.8, size = 0.8) +
scale_fill_manual(values = stage_pal[2:3]) +
scale_color_manual(values = stage_pal[2:3]) +
guides(color = "none") +
labs(y = "Length [cm]", x = "Age [yrs]", fill = "Temperature [°C]") +
geom_text(data = f_labels, aes(label = label), inherit.aes = FALSE,
x = -Inf, y = Inf, size = 3, hjust = -0.1, vjust = 1.1) +
theme(legend.position = c(1, 0),
legend.justification = c(1.4, -0.2)) +
NULL
#> filter: removed 4,192 rows (50%), 4,199 rows remaining
#> filter: removed 4,199 rows (50%), 4,192 rows remaining
Fig. 2: Predicted length at age for two temperatures (16°C in blue and 19°C in red). Bands indicate 90% credible intervals.
ggsave("figures/Fig4_pred_data_temp.png", width = 6.5, height = 6.5, dpi = 600)
ggplot(all_pred_temp, aes(x = age, y = l, color = factor(lifetime_avg_temp),
fill = factor(lifetime_avg_temp))) +
facet_wrap(~species, scales = "free", ncol = 2) +
stat_lineribbon(aes(y = .prediction), .width = c(0), alpha = 0.8, size = 0.8) +
NULL
Fig. 2: Predicted length at age for two temperatures (16°C in blue and 19°C in red). Bands indicate 90% credible intervals.
Same for birth-year temperature
# I want to predict at 16 and 19 C, because that's the range for most species
roach_pred_tempb <- roach %>%
data_grid(age_ct = seq_range(age_ct, by = 0.5),
birth_temp_ct = c(16 - mean(roach$birth_temp), 19 - mean(roach$birth_temp))) %>%
mutate(age_sq = age_ct*age_ct,
age = age_ct + mean(roach$age),
birth_temp = birth_temp_ct + mean(roach$birth_temp)) %>%
add_predicted_draws(roach_tempb, re_formula = NA) %>%
mutate(species = "roach")
#> mutate: new variable 'age_sq' (double) with 31 unique values and 0% NA
#> new variable 'age' (double) with 31 unique values and 0% NA
#> new variable 'birth_temp' (double) with 2 unique values and 0% NA
#> mutate (grouped): new variable 'species' (character) with one unique value and 0% NA
pike_pred_tempb <- pike %>%
data_grid(age_ct = seq_range(age_ct, by = 0.5),
birth_temp_ct = c(16 - mean(pike$birth_temp), 19 - mean(pike$birth_temp))) %>%
mutate(age_sq = age_ct*age_ct,
age = age_ct + mean(pike$age),
birth_temp = birth_temp_ct + mean(pike$birth_temp)) %>%
add_predicted_draws(pike_tempb, re_formula = NA) %>%
mutate(species = "pike")
#> mutate: new variable 'age_sq' (double) with 19 unique values and 0% NA
#> new variable 'age' (double) with 19 unique values and 0% NA
#> new variable 'birth_temp' (double) with 2 unique values and 0% NA
#> mutate (grouped): new variable 'species' (character) with one unique value and 0% NA
perch_pred_tempb <- perch %>%
data_grid(age_ct = seq_range(age_ct, by = 0.5),
birth_temp_ct = c(16 - mean(perch$birth_temp), 19 - mean(perch$birth_temp))) %>%
mutate(age_sq = age_ct*age_ct,
age = age_ct + mean(perch$age),
birth_temp = birth_temp_ct + mean(perch$birth_temp)) %>%
add_predicted_draws(perch_tempb, re_formula = NA) %>%
mutate(species = "perch")
#> mutate: new variable 'age_sq' (double) with 33 unique values and 0% NA
#> new variable 'age' (double) with 33 unique values and 0% NA
#> new variable 'birth_temp' (double) with 2 unique values and 0% NA
#> mutate (grouped): new variable 'species' (character) with one unique value and 0% NA
bream_pred_tempb <- bream %>%
data_grid(age_ct = seq_range(age_ct, by = 0.5),
birth_temp_ct = c(16 - mean(bream$birth_temp), 19 - mean(bream$birth_temp))) %>%
mutate(age_sq = age_ct*age_ct,
age = age_ct + mean(bream$age),
birth_temp = birth_temp_ct + mean(bream$birth_temp)) %>%
add_predicted_draws(bream_tempb, re_formula = NA) %>%
mutate(species = "bream")
#> mutate: new variable 'age_sq' (double) with 33 unique values and 0% NA
#> new variable 'age' (double) with 33 unique values and 0% NA
#> new variable 'birth_temp' (double) with 2 unique values and 0% NA
#> mutate (grouped): new variable 'species' (character) with one unique value and 0% NA
cisco_pred_tempb <- cisco %>%
data_grid(age_ct = seq_range(age_ct, by = 0.5),
birth_temp_ct = c(16 - mean(cisco$birth_temp), 19 - mean(cisco$birth_temp))) %>%
mutate(age_sq = age_ct*age_ct,
age = age_ct + mean(cisco$age),
birth_temp = birth_temp_ct + mean(cisco$birth_temp)) %>%
add_predicted_draws(cisco_tempb, re_formula = NA) %>%
mutate(species = "cisco")
#> mutate: new variable 'age_sq' (double) with 11 unique values and 0% NA
#> new variable 'age' (double) with 11 unique values and 0% NA
#> new variable 'birth_temp' (double) with 2 unique values and 0% NA
#> mutate (grouped): new variable 'species' (character) with one unique value and 0% NA
all_pred_tempb <- bind_rows(bream_pred_tempb, cisco_pred_tempb, perch_pred_tempb,
pike_pred_tempb, roach_pred_tempb)
# Sample size facets
f_labels <- data.frame(species = c("bream", "cisco", "perch", "pike", "roach"),
label = c(paste("n =", nrow(bream)),
paste("n =", nrow(cisco)),
paste("n =", nrow(perch)),
paste("n =", nrow(pike)),
paste("n =", nrow(roach))))
stage_pal <- RColorBrewer::brewer.pal(n = 9, name = "Set1")[c(9, 2, 1)]
ggplot(all_pred_tempb, aes(x = age, y = l, color = factor(birth_temp), fill = factor(birth_temp))) +
geom_jitter(data = filter(druks_df4, birth_temp_ct < 0), aes(age, l),
inherit.aes = FALSE, width = 0.2, height = 0,
alpha = 0.6, size = 1.2, shape = 21, fill = stage_pal[2], color = "white") +
geom_jitter(data = filter(druks_df4, birth_temp_ct > 0), aes(age, l),
inherit.aes = FALSE, width = 0.2, height = 0,
alpha = 0.6, size = 1.2, shape = 21, fill = stage_pal[3], color = "white") +
facet_wrap(~species, scales = "free", ncol = 2) +
stat_lineribbon(aes(y = .prediction), .width = c(0.9), alpha = 0.2, size = 0.8) +
stat_lineribbon(aes(y = .prediction), .width = c(0), alpha = 0.8, size = 0.8) +
scale_fill_manual(values = stage_pal[2:3]) +
scale_color_manual(values = stage_pal[2:3]) +
guides(color = "none") +
labs(y = "Length [cm]", x = "Age [yrs]", fill = "Temperature [°C]") +
geom_text(data = f_labels, aes(label = label), inherit.aes = FALSE,
x = -Inf, y = Inf, size = 3, hjust = -0.1, vjust = 1.1) +
theme(legend.position = c(1, 0),
legend.justification = c(1.4, -0.2)) +
NULL
#> filter: removed 4,156 rows (50%), 4,235 rows remaining
#> filter: removed 4,235 rows (50%), 4,156 rows remaining
Fig. 2: Predicted length at age for two temperatures (16°C in blue and 20°C in red). Bands indicate 90% credible intervals.
ggsave("figures/supp/Fig4b_pred_data_temp.png", width = 6.5, height = 6.5, dpi = 600)all_pred_tempb %>%
ungroup() %>%
dplyr::select(age, birth_temp, .prediction, species, .draw) %>%
as.data.frame() %>%
pivot_wider(names_from = birth_temp, values_from = .prediction) %>%
mutate(diff = `19` - `16`) %>%
group_by(species, age) %>%
summarise(median_diff = median(diff),
upr_diff = quantile(diff, probs = 0.25),
lwr_diff = quantile(diff, probs = 0.95)) %>%
ggplot(aes(x = age, y = median_diff)) +
#facet_wrap(~species, scales = "free_x", ncol = 3) +
facet_wrap(~species, scales = "free", ncol = 3) +
geom_ribbon(aes(x = age, ymin = lwr_diff, ymax = upr_diff), alpha = 0.4) +
geom_line(alpha = 0.8, size = 0.75) +
scale_y_continuous(breaks = seq(-10, 10, by = 2)) +
geom_hline(yintercept = 0, linetype = 2, color = "tomato") +
labs(y = "Difference in length (warm - cold) [cm]", x = "Age [yrs]", fill = "Temperature [°C]") +
geom_text(data = f_labels, aes(label = label), inherit.aes = FALSE,
x = -Inf, y = Inf, size = 3, hjust = -0.1, vjust = 1.1) +
theme(legend.position = c(1, 0),
legend.justification = c(1.4, -0.2)) +
NULL
#> ungroup: no grouping variables
#> pivot_wider: reorganized (birth_temp, .prediction) into (16, 19) [was 2032000x5, now 1016000x5]
#> mutate: new variable 'diff' (double) with 1,016,000 unique values and 0% NA
#> group_by: 2 grouping variables (species, age)
#> summarise: now 127 rows and 5 columns, one group variable remaining (species)
ggsave("figures/Fig5_pred_diff.png", width = 6.5, height = 6.5, dpi = 600)
all_pred_tempb %>%
ungroup() %>%
dplyr::select(age, birth_temp, .prediction, species, .draw) %>%
as.data.frame() %>%
pivot_wider(names_from = birth_temp, values_from = .prediction) %>%
mutate(diff = `19` - `16`) %>%
group_by(species, age) %>%
summarise(median_diff = median(diff),
upr_diff = quantile(diff, probs = 0.2),
lwr_diff = quantile(diff, probs = 0.9)) %>%
ggplot(aes(x = age, y = median_diff, color = species, fill = species)) +
geom_ribbon(aes(x = age, ymin = lwr_diff, ymax = upr_diff), alpha = 0.2, color = NA) +
geom_line(alpha = 1, size = 0.75) +
scale_y_continuous(breaks = seq(-10, 10, by = 2)) +
geom_hline(yintercept = 0, linetype = 2, color = "black") +
labs(y = "Difference in length (warm - cold) [cm]", x = "Age [yrs]", fill = "Species") +
guides(color = "none") +
scale_color_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2") +
NULL
#> ungroup: no grouping variables
#> pivot_wider: reorganized (birth_temp, .prediction) into (16, 19) [was 2032000x5, now 1016000x5]
#> mutate: new variable 'diff' (double) with 1,016,000 unique values and 0% NA
#> group_by: 2 grouping variables (species, age)
#> summarise: now 127 rows and 5 columns, one group variable remaining (species)
Results and discussion
We find that only cisco has parameter combinations that would result in TSR-like growth; that is, a positive temperature coefficient (size is larger when it gets warmer), and a negative temperature*age coefficient (the positive effect of temperature diminishes with age) (Fig. 1-2). This leads to cisco younger than 2 years benefit from warming (larger size-at-age), while older fish become smaller relative to colder temperatures. The results from the remaining four species are more scattered, with species showing a combinations of higher and lower intercepts and age coefficients in the periods. This result is also evident in Fig. 3, which depicts data and predicted length-at-age for all time periods for all species. For instance, for bream we predict an opposite TSR (smaller size at age for small fish and larger size at age for larger fish in warm temperatures). Perch size-at-age is larger in cold temperatures for all ages (interaction coefficient close to 0), and roach has a larger size-at-age in warm temperatures for all sizes. For pike we cannot detect a temperature-dependence in size-at-age (both temperature coefficient and temperature-age interaction are near 0).
Supplementary analysis
Check fit
# Posterior predictive checks
pp_roach <- pp_check(roach_temp) + ggtitle("roach") + guides(color = "none")
#> Using 10 posterior draws for ppc type 'dens_overlay' by default.
pp_pike <- pp_check(pike_temp) + ggtitle("pike")
#> Using 10 posterior draws for ppc type 'dens_overlay' by default.
pp_perch <- pp_check(perch_temp) + ggtitle("perch") + guides(color = "none")
#> Using 10 posterior draws for ppc type 'dens_overlay' by default.
pp_bream <- pp_check(bream_temp) + ggtitle("bream") + guides(color = "none")
#> Using 10 posterior draws for ppc type 'dens_overlay' by default.
pp_cisco <- pp_check(cisco_temp) + ggtitle("cisco") + guides(color = "none")
#> Using 10 posterior draws for ppc type 'dens_overlay' by default.
pp_perch + pp_roach + pp_bream + pp_cisco + pp_pike +
plot_layout(ncol = 2) + plot_annotation(title = "Posterior predictive checks for temperature models (avg. life temp)")
Fig. S1: Posterior predictive checks for the temperature models
ggsave("figures/supp/pp_check_temp.png", width = 6.5, height = 6.5, dpi = 600)
pp_roach_mean <- pp_check(roach_temp, type = 'stat', stat = 'mean') + ggtitle("roach") + guides(color = "none")
#> Using all posterior draws for ppc type 'stat' by default.
pp_pike_mean <- pp_check(pike_temp, type = 'stat', stat = 'mean') + ggtitle("pike")
#> Using all posterior draws for ppc type 'stat' by default.
pp_perch_mean <- pp_check(perch_temp, type = 'stat', stat = 'mean') + ggtitle("perch") + guides(color = "none")
#> Using all posterior draws for ppc type 'stat' by default.
pp_bream_mean <- pp_check(bream_temp, type = 'stat', stat = 'mean') + ggtitle("bream") + guides(color = "none")
#> Using all posterior draws for ppc type 'stat' by default.
pp_cisco_mean <- pp_check(cisco_temp, type = 'stat', stat = 'mean') + ggtitle("cisco") + guides(color = "none")
#> Using all posterior draws for ppc type 'stat' by default.
pp_perch_mean + pp_roach_mean + pp_bream_mean + pp_cisco_mean + pp_pike_mean +
plot_layout(ncol = 2) + plot_annotation(title = "Posterior predictive checks for temperature models (mean) (avg. life temp)")
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Fig. S1: Posterior predictive checks for the temperature models
ggsave("figures/supp/pp_check_temp_mean.png", width = 6.5, height = 6.5, dpi = 600)
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.# Posterior predictive checks
pp_roachb <- pp_check(roach_tempb) + ggtitle("roach") + guides(color = "none")
#> Using 10 posterior draws for ppc type 'dens_overlay' by default.
pp_pikeb <- pp_check(pike_tempb) + ggtitle("pike")
#> Using 10 posterior draws for ppc type 'dens_overlay' by default.
pp_perchb <- pp_check(perch_tempb) + ggtitle("perch") + guides(color = "none")
#> Using 10 posterior draws for ppc type 'dens_overlay' by default.
pp_breamb <- pp_check(bream_tempb) + ggtitle("bream") + guides(color = "none")
#> Using 10 posterior draws for ppc type 'dens_overlay' by default.
pp_ciscob <- pp_check(cisco_tempb) + ggtitle("cisco") + guides(color = "none")
#> Using 10 posterior draws for ppc type 'dens_overlay' by default.
pp_perchb + pp_roachb + pp_breamb + pp_ciscob + pp_pikeb +
plot_layout(ncol = 2) + plot_annotation(title = "Posterior predictive checks for temperature models (birth temp)")
Fig. S2: Posterior predictive checks for the temperature models with birth year temperature
ggsave("figures/supp/pp_check_birth_temp.png", width = 6.5, height = 6.5, dpi = 600)
pp_roach_meanb <- pp_check(roach_tempb, type = 'stat', stat = 'mean') + ggtitle("roach") + guides(color = "none")
#> Using all posterior draws for ppc type 'stat' by default.
pp_pike_meanb <- pp_check(pike_tempb, type = 'stat', stat = 'mean') + ggtitle("pike")
#> Using all posterior draws for ppc type 'stat' by default.
pp_perch_meanb <- pp_check(perch_tempb, type = 'stat', stat = 'mean') + ggtitle("perch") + guides(color = "none")
#> Using all posterior draws for ppc type 'stat' by default.
pp_bream_meanb <- pp_check(bream_tempb, type = 'stat', stat = 'mean') + ggtitle("bream") + guides(color = "none")
#> Using all posterior draws for ppc type 'stat' by default.
pp_cisco_meanb <- pp_check(cisco_tempb, type = 'stat', stat = 'mean') + ggtitle("cisco") + guides(color = "none")
#> Using all posterior draws for ppc type 'stat' by default.
pp_perch_meanb + pp_roach_meanb + pp_bream_meanb + pp_cisco_meanb + pp_pike_meanb +
plot_layout(ncol = 2) + plot_annotation(title = "Posterior predictive checks for temperature models (mean) (birth temp)")
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Fig. S2: Posterior predictive checks for the temperature models with birth year temperature
ggsave("figures/supp/pp_check_birth_temp_mean.png", width = 6.5, height = 6.5, dpi = 600)
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.Chain convergence
pal_diag <- rev(brewer.pal(n = 4, name = "Dark2"))
pars_monitor = c("b_Intercept", "b_age_ct", "b_lifetime_avg_temp_ct", "b_age_sq",
"b_age_ct:lifetime_avg_temp_ct", "sd_month__Intercept", "sd_cohort__Intercept",
"sd_year__Intercept", "sigma")
# Chain convergence
posterior_roach_temp <- as.array(roach_temp)
posterior_pike_temp <- as.array(pike_temp)
posterior_perch_temp <- as.array(perch_temp)
posterior_cisco_temp <- as.array(bream_temp)
posterior_bream_temp <- as.array(bream_temp)
# Roach
mcmc_trace(posterior_roach_temp,
pars = pars_monitor,
facet_args = list(ncol = 2, strip.position = "left")) +
theme(text = element_text(size = 12),
strip.text = element_text(size = 6),
legend.position = "top") +
scale_color_manual(values = alpha(pal_diag, alpha = 0.8)) +
ggtitle("Roach temperature model")
#> Scale for 'colour' is already present. Adding another scale for 'colour',
#> which will replace the existing scale.
ggsave("figures/supp/roach_temp_trace.png", width = 6.5, height = 6.5, dpi = 600)
# Pike
mcmc_trace(posterior_pike_temp,
pars = pars_monitor,
facet_args = list(ncol = 2, strip.position = "left")) +
theme(text = element_text(size = 12),
strip.text = element_text(size = 6),
legend.position = "top") +
scale_color_manual(values = alpha(pal_diag, alpha = 0.8)) +
ggtitle("Pike temperature model")
#> Scale for 'colour' is already present. Adding another scale for 'colour',
#> which will replace the existing scale.
ggsave("figures/supp/pike_temp_trace.png", width = 6.5, height = 6.5, dpi = 600)
# Perch
mcmc_trace(posterior_perch_temp,
pars = pars_monitor,
facet_args = list(ncol = 2, strip.position = "left")) +
theme(text = element_text(size = 12),
strip.text = element_text(size = 6),
legend.position = "top") +
scale_color_manual(values = alpha(pal_diag, alpha = 0.8)) +
ggtitle("Perch temperature model")
#> Scale for 'colour' is already present. Adding another scale for 'colour',
#> which will replace the existing scale.
ggsave("figures/supp/perch_temp_trace.png", width = 6.5, height = 6.5, dpi = 600)
# Bream
mcmc_trace(posterior_bream_temp,
pars = pars_monitor,
facet_args = list(ncol = 2, strip.position = "left")) +
theme(text = element_text(size = 12),
strip.text = element_text(size = 6),
legend.position = "top") +
scale_color_manual(values = alpha(pal_diag, alpha = 0.8)) +
ggtitle("Bream temperature model")
#> Scale for 'colour' is already present. Adding another scale for 'colour',
#> which will replace the existing scale.
ggsave("figures/supp/bream_temp_trace.png", width = 6.5, height = 6.5, dpi = 600)
# Cisco
mcmc_trace(posterior_cisco_temp,
pars = pars_monitor,
facet_args = list(ncol = 2, strip.position = "left")) +
theme(text = element_text(size = 12),
strip.text = element_text(size = 6),
legend.position = "top") +
scale_color_manual(values = alpha(pal_diag, alpha = 0.8)) +
ggtitle("Cisco temperature model")
#> Scale for 'colour' is already present. Adding another scale for 'colour',
#> which will replace the existing scale.
ggsave("figures/supp/cisco_temp_trace.png", width = 6.5, height = 6.5, dpi = 600)pal_diag <- rev(brewer.pal(n = 4, name = "Dark2"))
pars_monitor = c("b_Intercept", "b_age_ct", "b_birth_temp_ct", "b_age_sq",
"b_age_ct:birth_temp_ct", "sd_month__Intercept", "sd_cohort__Intercept",
"sd_year__Intercept", "sigma")
# Chain convergence
posterior_roach_tempb <- as.array(roach_tempb)
posterior_pike_tempb <- as.array(pike_tempb)
posterior_perch_tempb <- as.array(perch_tempb)
posterior_cisco_tempb <- as.array(bream_tempb)
posterior_bream_tempb <- as.array(bream_tempb)
# Roach
mcmc_trace(posterior_roach_tempb,
pars = pars_monitor,
facet_args = list(ncol = 2, strip.position = "left")) +
theme(text = element_text(size = 12),
strip.text = element_text(size = 6),
legend.position = "top") +
scale_color_manual(values = alpha(pal_diag, alpha = 0.8)) +
ggtitle("Roach temperature model")
#> Scale for 'colour' is already present. Adding another scale for 'colour',
#> which will replace the existing scale.
ggsave("figures/supp/roach_b_temp_trace.png", width = 6.5, height = 6.5, dpi = 600)
# Pike
mcmc_trace(posterior_pike_tempb,
pars = pars_monitor,
facet_args = list(ncol = 2, strip.position = "left")) +
theme(text = element_text(size = 12),
strip.text = element_text(size = 6),
legend.position = "top") +
scale_color_manual(values = alpha(pal_diag, alpha = 0.8)) +
ggtitle("Pike temperature model")
#> Scale for 'colour' is already present. Adding another scale for 'colour',
#> which will replace the existing scale.
ggsave("figures/supp/pike_b_temp_trace.png", width = 6.5, height = 6.5, dpi = 600)
# Perch
mcmc_trace(posterior_perch_tempb,
pars = pars_monitor,
facet_args = list(ncol = 2, strip.position = "left")) +
theme(text = element_text(size = 12),
strip.text = element_text(size = 6),
legend.position = "top") +
scale_color_manual(values = alpha(pal_diag, alpha = 0.8)) +
ggtitle("Perch temperature model")
#> Scale for 'colour' is already present. Adding another scale for 'colour',
#> which will replace the existing scale.
ggsave("figures/supp/perch_b_temp_trace.png", width = 6.5, height = 6.5, dpi = 600)
# Bream
mcmc_trace(posterior_bream_tempb,
pars = pars_monitor,
facet_args = list(ncol = 2, strip.position = "left")) +
theme(text = element_text(size = 12),
strip.text = element_text(size = 6),
legend.position = "top") +
scale_color_manual(values = alpha(pal_diag, alpha = 0.8)) +
ggtitle("Bream temperature model")
#> Scale for 'colour' is already present. Adding another scale for 'colour',
#> which will replace the existing scale.
ggsave("figures/supp/bream_b_temp_trace.png", width = 6.5, height = 6.5, dpi = 600)
# Cisco
mcmc_trace(posterior_cisco_tempb,
pars = pars_monitor,
facet_args = list(ncol = 2, strip.position = "left")) +
theme(text = element_text(size = 12),
strip.text = element_text(size = 6),
legend.position = "top") +
scale_color_manual(values = alpha(pal_diag, alpha = 0.8)) +
ggtitle("Cisco temperature model")
#> Scale for 'colour' is already present. Adding another scale for 'colour',
#> which will replace the existing scale.
ggsave("figures/supp/cisco_b_temp_trace.png", width = 6.5, height = 6.5, dpi = 600)knitr::knit_exit()