Bayesian MAR(1) models and structural stability
Matthew A. Barbour
2022-01-24
Last updated: 2022-01-24
Checks: 7 0
Knit directory: genes-to-foodweb-stability/
This reproducible R Markdown analysis was created with workflowr (version 1.6.2). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.
Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.
Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.
The command set.seed(20200205) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.
Great job! Recording the operating system, R version, and package versions is critical for reproducibility.
Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.
Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.
Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.
The results in this page were generated with repository version 18d1722. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.
Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:
Ignored files:
Ignored: .Rhistory
Ignored: .Rproj.user/
Ignored: code/.Rhistory
Ignored: output/.Rapp.history
Untracked files:
Untracked: output/all.mar1.brm.adj.rds
Untracked: output/all.mar1.brm.unadj.ar2.lag.rds
Untracked: output/all.mar1.brm.unadj.noBRBRonLYER.rds
Untracked: output/all.mar1.brm.unadj.rds
Untracked: output/all.mar1.brm.unadj.xAOP2.rds
Untracked: output/initial.mar1.brm.adj.rds
Untracked: output/initial.mar1.brm.unadj.rds
Unstaged changes:
Modified: figures/AOP2-growth-no-insects.pdf
Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.
These are the previous versions of the repository in which changes were made to the R Markdown (analysis/structural-stability.Rmd) and HTML (docs/structural-stability.html) files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view the files as they were in that past version.
| File | Version | Author | Date | Message |
|---|---|---|---|---|
| html | c802852 | mabarbour | 2021-06-24 | Build website |
| Rmd | a054471 | mabarbour | 2021-06-24 | Update location of .rds files to Release v2.1 |
| Rmd | f467085 | mabarbour | 2021-06-24 | Revise analysis to address reviewer’s comments |
| Rmd | 399b240 | mabarbour | 2020-10-11 | remove old versions of analysis |
| Rmd | c1aaf8e | mabarbour | 2020-10-11 | remove old versions of analysis |
| html | 51fc18b | mabarbour | 2020-06-23 | Build site. |
| Rmd | 86116c8 | mabarbour | 2020-06-23 | bioRxiv version of code and data. |
| html | 86116c8 | mabarbour | 2020-06-23 | bioRxiv version of code and data. |
General priors
Intrinsic growth rates
# from Jahan et al. 2014, Journal of Insect Science
# Table 4 lambda (finite rate of increase, discrete time analogue of intrinsic growth rate)
# calculated on a per-day basis and not logged. This is why I multiply by 7 and then take the natural logarithm
Jahan.r.BRBR <- log(c(1.42, 1.36, 1.32, 1.35, 1.40, 1.33, 1.38, 1.37) * 7)
mean(Jahan.r.BRBR) # 2.26[1] 2.257713sd(Jahan.r.BRBR) # 0.02[1] 0.02468356# visualize prior
hist(rnorm(1000, mean(Jahan.r.BRBR), sd = 1))
prior.r.BRBR <- "normal(2.26, 1)"
# from Taghizadeh 2019, J. Agr. Sci. Tech.
# Table 2 lambda (finite rate of increase, discrete time analogue of intrinsic growth rate)
# calculated on a per-day basis and not logged. This is why I multiply by 7 and then take the natural logarithm
Tag.r.LYER <- log(c(1.35, 1.30, 1.26, 1.23) * 7)
mean(Tag.r.LYER) # 2.20[1] 2.196059sd(Tag.r.LYER) # 0.04[1] 0.04028153# visualize prior
hist(rnorm(1000, mean(Tag.r.LYER), sd = 1))
prior.r.LYER <- "normal(2.20, 1)"
# random effects prior based on variance among cultivars
# I'm just going to use this for all of them, including parasitoids
mean.r.sd <- mean(c(sd(Jahan.r.BRBR), sd(Tag.r.LYER)))
# visualize prior
hist(rnorm(1000, mean = mean.r.sd, sd = 0.5))
# not changing for now
prior.random.effects <- "normal(0.03, 0.5)" # mean of BRBR and LYER
# I don't have a great sense for the growth rate of the parasitoid, other than that it should be negative
# this seems like a reasonable starting point
# visualize prior
hist(rnorm(1000, mean = -1.5, sd = 1))
prior.r.Ptoid <- "normal(-1.5, 1)"Intra- and interspecific interactions
I assume they are weak, i.e. much less than \(|1|\). I also assume that all species exhibit intraspecific competition, aphids have negative interspecific effects with each other, but positive interspecific effects on the parasitoid. I also assume parasitoids have negative interspecific effects on the aphids.
## intraspecific, 1 = no density dependence. I would prefer to specify an offset first, so that 0 = no density dependence, like the other coefs, but I can't use the offset if I incorporate measurement error :-(
# visualize prior
hist(rnorm(1000, mean = 0.9, sd = 0.5))
prior.intra.BRBR <- "normal(0.9, 0.5)"
prior.intra.LYER <- "normal(0.9, 0.5)"
prior.intra.Ptoid <- "normal(0.9, 0.5)"
## negative interspecific, 0 = no interaction
# visualize prior
hist(rnorm(1000, mean = -0.1, sd = 0.5))
# most of these values are less than 1, which
# is indicative of saturating effects
prior.LYERonBRBR <- "normal(-0.1, 0.5)"
prior.PtoidonBRBR <- "normal(-0.1, 0.5)"
prior.BRBRonLYER <- "normal(-0.1, 0.5)"
prior.PtoidonLYER <- "normal(-0.1, 0.5)"
## positive interspecific
# visualize prior
hist(rnorm(1000, mean = 0.1, sd = 0.5))
# most of these values are less than 1, which
# is indicative of saturating effects
prior.BRBRonPtoid <- "normal(0.1, 0.5)"
prior.LYERonPtoid <- "normal(0.1, 0.5)"AOP2 effects
It was unclear to me a priori exactly how allelic differences at AOP2 would affect species’ growth rates or intra- and interspecific interactions. Therefore, I assumed these effects on specific rates could be positive or negative, but would likely be between -1 and 1 (i.e., not ridiculously large).
prior.AOP2 <- "normal(0, 0.5)"Temperature effects
As with AOP2 it was unclear to me a priori how warming would affect species’ growth rates or intra- and interspecific interactions; therefore, I used the same prior as for AOP2.
prior.temp <- "normal(0, 0.5)"Biomass effects
For both aphids, I thought that increasing biomass would increase their intrinsic growth rates, but only weakly, because I didn’t expect biomass to be limiting.
# visualize prior
hist(rnorm(1000, mean = 0.1, sd = 0.5))
prior.AphidBiomass <- "normal(0.1, 0.5)"For the parasitoid, it was unclear to me whether increasing biomass would have positive or negative effects. I could imagine both, as increasing biomass may increase the search effort of parasitoids, resulting in a negative effect on their growth rate. Alternatively, more biomass may increase the quality of aphids, which could increase the parasitoid’s growth rate. Therefore, I specified a normal prior with mean = 0 and SD = 0.5, so that most of the distribution lied between -1 and 1 (i.e. saturating negative or positive effects).
# visualize prior
hist(rnorm(1000, mean = 0, sd = 0.5))
prior.PtoidBiomass <- "normal(0, 0.5)"Model including all species
Formula
This follows equation 1 in the Supplementary Material.
# BRBR
all.BRBR.bf <- bf(log(BRBR_t1) ~ 0 + Intercept + (me(logBRBR_t, se_logBRBRt) + me(logLYER_t, se_logLYERt) + me(logPtoid_t, se_logPtoidt)) + aop2_genotypes + AOP2_genotypes + temp + log(Biomass_g_t1) + (1|p|Cage))
# LYER
all.LYER.bf <- bf(log(LYER_t1) ~ 0 + Intercept + (me(logLYER_t, se_logLYERt) + me(logBRBR_t, se_logBRBRt) + me(logPtoid_t, se_logPtoidt)) + aop2_genotypes + AOP2_genotypes + temp + log(Biomass_g_t1) + (1|p|Cage))
# Ptoid
all.Ptoid.bf <- bf(log(Ptoid_t1) ~ 0 + Intercept + me(logPtoid_t, se_logPtoidt) + me(logBRBR_t, se_logBRBRt) + me(logLYER_t, se_logLYERt) + aop2_genotypes + AOP2_genotypes + temp + log(Biomass_g_t1) + (1|p|Cage))Priors summarized in Table S3
all.mv.priors <- c(
# aop2 and AOP2 effects
set_prior(prior.AOP2, class = "b", coef = "aop2_genotypes", resp = "logBRBRt1"),
set_prior(prior.AOP2, class = "b", coef = "AOP2_genotypes", resp = "logBRBRt1"),
set_prior(prior.AOP2, class = "b", coef = "aop2_genotypes", resp = "logLYERt1"),
set_prior(prior.AOP2, class = "b", coef = "AOP2_genotypes", resp = "logLYERt1"),
set_prior(prior.AOP2, class = "b", coef = "aop2_genotypes", resp = "logPtoidt1"),
set_prior(prior.AOP2, class = "b", coef = "AOP2_genotypes", resp = "logPtoidt1"),
# biomass effects
set_prior(prior.AphidBiomass, class = "b", coef = "logBiomass_g_t1", resp = "logBRBRt1"),
set_prior(prior.AphidBiomass, class = "b", coef = "logBiomass_g_t1", resp = "logLYERt1"),
set_prior(prior.PtoidBiomass, class = "b", coef = "logBiomass_g_t1", resp = "logPtoidt1"),
# baseline growth rates
set_prior(prior.r.BRBR, class = "b", coef = "Intercept", resp = "logBRBRt1"),
set_prior(prior.r.LYER, class = "b", coef = "Intercept", resp = "logLYERt1"),
set_prior(prior.r.Ptoid, class = "b", coef = "Intercept", resp = "logPtoidt1"),
# intraspecific effects
set_prior(prior.intra.BRBR, class = "b", coef = "melogBRBR_tse_logBRBRt", resp = "logBRBRt1"),
set_prior(prior.intra.LYER, class = "b", coef = "melogLYER_tse_logLYERt", resp = "logLYERt1"),
set_prior(prior.intra.Ptoid, class = "b", coef = "melogPtoid_tse_logPtoidt", resp = "logPtoidt1"),
# negative interspecific effects
set_prior(prior.LYERonBRBR, class = "b", coef = "melogLYER_tse_logLYERt", resp = "logBRBRt1"),
set_prior(prior.BRBRonLYER, class = "b", coef = "melogBRBR_tse_logBRBRt", resp = "logLYERt1"),
set_prior(prior.PtoidonBRBR, class = "b", coef = "melogPtoid_tse_logPtoidt", resp = "logBRBRt1"),
set_prior(prior.PtoidonLYER, class = "b", coef = "melogPtoid_tse_logPtoidt", resp = "logLYERt1"),
# positive interspecific effects
set_prior(prior.BRBRonPtoid, class = "b", coef = "melogBRBR_tse_logBRBRt", resp = "logPtoidt1"),
set_prior(prior.LYERonPtoid, class = "b", coef = "melogLYER_tse_logLYERt", resp = "logPtoidt1"),
# temp effects
set_prior(prior.temp, class = "b", coef = "temp", resp = "logBRBRt1"),
set_prior(prior.temp, class = "b", coef = "temp", resp = "logLYERt1"),
set_prior(prior.temp, class = "b", coef = "temp", resp = "logPtoidt1"),
# random effects
set_prior(prior.random.effects, class = "sd", resp = "logBRBRt1"),
set_prior(prior.random.effects, class = "sd", resp = "logLYERt1"),
set_prior(prior.random.effects, class = "sd", resp = "logPtoidt1"))Fit Model
all.mar1.brm.unadj <- brm(
data = full_df,
formula = mvbf(all.BRBR.bf, all.LYER.bf, all.Ptoid.bf) + set_rescor(TRUE),
iter = 5000,
save_pars = save_pars(latent = TRUE, all = TRUE),
prior = all.mv.priors,
file = "output/all.mar1.brm.unadj.rds")
all.mar1.brm.unadj Family: MV(gaussian, gaussian, gaussian)
Links: mu = identity; sigma = identity
mu = identity; sigma = identity
mu = identity; sigma = identity
Formula: log(BRBR_t1) ~ 0 + Intercept + (me(logBRBR_t, se_logBRBRt) + me(logLYER_t, se_logLYERt) + me(logPtoid_t, se_logPtoidt)) + aop2_genotypes + AOP2_genotypes + temp + log(Biomass_g_t1) + (1 | p | Cage)
log(LYER_t1) ~ 0 + Intercept + (me(logLYER_t, se_logLYERt) + me(logBRBR_t, se_logBRBRt) + me(logPtoid_t, se_logPtoidt)) + aop2_genotypes + AOP2_genotypes + temp + log(Biomass_g_t1) + (1 | p | Cage)
log(Ptoid_t1) ~ 0 + Intercept + me(logPtoid_t, se_logPtoidt) + me(logBRBR_t, se_logBRBRt) + me(logLYER_t, se_logLYERt) + aop2_genotypes + AOP2_genotypes + temp + log(Biomass_g_t1) + (1 | p | Cage)
Data: full_df (Number of observations: 264)
Draws: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
total post-warmup draws = 10000
Group-Level Effects:
~Cage (Number of levels: 60)
Estimate Est.Error l-95% CI
sd(logBRBRt1_Intercept) 0.18 0.12 0.01
sd(logLYERt1_Intercept) 0.17 0.10 0.01
sd(logPtoidt1_Intercept) 0.07 0.06 0.00
cor(logBRBRt1_Intercept,logLYERt1_Intercept) -0.10 0.48 -0.89
cor(logBRBRt1_Intercept,logPtoidt1_Intercept) 0.08 0.50 -0.85
cor(logLYERt1_Intercept,logPtoidt1_Intercept) 0.01 0.50 -0.88
u-95% CI Rhat Bulk_ESS Tail_ESS
sd(logBRBRt1_Intercept) 0.46 1.00 2549 4640
sd(logLYERt1_Intercept) 0.39 1.00 2321 4183
sd(logPtoidt1_Intercept) 0.21 1.00 6017 4961
cor(logBRBRt1_Intercept,logLYERt1_Intercept) 0.81 1.00 4528 6645
cor(logBRBRt1_Intercept,logPtoidt1_Intercept) 0.91 1.00 12875 7451
cor(logLYERt1_Intercept,logPtoidt1_Intercept) 0.88 1.00 11683 8543
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat
logBRBRt1_Intercept 2.15 0.59 0.99 3.31 1.00
logBRBRt1_aop2_genotypes 0.03 0.13 -0.22 0.29 1.00
logBRBRt1_AOP2_genotypes 0.00 0.12 -0.23 0.23 1.00
logBRBRt1_temp -0.45 0.08 -0.61 -0.31 1.00
logBRBRt1_logBiomass_g_t1 0.02 0.25 -0.46 0.50 1.00
logLYERt1_Intercept 3.08 0.51 2.09 4.08 1.00
logLYERt1_aop2_genotypes 0.19 0.11 -0.01 0.40 1.00
logLYERt1_AOP2_genotypes -0.08 0.10 -0.27 0.11 1.00
logLYERt1_temp 0.00 0.06 -0.11 0.13 1.00
logLYERt1_logBiomass_g_t1 0.05 0.21 -0.35 0.45 1.00
logPtoidt1_Intercept -1.92 0.55 -2.98 -0.84 1.00
logPtoidt1_aop2_genotypes 0.22 0.11 0.01 0.44 1.00
logPtoidt1_AOP2_genotypes -0.11 0.10 -0.30 0.09 1.00
logPtoidt1_temp -0.16 0.06 -0.27 -0.04 1.00
logPtoidt1_logBiomass_g_t1 -1.12 0.22 -1.55 -0.68 1.00
logBRBRt1_melogBRBR_tse_logBRBRt 0.72 0.09 0.54 0.91 1.00
logBRBRt1_melogLYER_tse_logLYERt 0.04 0.12 -0.19 0.27 1.00
logBRBRt1_melogPtoid_tse_logPtoidt -0.73 0.06 -0.86 -0.61 1.00
logLYERt1_melogLYER_tse_logLYERt 0.34 0.10 0.14 0.54 1.00
logLYERt1_melogBRBR_tse_logBRBRt 0.24 0.08 0.09 0.40 1.00
logLYERt1_melogPtoid_tse_logPtoidt -0.68 0.05 -0.78 -0.59 1.00
logPtoidt1_melogPtoid_tse_logPtoidt 0.99 0.05 0.88 1.09 1.00
logPtoidt1_melogBRBR_tse_logBRBRt 0.07 0.07 -0.08 0.21 1.00
logPtoidt1_melogLYER_tse_logLYERt 0.45 0.10 0.26 0.64 1.00
Bulk_ESS Tail_ESS
logBRBRt1_Intercept 10107 7405
logBRBRt1_aop2_genotypes 13546 8138
logBRBRt1_AOP2_genotypes 13125 8587
logBRBRt1_temp 6602 7763
logBRBRt1_logBiomass_g_t1 11071 8333
logLYERt1_Intercept 9645 7397
logLYERt1_aop2_genotypes 12156 7630
logLYERt1_AOP2_genotypes 12192 7880
logLYERt1_temp 7232 6676
logLYERt1_logBiomass_g_t1 11013 8015
logPtoidt1_Intercept 9133 6787
logPtoidt1_aop2_genotypes 13611 7930
logPtoidt1_AOP2_genotypes 13266 8190
logPtoidt1_temp 10937 8428
logPtoidt1_logBiomass_g_t1 10686 8894
logBRBRt1_melogBRBR_tse_logBRBRt 6297 6948
logBRBRt1_melogLYER_tse_logLYERt 6038 6329
logBRBRt1_melogPtoid_tse_logPtoidt 11112 8659
logLYERt1_melogLYER_tse_logLYERt 5820 7142
logLYERt1_melogBRBR_tse_logBRBRt 5761 6059
logLYERt1_melogPtoid_tse_logPtoidt 11829 7963
logPtoidt1_melogPtoid_tse_logPtoidt 13006 8341
logPtoidt1_melogBRBR_tse_logBRBRt 9621 7257
logPtoidt1_melogLYER_tse_logLYERt 9713 7690
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma_logBRBRt1 1.23 0.06 1.13 1.35 1.00 12496 7477
sigma_logLYERt1 0.95 0.05 0.87 1.05 1.00 10246 7704
sigma_logPtoidt1 1.02 0.05 0.94 1.12 1.00 15968 8374
Residual Correlations:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
rescor(logBRBRt1,logLYERt1) 0.48 0.05 0.38 0.58 1.00 9085
rescor(logBRBRt1,logPtoidt1) -0.22 0.06 -0.34 -0.10 1.00 14309
rescor(logLYERt1,logPtoidt1) -0.06 0.06 -0.18 0.07 1.00 13786
Tail_ESS
rescor(logBRBRt1,logLYERt1) 8516
rescor(logBRBRt1,logPtoidt1) 8007
rescor(logLYERt1,logPtoidt1) 7635
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).Carrying capacity check
# BRBR carrying capacity
fixef(all.mar1.brm.unadj)["logBRBRt1_Intercept","Estimate"] / (1 - fixef(all.mar1.brm.unadj)["logBRBRt1_melogBRBR_tse_logBRBRt","Estimate"]) [1] 7.81221max(log(full_df$BRBR_t1)) # max observed in experiment[1] 6.727432# LYER carrying capacity
fixef(all.mar1.brm.unadj)["logLYERt1_Intercept","Estimate"] / (1 - fixef(all.mar1.brm.unadj)["logLYERt1_melogLYER_tse_logLYERt","Estimate"]) # too low[1] 4.675535Reproduce Fig. S5
# intrinsic growth rates
r_plots <- mcmc_intervals_data(all.mar1.brm.unadj, pars = c("b_logBRBRt1_Intercept","b_logLYERt1_Intercept","b_logPtoidt1_Intercept"), prob = 0.66, prob_outer = 0.95) %>%
mutate(parameter = factor(parameter, levels = c("b_logPtoidt1_Intercept","b_logLYERt1_Intercept","b_logBRBRt1_Intercept"))) %>%
ggplot(aes(x = parameter, y = m)) +
geom_hline(yintercept = 0, linetype = "dotted") +
geom_linerange(aes(ymin = ll, ymax = hh), size = 0.5) +
geom_linerange(aes(ymin = l, ymax = h), size = 1) +
geom_point(size = 1.5) +
ylab(expression(paste("Intrinsic growth rates (", italic("r"),")"))) +
scale_x_discrete(name = "", labels = c("D","L","B")) +
coord_flip() +
theme_cowplot(font_size = 8)
# AOP2 effect on intrinsic growth rates
delta_r_plots <- mcmc_intervals_data(all.mar1.brm.unadj, regex_pars = c("_aop2_genotypes","_AOP2_genotypes")) %>%
separate(col = "parameter", into = c("b","species","allele","genotypes")) %>%
mutate(species = factor(species, levels = c("logPtoidt1","logLYERt1","logBRBRt1"))) %>%
ggplot(aes(x = species, y = m, color = allele)) +
geom_hline(yintercept = 0, linetype = "dotted", color = "grey") +
geom_point(position = position_dodge(width = 0.4), size = 1.5) +
geom_linerange(aes(ymin = l, ymax = h), size = 1, position = position_dodge(width = 0.4)) +
geom_linerange(aes(ymin = ll, ymax = hh), size = 0.5, position = position_dodge(width = 0.4)) +
scale_color_grey(name = "", labels = c("AOP2\u2013", "AOP2+")) + # viridis_d
ylab(expression(paste("Allele effect (", Delta*"r", italic("G"),")"))) +
scale_x_discrete(name = "", labels = c("D","L","B")) + #c(expression(italic("D. rapae")), expression(italic("L. erysimi")), expression(italic("B. brassicae")))) +
theme_cowplot(font_size = 8) +
coord_flip()
# Interaction matrix
b_data <- mcmc_intervals_data(
all.mar1.brm.unadj,
pars = c("bsp_logBRBRt1_melogBRBR_tse_logBRBRt",
"bsp_logBRBRt1_melogLYER_tse_logLYERt",
"bsp_logBRBRt1_melogPtoid_tse_logPtoidt",
"bsp_logLYERt1_melogBRBR_tse_logBRBRt",
"bsp_logLYERt1_melogLYER_tse_logLYERt",
"bsp_logLYERt1_melogPtoid_tse_logPtoidt",
"bsp_logPtoidt1_melogBRBR_tse_logBRBRt",
"bsp_logPtoidt1_melogLYER_tse_logLYERt",
"bsp_logPtoidt1_melogPtoid_tse_logPtoidt"),
prob = 0.66, prob_outer = 0.95) %>%
mutate(parameter = factor(parameter,
labels = c("B \u2b62 B","L \u2b62 B","D \u2b62 B","B \u2b62 L","L \u2b62 L","D \u2b62 L","B \u2b62 D","L \u2b62 D","D \u2b62 D")),
baseline = c(1,0,0,
0,1,0,
0,0,1))
b_plots <- ggplot(b_data, aes(x = parameter, y = m, color = factor(baseline))) +
geom_hline(aes(yintercept = baseline, color = factor(baseline)), linetype = "dotted") +
geom_linerange(aes(ymin = ll, ymax = hh), size = 0.5) +
geom_linerange(aes(ymin = l, ymax = h), size = 1) +
geom_point(size = 1.5) +
ylab(expression(paste("Interactions (", italic("b"),")"))) +
scale_x_discrete(name = "", limits = rev) +
theme_cowplot(font_size = 8) +
coord_flip() +
scale_color_grey(name = "", labels = c("Inter","Intra")) # viridis_d
# biomass effect on intrinsic growth rates
p_plots <- mcmc_intervals_data(all.mar1.brm.unadj, regex_pars = "Biomass_g_t1$", prob = 0.66, prob_outer = 0.95) %>%
mutate(parameter = factor(parameter, levels = c("b_logPtoidt1_logBiomass_g_t1","b_logLYERt1_logBiomass_g_t1","b_logBRBRt1_logBiomass_g_t1"))) %>%
ggplot(aes(x = parameter, y = m)) +
geom_hline(yintercept = 0, linetype = "dotted") +
geom_linerange(aes(ymin = ll, ymax = hh), size = 0.5) +
geom_linerange(aes(ymin = l, ymax = h), size = 1) +
geom_point(size = 1.5) +
ylab(expression(paste("Plant biomass effect (", italic("p"),")"))) +
scale_x_discrete(name = "", labels = c("D","L","B")) +
coord_flip() +
theme_cowplot(font_size = 8)
# temperature effect on intrinsic growth rates
t_plots <- mcmc_intervals_data(all.mar1.brm.unadj, regex_pars = "temp$", prob = 0.66, prob_outer = 0.95) %>%
mutate(parameter = factor(parameter, levels = c("b_logPtoidt1_temp","b_logLYERt1_temp","b_logBRBRt1_temp"))) %>%
ggplot(aes(x = parameter, y = m)) +
geom_hline(yintercept = 0, linetype = "dotted") +
geom_linerange(aes(ymin = ll, ymax = hh), size = 0.5) +
geom_linerange(aes(ymin = l, ymax = h), size = 1) +
geom_point(size = 1.5) +
ylab(expression(paste("Temperature effect (", Delta*"r", italic("T"),")"))) +
scale_x_discrete(name = "", labels = c("D","L","B")) +
coord_flip() +
theme_cowplot(font_size = 8)
# combine plots
MAR1_parameter_plot <- plot_grid(r_plots, delta_r_plots, t_plots, b_plots, p_plots, nrow = 3, ncol = 2, align = 'v', axis = 'l', rel_widths = c(0.8,1))
x11(); MAR1_parameter_plot
#save_plot(filename = "figures/MAR1-parameter-plot.pdf", plot = MAR1_parameter_plot, device=cairo_pdf, base_asp = 1.1)
| Version | Author | Date |
|---|---|---|
| c802852 | mabarbour | 2021-06-24 |
Structural stability
stability_all.mar1.brm.unadj <- aop2_vs_AOP2_posterior_samples_unadj(all.mar1.brm.unadj, n.geno = 2, temp.value = 0, logbiomass.value = 0)
stability_all.mar1.brm.unadj$aop2_SS_LP_BayesP # clearly stabilizing[1] 0.9966stability_all.mar1.brm.unadj$aop2_SS_LP_effect[1] 13.0206stability_all.mar1.brm.unadj$aop2_SS_LP_95CI 2.5% 97.5%
3.873221 23.478849 Aphid and parasitoid intrinsic growth rates
stability_all.mar1.brm.unadj$aop2_r_LYER_effect[1] 0.5389944stability_all.mar1.brm.unadj$aop2_r_LYER_95CI 2.5% 97.5%
-0.04038828 1.12561843 stability_all.mar1.brm.unadj$aop2_r_Ptoid_effect[1] 0.6588181stability_all.mar1.brm.unadj$aop2_r_Ptoid_95CI 2.5% 97.5%
0.0859789 1.2534024 stability_all.mar1.brm.unadj_n.geno1 <- aop2_vs_AOP2_posterior_samples_unadj(all.mar1.brm.unadj, n.geno = 1, temp.value = 0, logbiomass.value = 0)
stability_all.mar1.brm.unadj_n.geno1$aop2_r_LYER_effect[1] 0.2694972stability_all.mar1.brm.unadj_n.geno1$aop2_r_LYER_95CI 2.5% 97.5%
-0.02019414 0.56280922 Reproduce Fig. S6 and S7
# merge full data set and that with remaining
complete_df <- bind_rows(aphids_only_df, full_df, LP_df, L_df, P_df) %>%
arrange(Cage, Week)
predict_data <- complete_df
# get predictions
BRBR_predict_df <- predict_data %>% filter(BRBR_Survival == 1) %>% mutate(Abundance = BRBR_t1)
predict_BRBRt1 <- bind_cols(BRBR_predict_df, data.frame(predict(all.mar1.brm.unadj, newdata = BRBR_predict_df, resp = "logBRBRt1"))) %>% mutate(Species = "BRBR")
LYER_predict_df <- predict_data %>% filter(LYER_Survival == 1) %>% mutate(Abundance = LYER_t1)
predict_LYERt1 <- bind_cols(LYER_predict_df, data.frame(predict(all.mar1.brm.unadj, newdata = LYER_predict_df, resp = "logLYERt1"))) %>% mutate(Species = "LYER")
Ptoid_predict_df <- predict_data %>% filter(Mummy_Ptoids_Survival == 1) %>% mutate(Abundance = Ptoid_t1)
predict_Ptoidt1 <- bind_cols(Ptoid_predict_df, data.frame(predict(all.mar1.brm.unadj, newdata = Ptoid_predict_df, resp = "logPtoidt1"))) %>% mutate(Species = "Ptoid")
# combine data
predict_all <- bind_rows(predict_BRBRt1, predict_LYERt1, predict_Ptoidt1) %>%
mutate(Species = factor(Species, levels = c("BRBR","LYER","Ptoid"), labels = c("B. brassicae", "L. erysimi", "D. rapae")))
# get week when first species went extinct (always BRBR at least)
full_df_last_week <- full_df %>%
group_by(Cage) %>%
summarise(last_week = last(Week))
# would be nice to color code facets by number of aop2 genotypes e.g.
# something like? https://stackoverflow.com/questions/19440069/ggplot2-facet-wrap-strip-color-based-on-variable-in-data-set
# 20 C cages
plot_cage_dynamics20C <- filter(predict_all, Cage %in% 1:30) %>%
left_join(., full_df_last_week) %>%
ggplot(aes(x = Week, y = log(Abundance), group = interaction(Cage, Species))) +
geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5, fill = Species), alpha = 0.25) +
geom_line(aes(y = Estimate, color = Species)) +
geom_jitter(aes(color = Species), size = 0.5) +
facet_wrap(~Cage, nrow = 5, ncol = 6) +
scale_color_viridis_d() +
scale_fill_viridis_d() +
theme_cowplot(font_size = 10) +
geom_vline(xintercept = 1.5, linetype = "dotted", color = "grey") +
geom_vline(aes(xintercept = last_week), linetype = "dotted", color = "grey") +
coord_cartesian(ylim = c(0,8)) +
ylab(expression(log(Abundance[t+1])))
x11(); plot_cage_dynamics20C
#ggsave(plot = plot_cage_dynamics20C, filename = "figures/cage-dynamics-20C.pdf", height = 5, width = 6, device=cairo_pdf)
# 23 C cages
plot_cage_dynamics23C <- filter(predict_all, Cage %in% 31:60) %>%
left_join(., full_df_last_week) %>%
ggplot(aes(x = Week, y = log(Abundance), group = interaction(Cage, Species))) +
geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5, fill = Species), alpha = 0.25) +
geom_line(aes(y = Estimate, color = Species)) +
geom_point(aes(color = Species), size = 0.5) +
facet_wrap(~Cage, nrow = 5, ncol = 6) +
scale_color_viridis_d() +
scale_fill_viridis_d() +
theme_cowplot(font_size = 10) +
geom_vline(xintercept = 1.5, linetype = "dotted", color = "grey") +
geom_vline(aes(xintercept = last_week), linetype = "dotted", color = "grey") +
coord_cartesian(ylim = c(0,8)) +
ylab(expression(log(Abundance[t+1])))
| Version | Author | Date |
|---|---|---|
| c802852 | mabarbour | 2021-06-24 |
x11(); plot_cage_dynamics23C
#ggsave(plot = plot_cage_dynamics23C, filename = "figures/cage-dynamics-23C.pdf", height = 5, width = 6, device=cairo_pdf)
| Version | Author | Date |
|---|---|---|
| c802852 | mabarbour | 2021-06-24 |
MAR(1) Bayesian R2 in Table S4
bayes_R2(all.mar1.brm.unadj, newdata = BRBR_predict_df, resp = "logBRBRt1") Estimate Est.Error Q2.5 Q97.5
R2logBRBRt1 0.6459414 0.0240262 0.5964847 0.6889248bayes_R2(all.mar1.brm.unadj, newdata = LYER_predict_df, resp = "logLYERt1") Estimate Est.Error Q2.5 Q97.5
R2logLYERt1 0.4868861 0.03244113 0.4238835 0.5504909bayes_R2(all.mar1.brm.unadj, newdata = Ptoid_predict_df, resp = "logPtoidt1") Estimate Est.Error Q2.5 Q97.5
R2logPtoidt1 0.6301998 0.01178371 0.6018383 0.6475559Reproduce core of Fig. 4
Plot the effect of AOP2 gene on structural stability of three-species food chain.
median_all.mar1.brm.unadj <- aop2_vs_AOP2_median_effects_unadj(all.mar1.brm.unadj, n.geno = 2, temp.value = 0, logbiomass.value = 0)
# get raw data for manually making plot
get_FD.2sp <- FeasibilityDomain2sp(A = list(median_all.mar1.brm.unadj$aop2.mat[2:3,2:3],
median_all.mar1.brm.unadj$AOP2.mat[2:3,2:3]),
r = list(median_all.mar1.brm.unadj$aop2.IGR[2:3],
median_all.mar1.brm.unadj$AOP2.IGR[2:3]),
labels = c("aop2", "AOP2"),
normalize = TRUE) %>%
rename(aop2 = A_ID)
# Draw polygon for feasibility domain
# from: https://stackoverflow.com/questions/12794596/how-fill-part-of-a-circle-using-ggplot2
# define the circle; add a point at the center if the 'pie slice' if the shape is to be filled
circleFun <- function(center=c(0,0), diameter=1, npoints=100, start=0, end=2, filled=TRUE){
tt <- seq(start*pi, end*pi, length.out=npoints)
df <- data.frame(
x = center[1] + diameter / 2 * cos(tt),
y = center[2] + diameter / 2 * sin(tt)
)
if(filled==TRUE) { #add a point at the center so the whole 'pie slice' is filled
df <- rbind(df, center)
}
return(df)
}
## plot figure 4
# alpha_level <- 0.05 # very low so use better looking arrows with keynote
plot_fig_4 <- ggplot(filter(get_FD.2sp, Type == "r"), aes(x = Sp_1, y = Sp_2)) +
# draw intrinsic growth rate vectors
geom_segment(aes(x = 0, y = 0, xend = Sp_1, yend = Sp_2, linetype = aop2), # alpha for manuscript
arrow = arrow(type = 'open', length = unit(0.1,"cm"))) +
# draw critical boundary (remove for final plot after adjusting geom_polygon)
# geom_segment(data = filter(get_FD.2sp, Type == "A")[c(1,3),], # just need one lower bound
# aes(x = 0, y = 0, xend = Sp_1, yend = Sp_2, alpha = aop2),
# linetype = "solid",
# size = 0.5) +
xlab("Aphid growth rate (normalized)") +
ylab("Parasitoid growth rate (normalized)") +
# illustrate circular nature of feasibility domain
coord_cartesian(xlim = c(0,1), ylim = c(-0.9,0), expand = F) +
# scale_alpha_manual(values = c(alpha_level, alpha_level), labels = c("AOP2\u2013","AOP2+"), name = "") + # used for manuscript plot to provide better looking arrows with keynote
scale_linetype_manual(values = c("solid", "dashed"), labels = c("AOP2\u2013","AOP2+"), name = "") +
# adjusted until critical boundary was correct, then removed critical boundary for final plot
geom_polygon(data=circleFun(c(0,0), diameter=2, start=0, end=-0.192, npoints=100, filled=TRUE),
aes(x,y), alpha = 0.1, inherit.aes = F) +
theme_cowplot(font_size = 18, line_size = 1)
x11(); plot_fig_4
# changes for final version
# alpha = aop2 in geom_segment(aes())
# uncomment scale_alpha_manual
# comment out scale_linetype_manual
# remove line for critical boundary
# alpha version saved for keynote # ggsave(plot = plot_fig_4, filename = "figures/keystone-structural-stability-forkeynote.pdf", width = 8, height = 8, units = "in")
| Version | Author | Date |
|---|---|---|
| c802852 | mabarbour | 2021-06-24 |
I then used Keynote to manually add images and text to create the final version presented in Figure 4 of the main text.
Reproduce Fig. S8
# subsample 1/8 of the posterior to make it easier to visualize
rsamp <- sample(1:10000, size = 1000)
# plot
plot_MAR1_posterior_foodchain_AOP2 <- stability_all.mar1.brm.unadj$all.aop2_vs_AOP2_stability.df %>%
#filter(r_Ptoid < 0) %>%
filter(posterior_sample %in% rsamp) %>%
mutate(n.allele = as.numeric(as.factor(aop2_vs_AOP2))) %>%
ggplot(., aes(x = n.allele, y = FeasibilityBoundaryLYER.Ptoid)) +
geom_line(aes(group = posterior_sample), alpha = 0.1) +
stat_summary(fun.y = mean, geom = "line", color = "firebrick1", size = 1) +
stat_summary(fun.y = mean, geom = "point", color = "firebrick1", size = 1.5) +
theme_minimal_hgrid() +
scale_x_continuous(name = "Allele", breaks = c(1,2), labels = c("AOP2+","AOP2\u2013"), expand = c(0.1,0.1)) +
ylab("Normalized angle from critical boundary")
#ggsave(filename = "figures/MAR1-posterior-foodchain-AOP2.pdf", width = 6, height = 5, device=cairo_pdf)Reproduce Fig. S9
The above plot illustrates the effect of AOP2 on the structural stability of the equilibrium abundances of species. I can explore whether our results hold in a non-equilibrium scenario that better characterizes our observational data.
To do this, I look at the the effect of AOP2 across a range of initial conditions for the abundances of LYER and Ptoid. I get this data by simulating community dynamics with the observed effects of AOP2 across a range of initial conditions. I restricted our simulation to 10 time steps, as BRBR went extinct commonly at week 7 (experiment was 17 weeks long). I also set an extinction threshold of 2 individuals.
LP_duration <- 10
LP_threshold <- log(2) # set threshold of two individuals in the populations
res <- 0.1
LP_test_df <- expand.grid(LYER = seq(1, 6, by = res), Ptoid = seq(1, 6, by = res))
## simulate population dynamics and determine which species goes extinct
# aop2
FE_LP_aop2 <- list()
for(i in 1:length(LP_test_df$LYER)){
FE_LP_aop2[[i]] <- first_extinction_2sp(Initial.States = c(LYER = LP_test_df[i,"LYER"], Ptoid = LP_test_df[i,"Ptoid"]),
Interaction.Matrix = median_all.mar1.brm.unadj$aop2.mat[2:3,2:3] + diag(2),
IGR.Vector = median_all.mar1.brm.unadj$aop2.IGR[2:3],
Duration = LP_duration,
threshold = LP_threshold)
}
FE_LP_aop2_df <- bind_cols(LP_test_df, plyr::ldply(FE_LP_aop2)) %>%
mutate(allele = "aop2")
# AOP2
FE_LP_AOP2 <- list()
for(i in 1:length(LP_test_df$LYER)){
FE_LP_AOP2[[i]] <- first_extinction_2sp(Initial.States = c(LYER = LP_test_df[i,"LYER"], Ptoid = LP_test_df[i,"Ptoid"]),
Interaction.Matrix = median_all.mar1.brm.unadj$AOP2.mat[2:3,2:3] + diag(2),
IGR.Vector = median_all.mar1.brm.unadj$AOP2.IGR[2:3],
Duration = LP_duration,
threshold = LP_threshold)
}
FE_LP_AOP2_df <- bind_cols(LP_test_df, plyr::ldply(FE_LP_AOP2)) %>%
mutate(allele = "AOP2")
# get observed data on initial abundances of LYER and Ptoid after BRBR went extinct
cage_type <- LP_df %>%
distinct(Cage, aop2_vs_AOP2) %>%
mutate(allele = ifelse(aop2_vs_AOP2 > 0, 1,
ifelse(aop2_vs_AOP2 < 0, -1, NA)))
LP_actual_df <- LP_df %>%
group_by(Cage) %>% # rich,
summarise_at(vars(LYER_t, Ptoid_t), list(first = first)) %>%
ungroup() %>%
mutate(log_LYER_t_first = log(LYER_t_first),
log_Ptoid_t_first = log(Ptoid_t_first)) %>%
as.data.frame() %>%
left_join(., cage_type) %>%
drop_na() %>% # don't use this in order to plot all cages
mutate(allele = factor(allele, levels = c(1,-1), labels = c("AOP2\u2013","AOP2+")))The graph below shows a couple of useful things. First, our predictions hold for outside of equilibrium. That is, there is a greater likelihood of LYER-Ptoid persistence when there are genotypes with the null AOP2\(-\) allele in the population.
It’s also important to note that there is a region of parameter space where LYER goes extinct before Ptoid, which would eventually lead to the collapse of the Ptoid since it has lost its resource. This is not possible if I were to assume the community is at equilibrium.
cbPalette <- viridis::viridis(4)
with(bind_rows(FE_LP_aop2_df, FE_LP_AOP2_df), table(species))species
LYER Ptoid
1681 1766 plot_fig_S9 <- bind_rows(FE_LP_aop2_df, FE_LP_AOP2_df) %>%
mutate(allele = factor(allele, labels = c("AOP2\u2013","AOP2+"))) %>%
mutate(species = ifelse(is.na(species) == T, "Food chain persists", species),
fspecies = factor(species, levels = c("LYER","Ptoid","Food chain persists"), labels = c("Arabidopsis only","Aphid only","Food chain persists"))) %>%
ggplot(., aes(x = LYER, y = Ptoid)) + # fspecies
geom_tile(aes(fill = fspecies)) +
#geom_point(data = LP_actual_df, aes(x = log_LYER_t_adj_first, y = log_Ptoid_t_adj_first)) +
facet_grid(~allele) +
#scale_fill_viridis_d(name = "Critical transition") +
scale_fill_manual(name = "Food-web transition", values = cbPalette[1:3]) +
coord_cartesian(xlim = c(1, max(LP_actual_df$log_LYER_t_first)),
ylim = c(1, max(LP_actual_df$log_Ptoid_t_first))) +
xlab("Aphid initial abundance (log scale)") +
ylab("Parasitoid initial abundance (log scale)") +
theme_cowplot() +
theme(strip.background = element_blank())
x11(); plot_fig_S9
#ggsave(plot = plot_fig_S9, filename = "figures/MAR1-nonequilibrium-foodchain-AOP2.pdf", width = 6, height = 5, device=cairo_pdf)
Alternative models
Drop BRBR -> LYER
# update LYER formula
all.LYER.bf.noBRBRonLYER <- update(all.LYER.bf, .~. -me(logBRBR_t, se_logBRBRt))
all.mar1.brm.unadj.noBRBRonLYER <- brm(
data = full_df,
formula = mvbf(all.BRBR.bf, all.LYER.bf.noBRBRonLYER, all.Ptoid.bf) + set_rescor(TRUE),
iter = 5000,
save_pars = save_pars(latent = TRUE, all = TRUE),
prior = all.mv.priors[-17,],
file = "output/all.mar1.brm.unadj.noBRBRonLYER.rds")
all.mar1.brm.unadj.noBRBRonLYER Family: MV(gaussian, gaussian, gaussian)
Links: mu = identity; sigma = identity
mu = identity; sigma = identity
mu = identity; sigma = identity
Formula: log(BRBR_t1) ~ 0 + Intercept + (me(logBRBR_t, se_logBRBRt) + me(logLYER_t, se_logLYERt) + me(logPtoid_t, se_logPtoidt)) + aop2_genotypes + AOP2_genotypes + temp + log(Biomass_g_t1) + (1 | p | Cage)
log(LYER_t1) ~ Intercept + me(logLYER_t, se_logLYERt) + me(logPtoid_t, se_logPtoidt) + aop2_genotypes + AOP2_genotypes + temp + log(Biomass_g_t1) + (1 | p | Cage) - 1
log(Ptoid_t1) ~ 0 + Intercept + me(logPtoid_t, se_logPtoidt) + me(logBRBR_t, se_logBRBRt) + me(logLYER_t, se_logLYERt) + aop2_genotypes + AOP2_genotypes + temp + log(Biomass_g_t1) + (1 | p | Cage)
Data: full_df (Number of observations: 264)
Draws: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
total post-warmup draws = 10000
Group-Level Effects:
~Cage (Number of levels: 60)
Estimate Est.Error l-95% CI
sd(logBRBRt1_Intercept) 0.19 0.13 0.01
sd(logLYERt1_Intercept) 0.09 0.07 0.00
sd(logPtoidt1_Intercept) 0.08 0.06 0.00
cor(logBRBRt1_Intercept,logLYERt1_Intercept) 0.04 0.49 -0.86
cor(logBRBRt1_Intercept,logPtoidt1_Intercept) 0.08 0.49 -0.85
cor(logLYERt1_Intercept,logPtoidt1_Intercept) 0.03 0.50 -0.87
u-95% CI Rhat Bulk_ESS Tail_ESS
sd(logBRBRt1_Intercept) 0.48 1.00 1726 3382
sd(logLYERt1_Intercept) 0.25 1.00 3642 4521
sd(logPtoidt1_Intercept) 0.21 1.00 5652 5130
cor(logBRBRt1_Intercept,logLYERt1_Intercept) 0.88 1.00 7182 6805
cor(logBRBRt1_Intercept,logPtoidt1_Intercept) 0.90 1.00 8649 6516
cor(logLYERt1_Intercept,logPtoidt1_Intercept) 0.88 1.00 7696 7899
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat
logBRBRt1_Intercept 2.33 0.59 1.18 3.48 1.00
logBRBRt1_aop2_genotypes 0.04 0.13 -0.22 0.30 1.00
logBRBRt1_AOP2_genotypes 0.01 0.12 -0.23 0.24 1.00
logBRBRt1_temp -0.52 0.08 -0.67 -0.37 1.00
logBRBRt1_logBiomass_g_t1 -0.02 0.25 -0.50 0.46 1.00
logLYERt1_Intercept 3.49 0.49 2.52 4.45 1.00
logLYERt1_aop2_genotypes 0.20 0.10 0.01 0.40 1.00
logLYERt1_AOP2_genotypes -0.08 0.09 -0.26 0.10 1.00
logLYERt1_temp -0.12 0.04 -0.21 -0.04 1.00
logLYERt1_logBiomass_g_t1 -0.07 0.21 -0.47 0.33 1.00
logPtoidt1_Intercept -1.94 0.53 -2.98 -0.90 1.00
logPtoidt1_aop2_genotypes 0.22 0.11 0.01 0.43 1.00
logPtoidt1_AOP2_genotypes -0.11 0.10 -0.30 0.08 1.00
logPtoidt1_temp -0.15 0.06 -0.27 -0.03 1.00
logPtoidt1_logBiomass_g_t1 -1.12 0.22 -1.54 -0.68 1.00
logBRBRt1_melogBRBR_tse_logBRBRt 0.60 0.09 0.43 0.77 1.00
logBRBRt1_melogLYER_tse_logLYERt 0.13 0.12 -0.09 0.36 1.00
logBRBRt1_melogPtoid_tse_logPtoidt -0.75 0.06 -0.87 -0.63 1.00
logLYERt1_melogLYER_tse_logLYERt 0.52 0.08 0.37 0.68 1.00
logLYERt1_melogPtoid_tse_logPtoidt -0.71 0.05 -0.81 -0.62 1.00
logPtoidt1_melogPtoid_tse_logPtoidt 0.99 0.05 0.89 1.09 1.00
logPtoidt1_melogBRBR_tse_logBRBRt 0.08 0.08 -0.07 0.22 1.00
logPtoidt1_melogLYER_tse_logLYERt 0.45 0.10 0.26 0.63 1.00
Bulk_ESS Tail_ESS
logBRBRt1_Intercept 4900 6370
logBRBRt1_aop2_genotypes 6588 7044
logBRBRt1_AOP2_genotypes 7221 7139
logBRBRt1_temp 4242 4897
logBRBRt1_logBiomass_g_t1 5461 7121
logLYERt1_Intercept 3949 5754
logLYERt1_aop2_genotypes 7248 7478
logLYERt1_AOP2_genotypes 7640 7223
logLYERt1_temp 10199 7636
logLYERt1_logBiomass_g_t1 5320 5978
logPtoidt1_Intercept 4843 6482
logPtoidt1_aop2_genotypes 8135 7125
logPtoidt1_AOP2_genotypes 9386 7823
logPtoidt1_temp 6076 6765
logPtoidt1_logBiomass_g_t1 6296 6744
logBRBRt1_melogBRBR_tse_logBRBRt 3616 4579
logBRBRt1_melogLYER_tse_logLYERt 3274 5075
logBRBRt1_melogPtoid_tse_logPtoidt 6194 7109
logLYERt1_melogLYER_tse_logLYERt 4190 5881
logLYERt1_melogPtoid_tse_logPtoidt 6360 6526
logPtoidt1_melogPtoid_tse_logPtoidt 7587 6935
logPtoidt1_melogBRBR_tse_logBRBRt 5280 6362
logPtoidt1_melogLYER_tse_logLYERt 4626 6567
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma_logBRBRt1 1.24 0.06 1.13 1.36 1.00 8048 7011
sigma_logLYERt1 0.98 0.05 0.90 1.07 1.00 9923 7548
sigma_logPtoidt1 1.02 0.05 0.94 1.12 1.00 11586 7535
Residual Correlations:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
rescor(logBRBRt1,logLYERt1) 0.49 0.05 0.38 0.58 1.00 7244
rescor(logBRBRt1,logPtoidt1) -0.22 0.06 -0.34 -0.10 1.00 9870
rescor(logLYERt1,logPtoidt1) -0.06 0.07 -0.19 0.07 1.00 10875
Tail_ESS
rescor(logBRBRt1,logLYERt1) 7667
rescor(logBRBRt1,logPtoidt1) 7561
rescor(logLYERt1,logPtoidt1) 7909
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).Carrying capacity check
# BRBR carrying capacity
fixef(all.mar1.brm.unadj.noBRBRonLYER)["logBRBRt1_Intercept","Estimate"] / (1 - fixef(all.mar1.brm.unadj.noBRBRonLYER)["logBRBRt1_melogBRBR_tse_logBRBRt","Estimate"]) # a bit low I think[1] 5.876553max(log(full_df$BRBR_t1))[1] 6.727432# LYER carrying capacity
fixef(all.mar1.brm.unadj.noBRBRonLYER)["logLYERt1_Intercept","Estimate"] / (1 - fixef(all.mar1.brm.unadj.noBRBRonLYER)["logLYERt1_melogLYER_tse_logLYERt","Estimate"])[1] 7.285828max(log(full_df$LYER_t1))[1] 6.975414Structural stability check
stability_all.mar1.brm.unadj.noBRBRonLYER <- aop2_vs_AOP2_posterior_samples_unadj(all.mar1.brm.unadj.noBRBRonLYER, n.geno = 2, temp.value = 0, logbiomass.value = 0)
stability_all.mar1.brm.unadj.noBRBRonLYER$aop2_SS_LP_BayesP # same clear inference[1] 0.9974AOP2 alters interaction matrix
## Update formula ----
# BRBR
all.BRBR.bf.xAOP2 <- update(all.BRBR.bf, .~. + (me(logBRBR_t, se_logBRBRt) + me(logLYER_t, se_logLYERt) + me(logPtoid_t, se_logPtoidt)):(aop2_genotypes + AOP2_genotypes))
# LYER
all.LYER.bf.xAOP2 <- update(all.LYER.bf, .~. + (me(logBRBR_t, se_logBRBRt) + me(logLYER_t, se_logLYERt) + me(logPtoid_t, se_logPtoidt)):(aop2_genotypes + AOP2_genotypes))
# Ptoid
all.Ptoid.bf.xAOP2 <- update(all.Ptoid.bf, .~. + (me(logBRBR_t, se_logBRBRt) + me(logLYER_t, se_logLYERt) + me(logPtoid_t, se_logPtoidt)):(aop2_genotypes + AOP2_genotypes))
## Update priors ----
all.mv.xAOP2.priors <- c(
# same priors as before
all.mv.priors,
# aop2 effects on interaction matrix
set_prior(prior.AOP2, class = "b", coef = "melogBRBR_tse_logBRBRt:aop2_genotypes", resp = "logBRBRt1"),
set_prior(prior.AOP2, class = "b", coef = "melogLYER_tse_logLYERt:aop2_genotypes", resp = "logLYERt1"),
set_prior(prior.AOP2, class = "b", coef = "melogPtoid_tse_logPtoidt:aop2_genotypes", resp = "logPtoidt1"),
set_prior(prior.AOP2, class = "b", coef = "melogLYER_tse_logLYERt:aop2_genotypes", resp = "logBRBRt1"),
set_prior(prior.AOP2, class = "b", coef = "melogBRBR_tse_logBRBRt:aop2_genotypes", resp = "logLYERt1"),
set_prior(prior.AOP2, class = "b", coef = "melogPtoid_tse_logPtoidt:aop2_genotypes", resp = "logBRBRt1"),
set_prior(prior.AOP2, class = "b", coef = "melogPtoid_tse_logPtoidt:aop2_genotypes", resp = "logLYERt1"),
set_prior(prior.AOP2, class = "b", coef = "melogBRBR_tse_logBRBRt:aop2_genotypes", resp = "logPtoidt1"),
set_prior(prior.AOP2, class = "b", coef = "melogLYER_tse_logLYERt:aop2_genotypes", resp = "logPtoidt1"),
# AOP2 effects on interaction matrix
set_prior(prior.AOP2, class = "b", coef = "melogBRBR_tse_logBRBRt:AOP2_genotypes", resp = "logBRBRt1"),
set_prior(prior.AOP2, class = "b", coef = "melogLYER_tse_logLYERt:AOP2_genotypes", resp = "logLYERt1"),
set_prior(prior.AOP2, class = "b", coef = "melogPtoid_tse_logPtoidt:AOP2_genotypes", resp = "logPtoidt1"),
set_prior(prior.AOP2, class = "b", coef = "melogLYER_tse_logLYERt:AOP2_genotypes", resp = "logBRBRt1"),
set_prior(prior.AOP2, class = "b", coef = "melogBRBR_tse_logBRBRt:AOP2_genotypes", resp = "logLYERt1"),
set_prior(prior.AOP2, class = "b", coef = "melogPtoid_tse_logPtoidt:AOP2_genotypes", resp = "logBRBRt1"),
set_prior(prior.AOP2, class = "b", coef = "melogPtoid_tse_logPtoidt:AOP2_genotypes", resp = "logLYERt1"),
set_prior(prior.AOP2, class = "b", coef = "melogBRBR_tse_logBRBRt:AOP2_genotypes", resp = "logPtoidt1"),
set_prior(prior.AOP2, class = "b", coef = "melogLYER_tse_logLYERt:AOP2_genotypes", resp = "logPtoidt1"))
## Fit model ----
all.mar1.brm.unadj.xAOP2 <- brm(
data = full_df,
formula = mvbf(all.BRBR.bf.xAOP2, all.LYER.bf.xAOP2, all.Ptoid.bf.xAOP2) + set_rescor(TRUE),
iter = 6000, # increased to avoid Bulk ESS warnings
save_pars = save_pars(latent = TRUE, all = TRUE),
prior = all.mv.xAOP2.priors,
file = "output/all.mar1.brm.unadj.xAOP2.rds")
all.mar1.brm.unadj.xAOP2 Family: MV(gaussian, gaussian, gaussian)
Links: mu = identity; sigma = identity
mu = identity; sigma = identity
mu = identity; sigma = identity
Formula: log(BRBR_t1) ~ Intercept + me(logBRBR_t, se_logBRBRt) + me(logLYER_t, se_logLYERt) + me(logPtoid_t, se_logPtoidt) + aop2_genotypes + AOP2_genotypes + temp + log(Biomass_g_t1) + (1 | p | Cage) + me(logBRBR_t, se_logBRBRt):aop2_genotypes + me(logBRBR_t, se_logBRBRt):AOP2_genotypes + me(logLYER_t, se_logLYERt):aop2_genotypes + me(logLYER_t, se_logLYERt):AOP2_genotypes + me(logPtoid_t, se_logPtoidt):aop2_genotypes + me(logPtoid_t, se_logPtoidt):AOP2_genotypes - 1
log(LYER_t1) ~ Intercept + me(logLYER_t, se_logLYERt) + me(logBRBR_t, se_logBRBRt) + me(logPtoid_t, se_logPtoidt) + aop2_genotypes + AOP2_genotypes + temp + log(Biomass_g_t1) + (1 | p | Cage) + me(logBRBR_t, se_logBRBRt):aop2_genotypes + me(logBRBR_t, se_logBRBRt):AOP2_genotypes + me(logLYER_t, se_logLYERt):aop2_genotypes + me(logLYER_t, se_logLYERt):AOP2_genotypes + me(logPtoid_t, se_logPtoidt):aop2_genotypes + me(logPtoid_t, se_logPtoidt):AOP2_genotypes - 1
log(Ptoid_t1) ~ Intercept + me(logPtoid_t, se_logPtoidt) + me(logBRBR_t, se_logBRBRt) + me(logLYER_t, se_logLYERt) + aop2_genotypes + AOP2_genotypes + temp + log(Biomass_g_t1) + (1 | p | Cage) + me(logBRBR_t, se_logBRBRt):aop2_genotypes + me(logBRBR_t, se_logBRBRt):AOP2_genotypes + me(logLYER_t, se_logLYERt):aop2_genotypes + me(logLYER_t, se_logLYERt):AOP2_genotypes + me(logPtoid_t, se_logPtoidt):aop2_genotypes + me(logPtoid_t, se_logPtoidt):AOP2_genotypes - 1
Data: full_df (Number of observations: 264)
Draws: 4 chains, each with iter = 6000; warmup = 3000; thin = 1;
total post-warmup draws = 12000
Group-Level Effects:
~Cage (Number of levels: 60)
Estimate Est.Error l-95% CI
sd(logBRBRt1_Intercept) 0.19 0.13 0.01
sd(logLYERt1_Intercept) 0.16 0.10 0.01
sd(logPtoidt1_Intercept) 0.07 0.06 0.00
cor(logBRBRt1_Intercept,logLYERt1_Intercept) -0.09 0.49 -0.90
cor(logBRBRt1_Intercept,logPtoidt1_Intercept) 0.10 0.51 -0.86
cor(logLYERt1_Intercept,logPtoidt1_Intercept) -0.03 0.50 -0.89
u-95% CI Rhat Bulk_ESS Tail_ESS
sd(logBRBRt1_Intercept) 0.46 1.00 2421 5049
sd(logLYERt1_Intercept) 0.38 1.00 2268 4082
sd(logPtoidt1_Intercept) 0.21 1.00 7963 5988
cor(logBRBRt1_Intercept,logLYERt1_Intercept) 0.84 1.00 3997 6271
cor(logBRBRt1_Intercept,logPtoidt1_Intercept) 0.91 1.00 11321 8313
cor(logLYERt1_Intercept,logPtoidt1_Intercept) 0.86 1.00 11121 9465
Population-Level Effects:
Estimate Est.Error l-95% CI
logBRBRt1_Intercept 2.35 0.67 1.02
logBRBRt1_aop2_genotypes -0.00 0.39 -0.78
logBRBRt1_AOP2_genotypes -0.04 0.38 -0.79
logBRBRt1_temp -0.44 0.08 -0.60
logBRBRt1_logBiomass_g_t1 -0.03 0.25 -0.53
logLYERt1_Intercept 2.85 0.62 1.64
logLYERt1_aop2_genotypes 0.09 0.35 -0.60
logLYERt1_AOP2_genotypes 0.53 0.35 -0.15
logLYERt1_temp 0.01 0.06 -0.11
logLYERt1_logBiomass_g_t1 -0.05 0.21 -0.45
logPtoidt1_Intercept -1.72 0.65 -3.00
logPtoidt1_aop2_genotypes 0.00 0.38 -0.73
logPtoidt1_AOP2_genotypes -0.02 0.37 -0.76
logPtoidt1_temp -0.15 0.06 -0.27
logPtoidt1_logBiomass_g_t1 -1.15 0.22 -1.58
logBRBRt1_melogBRBR_tse_logBRBRt 0.81 0.16 0.48
logBRBRt1_melogLYER_tse_logLYERt -0.04 0.17 -0.37
logBRBRt1_melogPtoid_tse_logPtoidt -0.81 0.12 -1.03
logBRBRt1_melogBRBR_tse_logBRBRt:aop2_genotypes -0.01 0.10 -0.20
logBRBRt1_melogBRBR_tse_logBRBRt:AOP2_genotypes -0.07 0.09 -0.25
logBRBRt1_melogLYER_tse_logLYERt:aop2_genotypes -0.03 0.10 -0.22
logBRBRt1_melogLYER_tse_logLYERt:AOP2_genotypes 0.08 0.10 -0.11
logBRBRt1_melogPtoid_tse_logPtoidt:aop2_genotypes 0.10 0.07 -0.04
logBRBRt1_melogPtoid_tse_logPtoidt:AOP2_genotypes -0.03 0.07 -0.17
logLYERt1_melogLYER_tse_logLYERt 0.12 0.15 -0.17
logLYERt1_melogBRBR_tse_logBRBRt 0.55 0.13 0.29
logLYERt1_melogPtoid_tse_logPtoidt -0.67 0.09 -0.85
logLYERt1_melogBRBR_tse_logBRBRt:aop2_genotypes -0.13 0.08 -0.28
logLYERt1_melogBRBR_tse_logBRBRt:AOP2_genotypes -0.15 0.07 -0.29
logLYERt1_melogLYER_tse_logLYERt:aop2_genotypes 0.10 0.08 -0.06
logLYERt1_melogLYER_tse_logLYERt:AOP2_genotypes 0.06 0.08 -0.10
logLYERt1_melogPtoid_tse_logPtoidt:aop2_genotypes 0.09 0.06 -0.03
logLYERt1_melogPtoid_tse_logPtoidt:AOP2_genotypes -0.11 0.06 -0.22
logPtoidt1_melogPtoid_tse_logPtoidt 0.98 0.10 0.78
logPtoidt1_melogBRBR_tse_logBRBRt 0.33 0.14 0.06
logPtoidt1_melogLYER_tse_logLYERt 0.19 0.15 -0.11
logPtoidt1_melogBRBR_tse_logBRBRt:aop2_genotypes -0.13 0.08 -0.28
logPtoidt1_melogBRBR_tse_logBRBRt:AOP2_genotypes -0.11 0.08 -0.26
logPtoidt1_melogLYER_tse_logLYERt:aop2_genotypes 0.13 0.09 -0.03
logPtoidt1_melogLYER_tse_logLYERt:AOP2_genotypes 0.10 0.08 -0.06
logPtoidt1_melogPtoid_tse_logPtoidt:aop2_genotypes 0.05 0.06 -0.07
logPtoidt1_melogPtoid_tse_logPtoidt:AOP2_genotypes -0.05 0.06 -0.16
u-95% CI Rhat Bulk_ESS
logBRBRt1_Intercept 3.65 1.00 8250
logBRBRt1_aop2_genotypes 0.76 1.00 9742
logBRBRt1_AOP2_genotypes 0.70 1.00 9109
logBRBRt1_temp -0.29 1.00 6432
logBRBRt1_logBiomass_g_t1 0.47 1.00 11376
logLYERt1_Intercept 4.07 1.00 7045
logLYERt1_aop2_genotypes 0.79 1.00 9122
logLYERt1_AOP2_genotypes 1.22 1.00 7953
logLYERt1_temp 0.13 1.00 6592
logLYERt1_logBiomass_g_t1 0.36 1.00 10690
logPtoidt1_Intercept -0.44 1.00 8446
logPtoidt1_aop2_genotypes 0.74 1.00 9598
logPtoidt1_AOP2_genotypes 0.71 1.00 8896
logPtoidt1_temp -0.03 1.00 10675
logPtoidt1_logBiomass_g_t1 -0.71 1.00 12565
logBRBRt1_melogBRBR_tse_logBRBRt 1.12 1.00 4383
logBRBRt1_melogLYER_tse_logLYERt 0.30 1.00 4727
logBRBRt1_melogPtoid_tse_logPtoidt -0.58 1.00 6268
logBRBRt1_melogBRBR_tse_logBRBRt:aop2_genotypes 0.18 1.00 5734
logBRBRt1_melogBRBR_tse_logBRBRt:AOP2_genotypes 0.12 1.00 6401
logBRBRt1_melogLYER_tse_logLYERt:aop2_genotypes 0.17 1.00 5784
logBRBRt1_melogLYER_tse_logLYERt:AOP2_genotypes 0.27 1.00 6072
logBRBRt1_melogPtoid_tse_logPtoidt:aop2_genotypes 0.24 1.00 7314
logBRBRt1_melogPtoid_tse_logPtoidt:AOP2_genotypes 0.10 1.00 8784
logLYERt1_melogLYER_tse_logLYERt 0.42 1.00 4436
logLYERt1_melogBRBR_tse_logBRBRt 0.80 1.00 4427
logLYERt1_melogPtoid_tse_logPtoidt -0.49 1.00 6651
logLYERt1_melogBRBR_tse_logBRBRt:aop2_genotypes 0.02 1.00 5821
logLYERt1_melogBRBR_tse_logBRBRt:AOP2_genotypes -0.00 1.00 6166
logLYERt1_melogLYER_tse_logLYERt:aop2_genotypes 0.26 1.00 6202
logLYERt1_melogLYER_tse_logLYERt:AOP2_genotypes 0.21 1.00 5591
logLYERt1_melogPtoid_tse_logPtoidt:aop2_genotypes 0.20 1.00 7300
logLYERt1_melogPtoid_tse_logPtoidt:AOP2_genotypes -0.00 1.00 9162
logPtoidt1_melogPtoid_tse_logPtoidt 1.17 1.00 7012
logPtoidt1_melogBRBR_tse_logBRBRt 0.60 1.00 4759
logPtoidt1_melogLYER_tse_logLYERt 0.48 1.00 4981
logPtoidt1_melogBRBR_tse_logBRBRt:aop2_genotypes 0.03 1.00 6649
logPtoidt1_melogBRBR_tse_logBRBRt:AOP2_genotypes 0.03 1.00 6903
logPtoidt1_melogLYER_tse_logLYERt:aop2_genotypes 0.29 1.00 6382
logPtoidt1_melogLYER_tse_logLYERt:AOP2_genotypes 0.27 1.00 6471
logPtoidt1_melogPtoid_tse_logPtoidt:aop2_genotypes 0.18 1.00 8181
logPtoidt1_melogPtoid_tse_logPtoidt:AOP2_genotypes 0.07 1.00 9301
Tail_ESS
logBRBRt1_Intercept 8641
logBRBRt1_aop2_genotypes 8766
logBRBRt1_AOP2_genotypes 8472
logBRBRt1_temp 7380
logBRBRt1_logBiomass_g_t1 9052
logLYERt1_Intercept 7962
logLYERt1_aop2_genotypes 9173
logLYERt1_AOP2_genotypes 8537
logLYERt1_temp 7007
logLYERt1_logBiomass_g_t1 9333
logPtoidt1_Intercept 8150
logPtoidt1_aop2_genotypes 8694
logPtoidt1_AOP2_genotypes 8905
logPtoidt1_temp 9495
logPtoidt1_logBiomass_g_t1 9380
logBRBRt1_melogBRBR_tse_logBRBRt 6291
logBRBRt1_melogLYER_tse_logLYERt 6595
logBRBRt1_melogPtoid_tse_logPtoidt 7867
logBRBRt1_melogBRBR_tse_logBRBRt:aop2_genotypes 7981
logBRBRt1_melogBRBR_tse_logBRBRt:AOP2_genotypes 7482
logBRBRt1_melogLYER_tse_logLYERt:aop2_genotypes 7720
logBRBRt1_melogLYER_tse_logLYERt:AOP2_genotypes 7428
logBRBRt1_melogPtoid_tse_logPtoidt:aop2_genotypes 8702
logBRBRt1_melogPtoid_tse_logPtoidt:AOP2_genotypes 9321
logLYERt1_melogLYER_tse_logLYERt 6560
logLYERt1_melogBRBR_tse_logBRBRt 6765
logLYERt1_melogPtoid_tse_logPtoidt 8133
logLYERt1_melogBRBR_tse_logBRBRt:aop2_genotypes 8428
logLYERt1_melogBRBR_tse_logBRBRt:AOP2_genotypes 7545
logLYERt1_melogLYER_tse_logLYERt:aop2_genotypes 7801
logLYERt1_melogLYER_tse_logLYERt:AOP2_genotypes 7431
logLYERt1_melogPtoid_tse_logPtoidt:aop2_genotypes 8099
logLYERt1_melogPtoid_tse_logPtoidt:AOP2_genotypes 8485
logPtoidt1_melogPtoid_tse_logPtoidt 7899
logPtoidt1_melogBRBR_tse_logBRBRt 7047
logPtoidt1_melogLYER_tse_logLYERt 7167
logPtoidt1_melogBRBR_tse_logBRBRt:aop2_genotypes 8171
logPtoidt1_melogBRBR_tse_logBRBRt:AOP2_genotypes 8276
logPtoidt1_melogLYER_tse_logLYERt:aop2_genotypes 7798
logPtoidt1_melogLYER_tse_logLYERt:AOP2_genotypes 8083
logPtoidt1_melogPtoid_tse_logPtoidt:aop2_genotypes 8776
logPtoidt1_melogPtoid_tse_logPtoidt:AOP2_genotypes 8767
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma_logBRBRt1 1.23 0.06 1.12 1.35 1.00 11865 9993
sigma_logLYERt1 0.93 0.05 0.85 1.02 1.00 8899 8868
sigma_logPtoidt1 1.02 0.05 0.93 1.11 1.00 14446 9705
Residual Correlations:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
rescor(logBRBRt1,logLYERt1) 0.48 0.05 0.37 0.57 1.00 9700
rescor(logBRBRt1,logPtoidt1) -0.24 0.06 -0.36 -0.11 1.00 12106
rescor(logLYERt1,logPtoidt1) -0.10 0.06 -0.22 0.03 1.00 14457
Tail_ESS
rescor(logBRBRt1,logLYERt1) 9690
rescor(logBRBRt1,logPtoidt1) 9228
rescor(logLYERt1,logPtoidt1) 9204
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).Carrying capacity check
B. brassicae:
# baseline
r_BRBR.xAOP2.all <- fixef(all.mar1.brm.unadj.xAOP2)["logBRBRt1_Intercept","Estimate"]
intra_BRBR.xAOP2.all <- fixef(all.mar1.brm.unadj.xAOP2)["logBRBRt1_melogBRBR_tse_logBRBRt","Estimate"]
K_BRBR.xAOP2.all_base <- r_BRBR.xAOP2.all / (1 - intra_BRBR.xAOP2.all) # too high
# aop2 effect
aop2_r_BRBR.xAOP2.all <- fixef(all.mar1.brm.unadj.xAOP2)["logBRBRt1_aop2_genotypes","Estimate"]
aop2_intra_BRBR.xAOP2.all <- fixef(all.mar1.brm.unadj.xAOP2)["logBRBRt1_melogBRBR_tse_logBRBRt:aop2_genotypes","Estimate"]
K_BRBR.xAOP2.all_aop2 <- (r_BRBR.xAOP2.all + aop2_r_BRBR.xAOP2.all) / (1 - (intra_BRBR.xAOP2.all + aop2_intra_BRBR.xAOP2.all))
# AOP2 effect
AOP2_r_BRBR.xAOP2.all <- fixef(all.mar1.brm.unadj.xAOP2)["logBRBRt1_AOP2_genotypes","Estimate"]
AOP2_intra_BRBR.xAOP2.all <- fixef(all.mar1.brm.unadj.xAOP2)["logBRBRt1_melogBRBR_tse_logBRBRt:AOP2_genotypes","Estimate"]
K_BRBR.xAOP2.all_AOP2 <- (r_BRBR.xAOP2.all + AOP2_r_BRBR.xAOP2.all) / (1 - (intra_BRBR.xAOP2.all + AOP2_intra_BRBR.xAOP2.all))
# compare aop2 vs baseline carrying capacity
K_BRBR.xAOP2.all_base - K_BRBR.xAOP2.all_aop2 > 0 # base > aop2? Says it is, which doesn't make sense. [1] TRUE# compare AOP2 vs baseline
K_BRBR.xAOP2.all_base - K_BRBR.xAOP2.all_AOP2 > 0 # base > AOP2? Says it is, which is expected[1] TRUEL. erysimi:
# baseline
r_LYER.xAOP2.all <- fixef(all.mar1.brm.unadj.xAOP2)["logLYERt1_Intercept","Estimate"]
intra_LYER.xAOP2.all <- fixef(all.mar1.brm.unadj.xAOP2)["logLYERt1_melogLYER_tse_logLYERt","Estimate"]
K_LYER.xAOP2.all_base <- r_LYER.xAOP2.all / (1 - intra_LYER.xAOP2.all) # too low
# aop2 effect
aop2_r_LYER.xAOP2.all <- fixef(all.mar1.brm.unadj.xAOP2)["logLYERt1_aop2_genotypes","Estimate"]
aop2_intra_LYER.xAOP2.all <- fixef(all.mar1.brm.unadj.xAOP2)["logLYERt1_melogLYER_tse_logLYERt:aop2_genotypes","Estimate"]
K_LYER.xAOP2.all_aop2 <- (r_LYER.xAOP2.all + aop2_r_LYER.xAOP2.all) / (1 - (intra_LYER.xAOP2.all + aop2_intra_LYER.xAOP2.all))
# AOP2 effect
AOP2_r_LYER.xAOP2.all <- fixef(all.mar1.brm.unadj.xAOP2)["logLYERt1_AOP2_genotypes","Estimate"]
AOP2_intra_LYER.xAOP2.all <- fixef(all.mar1.brm.unadj.xAOP2)["logLYERt1_melogLYER_tse_logLYERt:AOP2_genotypes","Estimate"]
K_LYER.xAOP2.all_AOP2 <- (r_LYER.xAOP2.all + AOP2_r_LYER.xAOP2.all) / (1 - (intra_LYER.xAOP2.all + AOP2_intra_LYER.xAOP2.all))
# compare aop2 vs baseline carrying capacity
K_LYER.xAOP2.all_base - K_LYER.xAOP2.all_aop2 > 0 # base > aop2? Says it isn't, which makes sense. [1] FALSE# compare AOP2 vs baseline
K_LYER.xAOP2.all_base - K_LYER.xAOP2.all_AOP2 > 0 # base > AOP2? Says it isn't, but I would expect it to be[1] FALSEThis model gives unrealistic estimates for B. brassicae’s carrying capacity (ridiculously high, 176410 individuals) and L. erysimi’s carrying capacity (way too low, 25 individuals). Moreover, this model predicts that adding genotypes with a null AOP2\(-\) allele to the plant population will actually decrease the aphids carrying capacity, despite the documented positive effect on plant growth. All of these biological predictions are unreasonable, so we don’t consider this model further.
MAR(2) model
## update formula
all.BRBR.bf.ar2 <- update(all.BRBR.bf, .~. + me(logBRBR_t0, se_logBRBRt) + me(logLYER_t0, se_logLYERt) + me(logPtoid_t0, se_logPtoidt))
all.LYER.bf.ar2 <- update(all.LYER.bf, .~. + me(logBRBR_t0, se_logBRBRt) + me(logLYER_t0, se_logLYERt) + me(logPtoid_t0, se_logPtoidt))
all.Ptoid.bf.ar2 <- update(all.Ptoid.bf, .~. + me(logBRBR_t0, se_logBRBRt) + me(logLYER_t0, se_logLYERt) + me(logPtoid_t0, se_logPtoidt))
## fit model
all.mar1.brm.unadj.ar2 <- brm(
data = full_df,
formula = mvbf(all.BRBR.bf.ar2, all.LYER.bf.ar2, all.Ptoid.bf.ar2) + set_rescor(TRUE),
iter = 5000,
save_pars = save_pars(latent = TRUE, all = TRUE),
prior = all.mv.priors,
file = "output/all.mar1.brm.unadj.ar2.lag.rds")Carrying capacity check
# BRBR carrying capacity
fixef(all.mar1.brm.unadj.ar2)["logBRBRt1_Intercept","Estimate"] / (1 - fixef(all.mar1.brm.unadj.ar2)["logBRBRt1_melogBRBR_tse_logBRBRt","Estimate"]) # way too high[1] 12.56733# LYER carrying capacity
fixef(all.mar1.brm.unadj)["logLYERt1_Intercept","Estimate"] / (1 - fixef(all.mar1.brm.unadj)["logLYERt1_melogLYER_tse_logLYERt","Estimate"]) # way too low[1] 4.675535Structural stability check
stability_all.mar1.brm.unadj.ar2 <- aop2_vs_AOP2_posterior_samples_unadj(all.mar1.brm.unadj.ar2, n.geno = 2, temp.value = 0, logbiomass.value = 0)
stability_all.mar1.brm.unadj.ar2$aop2_SS_LP_BayesP # same inference![1] 0.9942Reproduce Bayesian R2 in Table S4 Bayesian R2
## MAR(1) model
bayes_R2(all.mar1.brm.unadj, newdata = BRBR_predict_df, resp = "logBRBRt1") Estimate Est.Error Q2.5 Q97.5
R2logBRBRt1 0.6458209 0.02397003 0.5954573 0.6884681bayes_R2(all.mar1.brm.unadj, newdata = LYER_predict_df, resp = "logLYERt1") Estimate Est.Error Q2.5 Q97.5
R2logLYERt1 0.4869248 0.03251046 0.4239504 0.5506908bayes_R2(all.mar1.brm.unadj, newdata = Ptoid_predict_df, resp = "logPtoidt1") Estimate Est.Error Q2.5 Q97.5
R2logPtoidt1 0.6301973 0.01173414 0.6023809 0.6473038## MAR(2) model
bayes_R2(all.mar1.brm.unadj.ar2, newdata = BRBR_predict_df, resp = "logBRBRt1") Estimate Est.Error Q2.5 Q97.5
R2logBRBRt1 0.7281811 0.01366882 0.6990006 0.7525384bayes_R2(all.mar1.brm.unadj.ar2, newdata = LYER_predict_df, resp = "logLYERt1") Estimate Est.Error Q2.5 Q97.5
R2logLYERt1 0.5115748 0.03210331 0.447848 0.5736911bayes_R2(all.mar1.brm.unadj.ar2, newdata = Ptoid_predict_df, resp = "logPtoidt1") Estimate Est.Error Q2.5 Q97.5
R2logPtoidt1 0.5987022 0.01475788 0.5667047 0.6245675Adjusted data
## update formula ----
# BRBR
all.BRBR.bf.adj <- bf(log(BRBR_t1) ~ 0 + Intercept + (me(logBRBR_t, se_logBRBRt) + me(logLYER_t_adj, se_logLYERt) + me(logPtoid_t_adj, se_logPtoidt)) + aop2_genotypes + AOP2_genotypes + temp + log(Biomass_g_t1) + (1|p|Cage))
# LYER
all.LYER.bf.adj <- bf(log(LYER_t1_adj) ~ 0 + Intercept + (me(logLYER_t_adj, se_logLYERt) + me(logBRBR_t, se_logBRBRt) + me(logPtoid_t_adj, se_logPtoidt)) + aop2_genotypes + AOP2_genotypes + temp + log(Biomass_g_t1) + (1|p|Cage))
# Ptoid
all.Ptoid.bf.adj <- bf(log(Ptoid_t1_adj) ~ 0 + Intercept + me(logPtoid_t_adj, se_logPtoidt) + me(logBRBR_t, se_logBRBRt) + me(logLYER_t_adj, se_logLYERt) + aop2_genotypes + temp + AOP2_genotypes + log(Biomass_g_t1) + (1|p|Cage))
## update priors ----
all.mv.priors.adj <- c(
# aop2 and AOP2 effects
set_prior(prior.AOP2, class = "b", coef = "aop2_genotypes", resp = "logBRBRt1"),
set_prior(prior.AOP2, class = "b", coef = "AOP2_genotypes", resp = "logBRBRt1"),
set_prior(prior.AOP2, class = "b", coef = "aop2_genotypes", resp = "logLYERt1adj"),
set_prior(prior.AOP2, class = "b", coef = "AOP2_genotypes", resp = "logLYERt1adj"),
set_prior(prior.AOP2, class = "b", coef = "aop2_genotypes", resp = "logPtoidt1adj"),
set_prior(prior.AOP2, class = "b", coef = "AOP2_genotypes", resp = "logPtoidt1adj"),
# biomass effects
set_prior(prior.AphidBiomass, class = "b", coef = "logBiomass_g_t1", resp = "logBRBRt1"),
set_prior(prior.AphidBiomass, class = "b", coef = "logBiomass_g_t1", resp = "logLYERt1adj"),
set_prior(prior.PtoidBiomass, class = "b", coef = "logBiomass_g_t1", resp = "logPtoidt1adj"),
# temp effects
set_prior(prior.temp, class = "b", coef = "temp", resp = "logBRBRt1"),
set_prior(prior.temp, class = "b", coef = "temp", resp = "logLYERt1adj"),
set_prior(prior.temp, class = "b", coef = "temp", resp = "logPtoidt1adj"),
# baseline growth rates
set_prior(prior.r.BRBR, class = "b", coef = "Intercept", resp = "logBRBRt1"),
set_prior(prior.r.LYER, class = "b", coef = "Intercept", resp = "logLYERt1adj"),
set_prior(prior.r.Ptoid, class = "b", coef = "Intercept", resp = "logPtoidt1adj"),
# intraspecific effects
set_prior(prior.intra.BRBR, class = "b", coef = "melogBRBR_tse_logBRBRt", resp = "logBRBRt1"),
set_prior(prior.intra.LYER, class = "b", coef = "melogLYER_t_adjse_logLYERt", resp = "logLYERt1adj"),
set_prior(prior.intra.Ptoid, class = "b", coef = "melogPtoid_t_adjse_logPtoidt", resp = "logPtoidt1adj"),
# negative interspecific effects
set_prior(prior.LYERonBRBR, class = "b", coef = "melogLYER_t_adjse_logLYERt", resp = "logBRBRt1"),
set_prior(prior.BRBRonLYER, class = "b", coef = "melogBRBR_tse_logBRBRt", resp = "logLYERt1adj"),
set_prior(prior.PtoidonBRBR, class = "b", coef = "melogPtoid_t_adjse_logPtoidt", resp = "logBRBRt1"),
set_prior(prior.PtoidonLYER, class = "b", coef = "melogPtoid_t_adjse_logPtoidt", resp = "logLYERt1adj"),
# positive interspecific effects
set_prior(prior.BRBRonPtoid, class = "b", coef = "melogBRBR_tse_logBRBRt", resp = "logPtoidt1adj"),
set_prior(prior.LYERonPtoid, class = "b", coef = "melogLYER_t_adjse_logLYERt", resp = "logPtoidt1adj"),
# random effects
set_prior(prior.random.effects, class = "sd", resp = "logBRBRt1"),
set_prior(prior.random.effects, class = "sd", resp = "logLYERt1adj"),
set_prior(prior.random.effects, class = "sd", resp = "logPtoidt1adj"))
## fit model
all.mar1.brm.adj <- brm(
data = full_df,
formula = mvbf(all.BRBR.bf.adj, all.LYER.bf.adj, all.Ptoid.bf.adj) + set_rescor(TRUE),
iter = 5000,
save_pars = save_pars(latent = TRUE, all = TRUE),
prior = all.mv.priors.adj,
file = "output/all.mar1.brm.adj.rds")Structural stability check
stability_all.mar1.brm.adj <- aop2_vs_AOP2_posterior_samples_adj(all.mar1.brm.adj, n.geno = 2, temp.value = 0, logbiomass.value = 0)
stability_all.mar1.brm.adj$aop2_SS_LP_BayesP # same inference[1] 0.979Aphid intrinsic growth rates
Confirmation with independent data (not used in model for all species) that AOP\(-\) increases the intrinsic growth rate of the aphids.
Unadjusted data
# BRBR
initial.BRBR.bf <- bf(log(BRBR_t1) ~ 0 + Intercept + offset(log(4)) + aop2_genotypes + AOP2_genotypes + temp)
# LYER
initial.LYER.bf <- bf(log(LYER_t1) ~ 0 + Intercept + offset(log(4)) + aop2_genotypes + AOP2_genotypes + temp)
# Prior
initial.mv.prior <- c(
# aop2 effects
set_prior(prior.AOP2, class = "b", coef = "aop2_genotypes", resp = "logBRBRt1"),
set_prior(prior.AOP2, class = "b", coef = "AOP2_genotypes", resp = "logBRBRt1"),
set_prior(prior.AOP2, class = "b", coef = "aop2_genotypes", resp = "logLYERt1"),
set_prior(prior.AOP2, class = "b", coef = "AOP2_genotypes", resp = "logLYERt1"),
# temp effects
set_prior(prior.temp, class = "b", coef = "temp", resp = "logBRBRt1"),
set_prior(prior.temp, class = "b", coef = "temp", resp = "logLYERt1"),
# growth rates
set_prior(prior.r.BRBR, class = "b", coef = "Intercept", resp = "logBRBRt1"),
set_prior(prior.r.LYER, class = "b", coef = "Intercept", resp = "logLYERt1"))
# fit model
initial.mar1.brm.unadj <- brm(
data = filter(aphids_only_df, Week == 0),
formula = mvbf(initial.BRBR.bf, initial.LYER.bf) + set_rescor(TRUE),
iter = 5000,
save_pars = save_pars(latent = TRUE, all = TRUE),
prior = initial.mv.prior,
file = "output/initial.mar1.brm.unadj.rds")
# get posteriors
ps_initial.mar1.brm.unadj <- posterior_samples(initial.mar1.brm.unadj, pars = c("b_logBRBRt1_aop2_genotypes","b_logBRBRt1_AOP2_genotypes","b_logLYERt1_aop2_genotypes","b_logLYERt1_AOP2_genotypes")) %>%
mutate(aop2_vs_AOP2_r_BRBR = b_logBRBRt1_aop2_genotypes - b_logBRBRt1_AOP2_genotypes,
aop2_vs_AOP2_r_LYER = b_logLYERt1_aop2_genotypes - b_logLYERt1_AOP2_genotypes)
# r BRBR, aop2 vs AOP2
median(ps_initial.mar1.brm.unadj$aop2_vs_AOP2_r_BRBR)[1] 0.3298602quantile(ps_initial.mar1.brm.unadj$aop2_vs_AOP2_r_BRBR, probs = c(0.025,0.975)) 2.5% 97.5%
0.06062928 0.58563572 # r LYER, aop2 vs AOP2
median(ps_initial.mar1.brm.unadj$aop2_vs_AOP2_r_LYER)[1] 1.014212quantile(ps_initial.mar1.brm.unadj$aop2_vs_AOP2_r_LYER, probs = c(0.025,0.975)) 2.5% 97.5%
0.5872244 1.4224264 Adjusted data
# LYER
initial.LYER.bf.adj <- bf(log(LYER_t1_adj) ~ 0 + Intercept + offset(log(4)) + aop2_genotypes + AOP2_genotypes + temp)
# Priors
initial.mv.prior.adj <- c(
# aop2 effects
set_prior(prior.AOP2, class = "b", coef = "aop2_genotypes", resp = "logBRBRt1"),
set_prior(prior.AOP2, class = "b", coef = "AOP2_genotypes", resp = "logBRBRt1"),
set_prior(prior.AOP2, class = "b", coef = "aop2_genotypes", resp = "logLYERt1adj"),
set_prior(prior.AOP2, class = "b", coef = "AOP2_genotypes", resp = "logLYERt1adj"),
# temp effects
set_prior(prior.temp, class = "b", coef = "temp", resp = "logBRBRt1"),
set_prior(prior.temp, class = "b", coef = "temp", resp = "logLYERt1adj"),
# growth rates
set_prior(prior.r.BRBR, class = "b", coef = "Intercept", resp = "logBRBRt1"),
set_prior(prior.r.LYER, class = "b", coef = "Intercept", resp = "logLYERt1adj"))
# fit model
initial.mar1.brm.adj <- brm(
data = filter(aphids_only_df, Week == 0),
formula = mvbf(initial.BRBR.bf, initial.LYER.bf.adj) + set_rescor(TRUE),
iter = 5000,
save_pars = save_pars(latent = TRUE, all = TRUE),
prior = initial.mv.prior.adj,
file = "output/initial.mar1.brm.adj.rds")
ps_initial.mar1.brm.adj <- posterior_samples(initial.mar1.brm.adj, pars = c("b_logBRBRt1_aop2_genotypes","b_logBRBRt1_AOP2_genotypes","b_logLYERt1adj_aop2_genotypes","b_logLYERt1adj_AOP2_genotypes")) %>%
mutate(aop2_vs_AOP2_r_BRBR = b_logBRBRt1_aop2_genotypes - b_logBRBRt1_AOP2_genotypes,
aop2_vs_AOP2_r_LYER = b_logLYERt1adj_aop2_genotypes - b_logLYERt1adj_AOP2_genotypes)
# r BRBR, aop2 vs AOP2
median(ps_initial.mar1.brm.adj$aop2_vs_AOP2_r_BRBR)[1] 0.3369524quantile(ps_initial.mar1.brm.adj$aop2_vs_AOP2_r_BRBR, probs = c(0.025,0.975)) 2.5% 97.5%
0.0749945 0.5938041 # r LYER, aop2 vs AOP2
median(ps_initial.mar1.brm.adj$aop2_vs_AOP2_r_LYER)[1] 0.4574658quantile(ps_initial.mar1.brm.adj$aop2_vs_AOP2_r_LYER, probs = c(0.025,0.975)) 2.5% 97.5%
0.2832546 0.6327860
sessionInfo()R version 4.1.2 (2021-11-01)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 16.04.7 LTS
Matrix products: default
BLAS: /usr/lib/libblas/libblas.so.3.6.0
LAPACK: /usr/lib/lapack/liblapack.so.3.6.0
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] bayesplot_1.8.1 rgl_0.106.8 knitr_1.37 tidybayes_2.3.1
[5] matlib_0.9.4 cowplot_1.1.1 forcats_0.5.1 stringr_1.4.0
[9] dplyr_1.0.7 purrr_0.3.4 readr_2.1.1 tidyr_1.1.4
[13] tibble_3.1.6 ggplot2_3.3.5 tidyverse_1.3.1 brms_2.16.3
[17] Rcpp_1.0.7 RCurl_1.98-1.3 MASS_7.3-54 workflowr_1.6.2
loaded via a namespace (and not attached):
[1] readxl_1.3.1 backports_1.4.1 plyr_1.8.6
[4] igraph_1.2.11 svUnit_1.0.6 crosstalk_1.2.0
[7] rstantools_2.1.1 inline_0.3.19 digest_0.6.29
[10] htmltools_0.5.2 viridis_0.6.1 rsconnect_0.8.25
[13] fansi_1.0.0 magrittr_2.0.1 checkmate_2.0.0
[16] tzdb_0.2.0 openxlsx_4.2.4 modelr_0.1.8
[19] RcppParallel_5.1.5 matrixStats_0.61.0 xts_0.12.1
[22] prettyunits_1.1.1 colorspace_2.0-2 rvest_1.0.2
[25] ggdist_2.4.1 haven_2.4.3 xfun_0.29
[28] callr_3.7.0 crayon_1.4.2 jsonlite_1.7.2
[31] zoo_1.8-9 glue_1.6.0 gtable_0.3.0
[34] distributional_0.2.2 car_3.0-10 pkgbuild_1.3.1
[37] rstan_2.21.3 abind_1.4-5 scales_1.1.1
[40] mvtnorm_1.1-3 DBI_1.1.2 miniUI_0.1.1.1
[43] viridisLite_0.4.0 xtable_1.8-3 foreign_0.8-81
[46] stats4_4.1.2 StanHeaders_2.21.0-7 DT_0.20
[49] htmlwidgets_1.5.4 httr_1.4.2 threejs_0.3.3
[52] arrayhelpers_1.1-0 posterior_1.2.0 ellipsis_0.3.2
[55] pkgconfig_2.0.3 loo_2.4.1 farver_2.1.0
[58] sass_0.4.0 dbplyr_2.1.1 utf8_1.2.2
[61] labeling_0.4.2 tidyselect_1.1.1 rlang_0.4.12
[64] manipulateWidget_0.11.0 reshape2_1.4.4 later_1.3.0
[67] munsell_0.5.0 cellranger_1.1.0 tools_4.1.2
[70] cli_3.1.0 generics_0.1.1 broom_0.7.11
[73] ggridges_0.5.3 evaluate_0.14 fastmap_1.1.0
[76] yaml_2.2.1 processx_3.5.2 fs_1.5.2
[79] zip_2.2.0 nlme_3.1-152 whisker_0.4
[82] mime_0.12 xml2_1.3.3 compiler_4.1.2
[85] shinythemes_1.2.0 rstudioapi_0.13 curl_4.3.2
[88] reprex_2.0.1 bslib_0.3.1 stringi_1.7.3
[91] highr_0.9 ps_1.6.0 Brobdingnag_1.2-6
[94] lattice_0.20-45 Matrix_1.4-0 markdown_1.1
[97] shinyjs_2.1.0 tensorA_0.36.2 vctrs_0.3.8
[100] pillar_1.6.4 lifecycle_1.0.1 jquerylib_0.1.4
[103] bridgesampling_1.1-2 data.table_1.14.2 bitops_1.0-7
[106] httpuv_1.6.5 R6_2.5.1 promises_1.2.0.1
[109] gridExtra_2.3 rio_0.5.26 codetools_0.2-18
[112] colourpicker_1.1.1 gtools_3.9.2 assertthat_0.2.1
[115] rprojroot_2.0.2 withr_2.4.3 shinystan_2.5.0
[118] parallel_4.1.2 hms_1.1.1 grid_4.1.2
[121] coda_0.19-4 rmarkdown_2.11 carData_3.0-4
[124] git2r_0.28.0 shiny_1.7.1 lubridate_1.8.0
[127] base64enc_0.1-3 dygraphs_1.1.1.6