Codes for: Disturbance regime changes leave long-lasting legacies on a microbial community’s composition and function.
H. Inamine, L. Lear, A. Miller, S. Roxburgh, A. Buckling, K. Shea
Updated: June 05, 2026
1 Preamble
The codes below import packages, define functions, and format the data.
R.version$version.string[1] "R version 4.5.2 (2025-10-31)"R version 4.5.2 (2025-10-31)# The codes in this section define functions and format the data.
rm(list=ls())
# Import libraries.
library(tidyverse)
library(scales)
library(vegan)
library(gridExtra)
library(ggfortify)
library(lemon)
library(ggpubr)
library(car)
library(pscl)
library(emmeans)
# Set theme for ggplot.
theme_set(theme_bw(base_size=12))
# Changing the default behavior of ggsave() so it creates a directory.
ggsave <- function(...) ggplot2::ggsave(..., create.dir = TRUE)
# This inserts a ggplot layer under existing ones.
`-.gg` <- function(plot, layer) {
if (missing(layer)) {
stop("Cannot use `-.gg()` with a single argument. Did you accidentally put - on a new line?")
}
if (!is.ggplot(plot)) {
stop('Need a plot on the left side')
}
plot$layers = c(layer, plot$layers)
plot
}
# This formats the scientific notation in a plot.
scientific_10 <- function(x) {
ifelse(x==0, "0",
parse(text=gsub("e\\+*", " %*% 10^", scales::scientific_format()(x)))
)
}
# These are labels for facet_wrap()
freq1.labels = c(
`0.125` = "Past: every 8 d",
`0.25` = "every 4 d",
`1` = "every 1 d")
freq2.labels = c(
`0.125` = "Novel: every 8 d",
`0.25` = "every 4 d",
`1` = "every 1 d")
num.total.labels = c(
`3` = "# disturbed: 3",
`9` = "9",
`10` = "10")
# This calculates the community dissimilarity index
ens.beta = function(data){
# Calculate the alpha diversity.
alpha = diversity(data, index="shannon") %>%
mean() %>% exp()
# Calculate the gamma diversity.
gamma = colMeans(data/rowSums(data)) %>%
diversity(., index="shannon") %>% exp()
# Calculate the beta diversity.
beta = gamma/alpha
# Calculate the community dissimilarity index.
turnover = (beta - 1) / (nrow(data) - 1)
return(turnover)
}
# Read and format the data from the samples before the regime changes.
data.past = read.csv("./data_preinvasion_pub.csv", comment.char="#") %>%
tibble() %>%
mutate(
# Divide cfu/plate by below to get cfu/ml
total.dilution = 10^(plate.dilution.log10) * glycerol.dilution * plated.volume.ul / 1000,
# Number of passages over the first 8 days.
num1 = 8/freq1,
# Calculate the frequency of passages, i.e., unit of [1/time]
freq1 = 1/freq1,
# Calculate the 1st order effective number of species (alpha diversity)
alpha1 = exp(diversity(select(., sm.raw, ws.raw, fs.raw), index="shannon"))
)
# Read and format the data from the samples after the regime changes.
data.novel = read.csv("./data_postinvasion_pub.csv", comment.char="#") %>%
tibble() %>%
mutate(
# Sample codes to match with the data.past samples.
preinv.sample.code = paste(freq1, replicate, sep=""),
# Divide cfu/plate by below to get cfu/ml
total.dilution = 10^(plate.dilution.log10) * glycerol.dilution * plated.volume.ul / 1000,
# Number of passages over the first and second 8 days, respectively.
num1 = 8/freq1,
num2 = 8/freq2,
num.total = num1 + num2, # total number of passages
# Calculate the frequency of passages, i.e., unit of [1/time]
freq1 = 1/freq1,
freq2 = 1/freq2,
# Did the frequency increase over time?
inc.freq = freq2 > freq1,
# Calculate the 1st order effective number of species (alpha diversity)
alpha1 = exp(diversity(select(., sm.raw, ws.raw, fs.raw), index="shannon")),
# Placeholder for the community dissimilarity index (beta diversity)
beta.past = NA # beta compared to the past regime.
)
# Calculate the community dissimilarity index between the samples before and after the regime changes.
for(i in 1:nrow(data.novel)){
# Match the pre- and post- regime change samples.
ind = which(data.past$sample.code == data.novel[i, ]$preinv.sample.code)
temp = rbind(
data.novel[i, c("sm.raw", "ws.raw", "fs.raw")],
data.past[ind, c("sm.raw", "ws.raw", "fs.raw")]) %>%
as.data.frame()
# Calculate the index.
data.novel[i, "beta.past"] = ens.beta(temp)
}
# New data used for model predictions.
freq = seq(1/8, 1, length.out=50) %>%
c(0.125, 0.25, 1) %>%
unique() %>%
sort()
newdat1 = data.frame(
freq1 = freq)
newdat2 = expand.grid(
freq1 = freq,
freq2 = freq,
invaded = c(F, T),
total.dilution = 1,
bac = 1)
newdat3 = expand.grid(
freq1 = freq,
freq2 = freq,
invaded = c(F, T),
alpha1 = seq(1, 2, by=0.01),
total.dilution = 1,
bac = 1)2 Effects of disturbance regimes on diversity
2.1 Diversity before the regime change
Here, we analyze how diversity changes with disturbance frequency before the regime change.
Hypothesis: Following previous studies, diversity should show a unimodal pattern with increasing disturbance frequency.
Results: Disturbance-diversity relationship before the regime
change has a unimodal relationship:
alpha1 ~ (freq1 + I(freq1^2))
We have positive freq1 and negative I(freq1^2)
coefficients, and no covariates should be dropped.
# Fit a model with unimodal pattern.
lm.past.alpha1 = lm(alpha1 ~ (freq1 + I(freq1^2)), data = data.past)
summary(lm.past.alpha1)
Call:
lm(formula = alpha1 ~ (freq1 + I(freq1^2)), data = data.past)
Residuals:
Min 1Q Median 3Q Max
-0.31964 -0.03517 -0.02359 0.12258 0.30576
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.9692 0.1854 5.227 0.000102 ***
freq1 2.7554 1.1617 2.372 0.031509 *
I(freq1^2) -2.6895 0.9901 -2.716 0.015923 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.1716 on 15 degrees of freedom
Multiple R-squared: 0.5843, Adjusted R-squared: 0.5289
F-statistic: 10.54 on 2 and 15 DF, p-value: 0.001383# Nothing should be dropped.
drop1(lm.past.alpha1, test="F")# Generate the model predictions.
pred.past.alpha1 = predict(lm.past.alpha1, newdata = newdat1, se.fit = T) %>%
magrittr::extract(c("fit", "se.fit")) %>%
as.data.frame() %>%
cbind(newdat1, .) %>%
transform(lwr = fit - se.fit, upr = fit + se.fit)
# Store the plot.
fig.past.alpha1 = ggplot(data = data.past, mapping = aes(x = freq1, y = alpha1)) +
scale_x_continuous(trans = "log2",
breaks=1/c(1,4,8), labels=c("1\nHigh", "4", "8\nLow")) +
geom_jitter(aes(color = as.factor(freq1)), width=0.02, alpha=0.7, size=2.5, shape = 1, stroke = 1.5) +
xlab("Past freq. (every x days)") +
ylab("Effective number of species") +
scale_color_brewer(palette = "Dark2", guide = NULL) -
geom_line(data = pred.past.alpha1, aes(y = fit), show.legend=F, size=0.75) +
geom_ribbon(data = pred.past.alpha1, aes(y = NULL, ymin = lwr, ymax=upr), show.legend=F, color=NA, alpha=0.1)2.2 Diversity after the regime change
Here, we analyze how diversity changes with disturbance frequency after the regime change.
Hypothesis: Similar to the results above, disturbance-diversity relationship after the novel regimes would be unimodal. Furthermore, the effects of the past regimes would have disappeared under the novel regimes.
Results: The best model is:
alpha1 ~ freq1 + freq2 + I(freq2^2) + freq1:freq2.
We have positive freq2, negative I(freq2^2),
and negative freq1:freq2. freq1 is
non-significant.
Step-wise model selection on the full model shows that:
- Invasions have no significant effect on the resident diversity.
- Past regimes affect diversity under the novel regimes, but only
through the interaction term
freq2:freq1.
add1on a model without anyfreq1term confirms thatfreq1should be included in the model.ANOVA on a model where the covariates are treated as factors results in qualitatively similar result:
freq1is important, butinvasionis not.
# Fit models with differing complexities.
lm.novel.alpha1.full = lm(alpha1 ~ (freq1 + I(freq1^2)) * (freq2 + I(freq2^2)) * invaded,
data = data.novel)
lm.novel.alpha1.v1 = update(lm.novel.alpha1.full,
~ freq1 * (freq2 + I(freq2^2)) * invaded)
lm.novel.alpha1.v2 = update(lm.novel.alpha1.full,
~ (freq1 + I(freq1^2)) * freq2 * invaded)
# The best model has a unimodal relationship with respect to the novel regimes.
AIC(lm.novel.alpha1.full, lm.novel.alpha1.v1, lm.novel.alpha1.v2)summary(lm.novel.alpha1.v1)
Call:
lm(formula = alpha1 ~ freq1 + freq2 + I(freq2^2) + invaded +
freq1:freq2 + freq1:I(freq2^2) + freq1:invaded + freq2:invaded +
I(freq2^2):invaded + freq1:freq2:invaded + freq1:I(freq2^2):invaded,
data = data.novel)
Residuals:
Min 1Q Median 3Q Max
-0.58971 -0.12573 0.01205 0.11493 0.64400
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.695781 0.240411 2.894 0.004706 **
freq1 -0.008093 0.401033 -0.020 0.983941
freq2 5.786348 1.506198 3.842 0.000219 ***
I(freq2^2) -4.696551 1.283642 -3.659 0.000414 ***
invadedTRUE 0.063377 0.339992 0.186 0.852518
freq1:freq2 -1.075665 2.512512 -0.428 0.669520
freq1:I(freq2^2) 0.286132 2.141263 0.134 0.893977
freq1:invadedTRUE -0.176342 0.567146 -0.311 0.756529
freq2:invadedTRUE 0.015445 2.130085 0.007 0.994230
I(freq2^2):invadedTRUE 0.018287 1.815344 0.010 0.991983
freq1:freq2:invadedTRUE 0.873597 3.553228 0.246 0.806315
freq1:I(freq2^2):invadedTRUE -0.786470 3.028203 -0.260 0.795640
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.2484 on 96 degrees of freedom
Multiple R-squared: 0.5896, Adjusted R-squared: 0.5426
F-statistic: 12.54 on 11 and 96 DF, p-value: 2.428e-14# Backward model selection: which terms can be dropped?
# step() shows that all invasion terms drop but some freq1 terms stay.
lm.novel.alpha1.step = step(lm.novel.alpha1.v1, test = "F")Start: AIC=-289.54
alpha1 ~ freq1 + freq2 + I(freq2^2) + invaded + freq1:freq2 +
freq1:I(freq2^2) + freq1:invaded + freq2:invaded + I(freq2^2):invaded +
freq1:freq2:invaded + freq1:I(freq2^2):invaded
Df Sum of Sq RSS AIC F value Pr(>F)
- freq1:freq2:invaded 1 0.0037300 5.9275 -291.47 0.0604 0.8063
- freq1:I(freq2^2):invaded 1 0.0041622 5.9280 -291.46 0.0675 0.7956
<none> 5.9238 -289.54
Step: AIC=-291.47
alpha1 ~ freq1 + freq2 + I(freq2^2) + invaded + freq1:freq2 +
freq1:I(freq2^2) + freq1:invaded + freq2:invaded + I(freq2^2):invaded +
freq1:I(freq2^2):invaded
Df Sum of Sq RSS AIC F value Pr(>F)
- freq1:I(freq2^2):invaded 1 0.0016770 5.9292 -293.44 0.0274 0.8688
- freq2:invaded 1 0.0056607 5.9332 -293.37 0.0926 0.7615
- freq1:freq2 1 0.0079793 5.9355 -293.33 0.1306 0.7186
<none> 5.9275 -291.47
Step: AIC=-293.44
alpha1 ~ freq1 + freq2 + I(freq2^2) + invaded + freq1:freq2 +
freq1:I(freq2^2) + freq1:invaded + freq2:invaded + I(freq2^2):invaded
Df Sum of Sq RSS AIC F value Pr(>F)
- freq1:I(freq2^2) 1 0.0003088 5.9295 -295.44 0.0051 0.9432
- I(freq2^2):invaded 1 0.0052770 5.9345 -295.35 0.0872 0.7684
- freq2:invaded 1 0.0056607 5.9349 -295.34 0.0936 0.7603
- freq1:freq2 1 0.0079793 5.9372 -295.30 0.1319 0.7173
- freq1:invaded 1 0.0138340 5.9430 -295.19 0.2287 0.6336
<none> 5.9292 -293.44
Step: AIC=-295.44
alpha1 ~ freq1 + freq2 + I(freq2^2) + invaded + freq1:freq2 +
freq1:invaded + freq2:invaded + I(freq2^2):invaded
Df Sum of Sq RSS AIC F value Pr(>F)
- I(freq2^2):invaded 1 0.00528 5.9348 -297.34 0.0881 0.7672
- freq2:invaded 1 0.00566 5.9352 -297.33 0.0945 0.7592
- freq1:invaded 1 0.01383 5.9433 -297.19 0.2310 0.6319
<none> 5.9295 -295.44
- freq1:freq2 1 1.40538 7.3349 -274.46 23.4645 4.706e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Step: AIC=-297.34
alpha1 ~ freq1 + freq2 + I(freq2^2) + invaded + freq1:freq2 +
freq1:invaded + freq2:invaded
Df Sum of Sq RSS AIC F value Pr(>F)
- freq2:invaded 1 0.0010 5.9358 -299.32 0.0173 0.8955
- freq1:invaded 1 0.0138 5.9486 -299.09 0.2331 0.6303
<none> 5.9348 -297.34
- freq1:freq2 1 1.4054 7.3402 -276.39 23.6804 4.249e-06 ***
- I(freq2^2) 1 4.0444 9.9792 -243.22 68.1478 6.403e-13 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Step: AIC=-299.32
alpha1 ~ freq1 + freq2 + I(freq2^2) + invaded + freq1:freq2 +
freq1:invaded
Df Sum of Sq RSS AIC F value Pr(>F)
- freq1:invaded 1 0.0138 5.9496 -301.07 0.2354 0.6286
<none> 5.9358 -299.32
- freq1:freq2 1 1.4054 7.3412 -278.37 23.9131 3.810e-06 ***
- I(freq2^2) 1 4.0444 9.9802 -245.21 68.8174 4.932e-13 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Step: AIC=-301.07
alpha1 ~ freq1 + freq2 + I(freq2^2) + invaded + freq1:freq2
Df Sum of Sq RSS AIC F value Pr(>F)
- invaded 1 0.0680 6.0176 -301.84 1.1655 0.2829
<none> 5.9496 -301.07
- freq1:freq2 1 1.4054 7.3550 -280.17 24.0937 3.492e-06 ***
- I(freq2^2) 1 4.0444 9.9941 -247.06 69.3371 3.974e-13 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Step: AIC=-301.84
alpha1 ~ freq1 + freq2 + I(freq2^2) + freq1:freq2
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 6.0176 -301.84
- freq1:freq2 1 1.4054 7.4230 -281.18 24.055 3.509e-06 ***
- I(freq2^2) 1 4.0444 10.0621 -248.32 69.226 3.877e-13 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# summary() shows no significant effect of freq1.
summary(lm.novel.alpha1.step)
Call:
lm(formula = alpha1 ~ freq1 + freq2 + I(freq2^2) + freq1:freq2,
data = data.novel)
Residuals:
Min 1Q Median 3Q Max
-0.62383 -0.12454 0.00766 0.14219 0.62215
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.71882 0.11489 6.257 9.11e-09 ***
freq1 -0.07739 0.09339 -0.829 0.409
freq2 5.85144 0.67178 8.710 5.38e-14 ***
I(freq2^2) -4.73650 0.56928 -8.320 3.88e-13 ***
freq1:freq2 -0.76403 0.15578 -4.905 3.51e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.2417 on 103 degrees of freedom
Multiple R-squared: 0.5831, Adjusted R-squared: 0.5669
F-statistic: 36.01 on 4 and 103 DF, p-value: < 2.2e-16# add1() shows that freq1 should nonetheless be included in the model.
update(lm.novel.alpha1.full, ~ freq2 + I(freq2^2)) %>%
add1(~ . + freq1, test="F")# Generate the model predictions.
pred.novel.alpha1.step = predict(lm.novel.alpha1.step, newdata = newdat2, se.fit = T) %>%
magrittr::extract(c("fit", "se.fit")) %>%
as.data.frame() %>%
cbind(newdat2, .) %>%
transform(lwr = fit - se.fit, upr = fit + se.fit)
# Store the plot.
fig.novel.alpha1_v02 = ggplot(data = data.novel,
aes(x = freq2, y = alpha1)) +
geom_point(aes(color = as.factor(freq1), shape = invaded),
position = position_jitterdodge(jitter.width = 0.01, jitter.height = 0,
dodge.width = 0.5, seed=1), alpha=0.5, size=2.5, stroke=1.1) +
scale_x_continuous(trans = "log2",
breaks=1/c(1,4,8), labels=c("1\nHigh", "4", "8\nLow")) +
scale_shape_manual(values=c(1, 19), name = "Invasion", label = c("Uninvaded", "Invaded")) +
scale_color_brewer(palette = "Dark2", name = "Past freq.",
label = c("8 Low", "4", "1 High")) +
xlab("Novel freq. (every x days)") +
ylab("Effective number of species") -
geom_line(data = subset(pred.novel.alpha1.step, freq1 %in% c(0.125, 0.25, 1)),
aes(y=fit, color=as.factor(freq1)), alpha=1, linewidth=1) -
geom_ribbon(data = subset(pred.novel.alpha1.step, freq1 %in% c(0.125, 0.25, 1)),
aes(y=NULL, ymin=lwr, ymax=upr, fill=as.factor(freq1)), alpha=0.2) +
scale_fill_brewer(palette = "Dark2", guide = NULL) +
guides(col = guide_legend(reverse = TRUE))2.3 Community dissimilarity before and after the regime change.
Here, we analyze how community dissimilarity index changes with disturbance regimes.
Results: The best model is:
beta.past ~ freq2 + I(freq2^2) + freq1 + I(freq1^2) + freq2:freq1
+ freq2:I(freq1^2) + I(freq2^2):freq1 +
I(freq22):I(freq12).
# Fit a model with unimodal relationships against past and novel regimes.
fit.novel.beta1 = lm(beta.past ~ (freq2 + I(freq2^2)) * (freq1 + I(freq1^2)) * invaded,
data = data.novel)
# Backward model selection. Which terms can be dropped?
fit.novel.beta1.step = step(fit.novel.beta1, test = "F")Start: AIC=-547.59
beta.past ~ (freq2 + I(freq2^2)) * (freq1 + I(freq1^2)) * invaded
Df Sum of Sq RSS AIC F value Pr(>F)
- freq2:I(freq1^2):invaded 1 4.1857e-05 0.48609 -549.58 0.0078 0.9300
- freq2:freq1:invaded 1 5.6786e-05 0.48611 -549.57 0.0105 0.9186
- I(freq2^2):I(freq1^2):invaded 1 1.2057e-04 0.48617 -549.56 0.0223 0.8816
- I(freq2^2):freq1:invaded 1 1.5384e-04 0.48621 -549.55 0.0285 0.8663
<none> 0.48605 -547.59
Step: AIC=-549.58
beta.past ~ freq2 + I(freq2^2) + freq1 + I(freq1^2) + invaded +
freq2:freq1 + freq2:I(freq1^2) + I(freq2^2):freq1 + I(freq2^2):I(freq1^2) +
freq2:invaded + I(freq2^2):invaded + freq1:invaded + I(freq1^2):invaded +
freq2:freq1:invaded + I(freq2^2):freq1:invaded + I(freq2^2):I(freq1^2):invaded
Df Sum of Sq RSS AIC F value Pr(>F)
- freq2:freq1:invaded 1 0.000147 0.48624 -551.54 0.0275 0.8686
- I(freq2^2):freq1:invaded 1 0.002215 0.48831 -551.09 0.4147 0.5212
- I(freq2^2):I(freq1^2):invaded 1 0.002535 0.48863 -551.01 0.4746 0.4926
<none> 0.48609 -549.58
- freq2:I(freq1^2) 1 0.109079 0.59517 -529.71 20.4202 1.864e-05
- freq2:freq1:invaded
- I(freq2^2):freq1:invaded
- I(freq2^2):I(freq1^2):invaded
<none>
- freq2:I(freq1^2) ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Step: AIC=-551.54
beta.past ~ freq2 + I(freq2^2) + freq1 + I(freq1^2) + invaded +
freq2:freq1 + freq2:I(freq1^2) + I(freq2^2):freq1 + I(freq2^2):I(freq1^2) +
freq2:invaded + I(freq2^2):invaded + freq1:invaded + I(freq1^2):invaded +
I(freq2^2):freq1:invaded + I(freq2^2):I(freq1^2):invaded
Df Sum of Sq RSS AIC F value Pr(>F)
- freq2:invaded 1 0.001461 0.48770 -553.22 0.2765 0.6003
- I(freq2^2):I(freq1^2):invaded 1 0.002535 0.48878 -552.98 0.4796 0.4903
- I(freq2^2):freq1:invaded 1 0.002955 0.48920 -552.89 0.5591 0.4565
<none> 0.48624 -551.54
- freq2:I(freq1^2) 1 0.109079 0.59532 -531.68 20.6384 1.680e-05
- freq2:freq1 1 0.123117 0.60936 -529.17 23.2945 5.502e-06
- freq2:invaded
- I(freq2^2):I(freq1^2):invaded
- I(freq2^2):freq1:invaded
<none>
- freq2:I(freq1^2) ***
- freq2:freq1 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Step: AIC=-553.22
beta.past ~ freq2 + I(freq2^2) + freq1 + I(freq1^2) + invaded +
freq2:freq1 + freq2:I(freq1^2) + I(freq2^2):freq1 + I(freq2^2):I(freq1^2) +
I(freq2^2):invaded + freq1:invaded + I(freq1^2):invaded +
I(freq2^2):freq1:invaded + I(freq2^2):I(freq1^2):invaded
Df Sum of Sq RSS AIC F value Pr(>F)
- I(freq2^2):I(freq1^2):invaded 1 0.002535 0.49024 -554.66 0.4834 0.4886
- I(freq2^2):freq1:invaded 1 0.002955 0.49066 -554.57 0.5634 0.4548
<none> 0.48770 -553.22
- freq2:I(freq1^2) 1 0.109079 0.59678 -533.42 20.8002 1.552e-05
- freq2:freq1 1 0.123117 0.61082 -530.91 23.4771 5.036e-06
- I(freq2^2):I(freq1^2):invaded
- I(freq2^2):freq1:invaded
<none>
- freq2:I(freq1^2) ***
- freq2:freq1 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Step: AIC=-554.66
beta.past ~ freq2 + I(freq2^2) + freq1 + I(freq1^2) + invaded +
freq2:freq1 + freq2:I(freq1^2) + I(freq2^2):freq1 + I(freq2^2):I(freq1^2) +
I(freq2^2):invaded + freq1:invaded + I(freq1^2):invaded +
I(freq2^2):freq1:invaded
Df Sum of Sq RSS AIC F value Pr(>F)
- I(freq1^2):invaded 1 0.000025 0.49026 -556.65 0.0048 0.9451
- I(freq2^2):freq1:invaded 1 0.002187 0.49243 -556.18 0.4193 0.5189
<none> 0.49024 -554.66
- I(freq2^2):I(freq1^2) 1 0.100983 0.59122 -536.43 19.3628 2.852e-05 ***
- freq2:I(freq1^2) 1 0.109079 0.59932 -534.96 20.9152 1.462e-05 ***
- freq2:freq1 1 0.123117 0.61336 -532.46 23.6068 4.711e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Step: AIC=-556.65
beta.past ~ freq2 + I(freq2^2) + freq1 + I(freq1^2) + invaded +
freq2:freq1 + freq2:I(freq1^2) + I(freq2^2):freq1 + I(freq2^2):I(freq1^2) +
I(freq2^2):invaded + freq1:invaded + I(freq2^2):freq1:invaded
Df Sum of Sq RSS AIC F value Pr(>F)
- I(freq2^2):freq1:invaded 1 0.002187 0.49245 -558.17 0.4238 0.5166
<none> 0.49026 -556.65
- I(freq2^2):I(freq1^2) 1 0.100983 0.59125 -538.43 19.5678 2.585e-05 ***
- freq2:I(freq1^2) 1 0.109079 0.59934 -536.96 21.1366 1.317e-05 ***
- freq2:freq1 1 0.123117 0.61338 -534.46 23.8568 4.193e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Step: AIC=-558.17
beta.past ~ freq2 + I(freq2^2) + freq1 + I(freq1^2) + invaded +
freq2:freq1 + freq2:I(freq1^2) + I(freq2^2):freq1 + I(freq2^2):I(freq1^2) +
I(freq2^2):invaded + freq1:invaded
Df Sum of Sq RSS AIC F value Pr(>F)
- freq1:invaded 1 0.000085 0.49253 -560.15 0.0165 0.8980
- I(freq2^2):invaded 1 0.003879 0.49633 -559.33 0.7563 0.3867
<none> 0.49245 -558.17
- I(freq2^2):I(freq1^2) 1 0.100983 0.59343 -540.03 19.6860 2.434e-05 ***
- freq2:I(freq1^2) 1 0.109079 0.60153 -538.56 21.2642 1.234e-05 ***
- I(freq2^2):freq1 1 0.114343 0.60679 -537.62 22.2905 7.983e-06 ***
- freq2:freq1 1 0.123117 0.61557 -536.07 24.0008 3.900e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Step: AIC=-560.15
beta.past ~ freq2 + I(freq2^2) + freq1 + I(freq1^2) + invaded +
freq2:freq1 + freq2:I(freq1^2) + I(freq2^2):freq1 + I(freq2^2):I(freq1^2) +
I(freq2^2):invaded
Df Sum of Sq RSS AIC F value Pr(>F)
- I(freq2^2):invaded 1 0.003879 0.49641 -561.31 0.764 0.3842
<none> 0.49253 -560.15
- I(freq2^2):I(freq1^2) 1 0.100983 0.59352 -542.01 19.888 2.210e-05 ***
- freq2:I(freq1^2) 1 0.109079 0.60161 -540.55 21.482 1.113e-05 ***
- I(freq2^2):freq1 1 0.114343 0.60688 -539.61 22.519 7.169e-06 ***
- freq2:freq1 1 0.123117 0.61565 -538.06 24.247 3.478e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Step: AIC=-561.31
beta.past ~ freq2 + I(freq2^2) + freq1 + I(freq1^2) + invaded +
freq2:freq1 + freq2:I(freq1^2) + I(freq2^2):freq1 + I(freq2^2):I(freq1^2)
Df Sum of Sq RSS AIC F value Pr(>F)
- invaded 1 0.000004 0.49642 -563.31 0.0008 0.978
<none> 0.49641 -561.31
- I(freq2^2):I(freq1^2) 1 0.100983 0.59740 -543.31 19.9356 2.145e-05 ***
- freq2:I(freq1^2) 1 0.109079 0.60549 -541.86 21.5339 1.078e-05 ***
- I(freq2^2):freq1 1 0.114343 0.61076 -540.92 22.5731 6.929e-06 ***
- freq2:freq1 1 0.123117 0.61953 -539.38 24.3052 3.352e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Step: AIC=-563.31
beta.past ~ freq2 + I(freq2^2) + freq1 + I(freq1^2) + freq2:freq1 +
freq2:I(freq1^2) + I(freq2^2):freq1 + I(freq2^2):I(freq1^2)
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 0.49642 -563.31
- I(freq2^2):I(freq1^2) 1 0.10098 0.59740 -545.31 20.139 1.947e-05 ***
- freq2:I(freq1^2) 1 0.10908 0.60550 -543.85 21.753 9.713e-06 ***
- I(freq2^2):freq1 1 0.11434 0.61076 -542.92 22.803 6.218e-06 ***
- freq2:freq1 1 0.12312 0.61953 -541.38 24.553 2.987e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Resulting model.
summary(fit.novel.beta1.step)
Call:
lm(formula = beta.past ~ freq2 + I(freq2^2) + freq1 + I(freq1^2) +
freq2:freq1 + freq2:I(freq1^2) + I(freq2^2):freq1 + I(freq2^2):I(freq1^2),
data = data.novel)
Residuals:
Min 1Q Median 3Q Max
-0.222236 -0.026480 -0.005259 0.014133 0.217084
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.7534 0.1432 -5.262 8.28e-07 ***
freq2 6.3613 0.8971 7.091 2.00e-10 ***
I(freq2^2) -5.3221 0.7645 -6.961 3.72e-10 ***
freq1 3.6059 0.8971 4.020 0.000114 ***
I(freq1^2) -2.9338 0.7645 -3.837 0.000219 ***
freq2:freq1 -27.8486 5.6202 -4.955 2.99e-06 ***
freq2:I(freq1^2) 22.3397 4.7898 4.664 9.71e-06 ***
I(freq2^2):freq1 22.8724 4.7898 4.775 6.22e-06 ***
I(freq2^2):I(freq1^2) -18.3186 4.0820 -4.488 1.95e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.07081 on 99 degrees of freedom
Multiple R-squared: 0.6195, Adjusted R-squared: 0.5887
F-statistic: 20.15 on 8 and 99 DF, p-value: < 2.2e-16# Generate the model predictions.
pred.novel.beta1.step = predict(fit.novel.beta1.step, newdata = newdat2, se.fit = T, type = "response") %>%
magrittr::extract(c("fit", "se.fit")) %>%
as.data.frame() %>%
cbind(newdat2, .) %>%
transform(lwr = fit - se.fit, upr = fit + se.fit)
# Store the plot.
fig.novel.beta1.v02 = ggplot(data = data.novel,
aes(x = freq2, y = beta.past)) +
geom_point(aes(color = as.factor(freq1), shape = invaded),
position = position_jitterdodge(jitter.width = 0.01, jitter.height = 0,
dodge.width = 0.5, seed=1), alpha=0.5, size=2.5, stroke=1.1) +
scale_x_continuous(trans = "log2",
breaks=1/c(1,4,8), labels=c("1\nHigh", "4", "8\nLow")) +
scale_shape_manual(values=c(1, 19), name = "Invasion", label = c("Uninvaded", "Invaded")) +
scale_color_brewer(palette = "Dark2", name = "Past freq.",
label = c("8 Low", "4", "1 High")) +
xlab("Novel freq. (every x days)") +
ylab("Community dissimilarity") -
geom_line(data = subset(pred.novel.beta1.step, freq1 %in% c(0.125, 0.25, 1)),
aes(y=fit, color=as.factor(freq1)), alpha=1, linewidth=1) -
geom_ribbon(data = subset(pred.novel.beta1.step, freq1 %in% c(0.125, 0.25, 1)),
aes(y=NULL, ymin=lwr, ymax=upr, fill=as.factor(freq1)), alpha=0.2) +
scale_fill_brewer(palette = "Dark2", guide = NULL) +
guides(col = guide_legend(reverse = TRUE))
# Plot and save all diversity analyses.
common.legend = get_legend(fig.novel.alpha1_v02)
ggarrange(
plotlist = list(
fig.past.alpha1 + ylim(0.92, 2.6),
fig.novel.alpha1_v02 + ylim(0.92, 2.6) + theme(legend.position = 'none'),
fig.novel.beta1.v02 + theme(legend.position = 'none'),
common.legend
),
labels = list("A", "B", "C", ""))
ggsave("./plot/Day17_beta_v02.png", width = 30, height = 25, units = "cm")3 Effects of disturbance regimes and diversity on invasibility
Here, we test the diversity-invasability relationship.
Hypothesis: Past regimes affect invasion success. Specifically, more diverse community would suppress invaders, so invasion success would follow similar pattern as the disturbance-diversity relationship.
Results: The best model is:
inv.raw ~ freq1 + I(freq1^2) + freq2 + alpha1.past +
I(freq1^2):freq2 + alpha1.past * (I(freq1^2) + freq2 + I(freq1^2):freq2)
+ offset(log(total.dilution)).
First, we tested whether both disturbance regimes are necessary to understand invasion outcomes. They are:
add1shows that bothfreq1andfreq2should be included, but not bothI(freq1^2)andI(freq2^2).Second, we tested whether including diversity helps explain invasions outcomes. It does. AICs show that including diversity always improves the model.
Third, we performed a stepwise backward model selection to assess how disturbance regimes and diversity jointly affect invasion outcomes. Our analyses show that both past and novel regimes qualitatively alter diversity-invasion relationship.
# We first test whether both disturbance regimes should be included in the model, using
# simple models WITHOUT interactions.
# Do we need both freq1 and freq2 terms? Sequential add1()'s suggest that both freq1
# and freq2 should be included, but not both I(freq1^2) and I(freq2^2).
data.legacy = data.past %>%
select(sample.code, alpha1) %>%
rename(preinv.sample.code = sample.code) %>%
merge(subset(data.novel, invaded == T & sample.day == 17), .,
by = c("preinv.sample.code"),
suffixes = c(".now", ".past"))
# This null model only has freq1 terms.
glm.inv.null.v1 = glm(
inv.raw ~ (freq1 + I(freq1^2)) + offset(log(total.dilution)),
data = data.novel,
# subset = (invaded==T & sample.day==17),
family = "poisson")
# add1() shows that freq2 should be included.
add1(glm.inv.null.v1, ~. + freq2, test="Chisq")glm.inv.null.v1.add = update(glm.inv.null.v1, ~. + freq2)
# add1() shows that I(freq2^2) need not be included.
add1(glm.inv.null.v1.add, ~. + I(freq2^2), test = "Chisq")# This null model only has freq2 terms.
glm.inv.null.v2 = glm(
inv.raw ~ (freq2 + I(freq2^2)) + offset(log(total.dilution)),
data = data.novel,
# subset = (invaded==T & sample.day==17),
family = "poisson")
# add1() shows that freq1 should be included.
add1(glm.inv.null.v2, ~. + freq1, test="Chisq")glm.inv.null.v2.add = update(glm.inv.null.v2, ~. + freq1)
# add1() shows that I(freq1^2) need not be included.
add1(glm.inv.null.v2.add, ~. + I(freq1^2), test = "Chisq")# Using more complex models that include both regimes, we compared how the resident diversity helps explain invasion outcomes. Resident diversity always increase model fit.
# Assume that all terms interact with one another.
glm.inv.full = glm(
inv.raw ~ (freq1 + I(freq1^2)) * (freq2 + I(freq2^2)) + offset(log(total.dilution)),
data = data.legacy,
# subset = (invaded==T & sample.day==17),
family = "poisson")
# ... and with diversity
glm.inv.pastdiv.full = update(glm.inv.full, ~ . * alpha1.past)
# Assume linear freq1
glm.inv.lin1 = update(glm.inv.full,
~ freq1 * (freq2 + I(freq2^2)) + offset(log(total.dilution)))
# ... and with diversity
glm.inv.pastdiv.lin1 = update(glm.inv.lin1, ~ . * alpha1.past )
# Assume linear freq2
glm.inv.lin2 = update(glm.inv.full,
~ (freq1 + I(freq1^2)) * freq2 + offset(log(total.dilution)))
# ... and with diversity
glm.inv.pastdiv.lin2 = update(glm.inv.lin2, ~ . * alpha1.past )
# Assume linear freq1 and freq2
glm.inv.lin12 = update(glm.inv.full,
~ freq1 * freq2 + offset(log(total.dilution)))
# ... and with diversity
glm.inv.pastdiv.lin12 = update(glm.inv.lin12, ~ . * alpha1.past )
# Diversity ALWAYS lowers the scores and improves the models.
AIC(glm.inv.full # (freq1 + I(freq1^2)) * (freq2 + I(freq2^2))
, glm.inv.pastdiv.full
, glm.inv.lin1 # freq1 * (freq2 + I(freq2^2))
, glm.inv.pastdiv.lin1
, glm.inv.lin2 # (freq1 + I(freq1^2)) * freq2
, glm.inv.pastdiv.lin2
, glm.inv.lin12 # freq1 * freq2
, glm.inv.pastdiv.lin12
)# Backward model selection
drop1(glm.inv.pastdiv.full, test = "Chisq")glm.inv.pastdiv.full.drop = update(glm.inv.pastdiv.full,
. ~ . - I(freq1^2):I(freq2^2):alpha1.past)
drop1(glm.inv.pastdiv.full.drop, test = "Chisq")glm.inv.pastdiv.full.drop = update(glm.inv.pastdiv.full.drop,
. ~ . - freq1:I(freq2^2):alpha1.past)
drop1(glm.inv.pastdiv.full.drop, test = "Chisq")glm.inv.pastdiv.full.drop = update(glm.inv.pastdiv.full.drop,
. ~ . - I(freq2^2):alpha1.past)
drop1(glm.inv.pastdiv.full.drop, test = "Chisq")# No significant difference, so still dropping.
glm.inv.pastdiv.full.drop = update(glm.inv.pastdiv.full.drop,
. ~ . - freq1:freq2:alpha1.past )
drop1(glm.inv.pastdiv.full.drop, test = "Chisq")glm.inv.pastdiv.full.drop = update(glm.inv.pastdiv.full.drop,
. ~ . - freq1:alpha1.past)
drop1(glm.inv.pastdiv.full.drop, test = "Chisq")glm.inv.pastdiv.full.drop = update(glm.inv.pastdiv.full.drop,
. ~ . - freq1:freq2)
drop1(glm.inv.pastdiv.full.drop, test = "Chisq")glm.inv.pastdiv.full.drop = update(glm.inv.pastdiv.full.drop,
. ~ . - I(freq1^2):I(freq2^2))
drop1(glm.inv.pastdiv.full.drop, test = "Chisq")glm.inv.pastdiv.full.drop = update(glm.inv.pastdiv.full.drop,
. ~ . - freq1:I(freq2^2))
drop1(glm.inv.pastdiv.full.drop, test = "Chisq")glm.inv.pastdiv.full.drop = update(glm.inv.pastdiv.full.drop,
. ~ . - I(freq2^2))
drop1(glm.inv.pastdiv.full.drop, test = "Chisq")summary(glm.inv.pastdiv.full.drop)
Call:
glm(formula = inv.raw ~ freq1 + I(freq1^2) + freq2 + alpha1.past +
I(freq1^2):freq2 + I(freq1^2):alpha1.past + freq2:alpha1.past +
I(freq1^2):freq2:alpha1.past + offset(log(total.dilution)),
family = "poisson", data = data.legacy)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 10.567 1.512 6.989 2.77e-12 ***
freq1 -9.102 4.590 -1.983 0.047376 *
I(freq1^2) 20.018 7.126 2.809 0.004970 **
freq2 13.393 2.669 5.019 5.21e-07 ***
alpha1.past 4.376 1.305 3.353 0.000801 ***
I(freq1^2):freq2 -26.975 8.045 -3.353 0.000799 ***
I(freq1^2):alpha1.past -11.244 5.559 -2.023 0.043088 *
freq2:alpha1.past -9.428 2.121 -4.445 8.78e-06 ***
I(freq1^2):freq2:alpha1.past 23.766 7.468 3.182 0.001461 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 107.939 on 53 degrees of freedom
Residual deviance: 53.762 on 45 degrees of freedom
AIC: 152.65
Number of Fisher Scoring iterations: 6# Calculate pseudo-R2 for Poisson regression.
pR2(glm.inv.pastdiv.full.drop)fitting null model for pseudo-r2 llh llhNull G2 McFadden r2ML r2CU
-67.3270533 -120.7720540 106.8900014 0.4425279 0.8618540 0.8718036 # Generate the model predictions.
pred.inv.pastdiv = predict.glm(glm.inv.pastdiv.full.drop,
newdata = subset(newdat3 %>% rename(alpha1.past = alpha1),
freq1 %in% c(0.125, 0.25, 1) & freq2 %in% c(0.125, 0.25, 1) & invaded==T),
type="link", se.fit=T) %>%
magrittr::extract(c("fit", "se.fit")) %>%
as.data.frame() %>%
cbind(subset(newdat3 %>% rename(alpha1.past = alpha1),
freq1 %in% c(0.125, 0.25, 1) & freq2 %in% c(0.125, 0.25, 1) & invaded==T), .) %>%
transform(upr=fit+se.fit, lwr=fit-se.fit) %>%
transform(r.fit=exp(fit), r.upr=exp(upr), r.lwr=exp(lwr))
# Store the plot for the data points.
fig.inv.pastdiv = ggplot(data=subset(data.legacy, invaded == T & sample.day == 17),
aes(y = inv.raw/total.dilution, x = alpha1.past, color = as.factor(freq1))) +
xlab("Effective number of species") +
ylab("Density of invaders (CFU/mL)") +
theme(axis.text.x = element_text(size = 10)) +
scale_y_continuous(label = scientific_10) +
coord_cartesian(ylim = c(0, 4.5e7)) +
scale_color_brewer(palette = "Dark2", name = "Past freq.",
breaks=1/c(8,4,1), label = c("8 Low", "4", "1 High"))+
scale_fill_brewer(palette = "Dark2", guide = NULL) +
geom_jitter(alpha=0.5, size=2.5, stroke=1.1, width = 0.05, height = 0) +
facet_wrap(~ freq2, labeller = as_labeller(freq2.labels))
# Plot and save the data points and the model predictions.
fig.inv.pastdiv -
geom_line(data=pred.inv.pastdiv, aes(y=r.fit),
size=0.75, linetype=1) -
geom_ribbon(data=pred.inv.pastdiv,
aes(y=NULL, ymin=r.lwr, ymax=r.upr, fill = as.factor(freq1)),
alpha=0.2, color=NA) +
guides(col = guide_legend(reverse = TRUE))
ggsave("./plot/Day17_inv_pastdiv.png", width = 30, units = "cm")Saving 30 x 12.7 cm image4 Counterfactuals
What happens if we don’t consider disturbance regime change?
# I(num1 + num2) is the total number of passages.
fit.counter.alpha1 = lm(alpha1 ~ as.factor(I(num1 + num2)) * inc.freq * invaded
, data = subset(data.novel, num1 != num2 & sample.day == 17))
Anova(fit.counter.alpha1)fit.counter.alpha1.v02 = update(fit.counter.alpha1, ~ as.factor(I(num1 + num2)) * invaded)
# Removing inc.freq makes the model worse.
AIC(fit.counter.alpha1, fit.counter.alpha1.v02)anova(fit.counter.alpha1, fit.counter.alpha1.v02, test = "F")rc.avg = subset(data.novel, num1 != num2 & sample.day == 17) %>%
group_by(num.total) %>%
summarize(avg = mean(alpha1, na.rm = T),
se = sd(alpha1, na.rm = T) / sqrt(n()))
fig.counter.alpha1 = ggplot(
# Focus on samples that underwent regime change.
data = subset(data.novel, num1 != num2 & sample.day == 17),
aes(x = paste(1/freq1, 1/freq2), y = alpha1, group = inc.freq)) +
facet_wrap(~ num.total , nrow = 1, scales = "free_x", labeller = as_labeller(num.total.labels)) +
scale_shape_manual(values = c(25, 24),
name = "Frequency change", labels = c("Decreased", "Increased")) +
scale_color_manual(values = c(1, 2),
name = "Invasion", labels = c("Uninvaded", "Invaded")) +
scale_fill_manual(guide = NULL, values = c("white", 2)) +
scale_x_discrete(
labels = c(
"4 8" = "4 \u2192 8",
"8 4" = "8 \u2192 4",
"1 8" = "1 \u2192 8",
"8 1" = "8 \u2192 1",
"1 4" = "1 \u2192 4",
"4 1" = "4 \u2192 1"
), drop = TRUE)+
xlab("Regime change") +
ylab("Effective number of species") -
# Individual data.
geom_point(aes(shape = inc.freq, color = invaded, fill = invaded),
position = position_jitterdodge(jitter.width = 0.1, jitter.height = 0, dodge.width = 0.4, seed=1),
alpha = 0.3, size = 1.75, stroke = 1) +
# Stats by frequency change direction.
geom_point(stat="summary", fun="mean", aes(shape = inc.freq),
position=position_dodge(width = 0.5), size = 2.5, stroke = 1, fill = 1) +
geom_errorbar(stat="summary", fun.data="mean_se", width = 0, size = 1,
position=position_dodge(width=0.5)) +
# Stats ignoring frequency change.
geom_rect(data = rc.avg,
aes(xmin = -Inf, xmax = Inf,
ymin = avg - se, ymax = avg + se),
inherit.aes = FALSE, alpha = 0.2) +
geom_hline(data = rc.avg, aes(yintercept = avg),
alpha = 1, linewidth = 1, inherit.aes = F) +
theme(
legend.title = element_text(size = 14),
legend.text = element_text(size = 14),
axis.text.x = element_text(size = 12))
fig.counter.alpha1
ggsave("./plot/rc_agnostic_alpha.png")Saving 7 x 5 in image5 Supporting analyses
5.1 ANOVA on resident diversity with unordered factors.
Here, we further test whether 1. invasion doesn’t affect diversity, and 2. past regimes improve model fit.
# anova() confirms that invasion really doesn't matter, but freq1 does.
aov.novel.alpha1 = aov(alpha1 ~ as.factor(freq1) * as.factor(freq2) * invaded,
data = data.novel)
anova(aov.novel.alpha1)# anova() on a simpler model without invasion again shows that freq1 matters:
aov.novel.alpha1.v2 = update(aov.novel.alpha1, . ~ as.factor(freq1) * as.factor(freq2))
anova(aov.novel.alpha1.v2)5.2 Sensitivity of the diversity-invasibility relationships to model structure.
Here, we assess how our results on diversity-invasibility relationships change when we analyze a simpler model.
summary(glm.inv.pastdiv.lin12)
Call:
glm(formula = inv.raw ~ freq1 + freq2 + alpha1.past + freq1:freq2 +
freq1:alpha1.past + freq2:alpha1.past + freq1:freq2:alpha1.past +
offset(log(total.dilution)), family = "poisson", data = data.legacy)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 8.617 2.370 3.636 0.000277 ***
freq1 13.193 6.310 2.091 0.036539 *
freq2 17.213 3.597 4.785 1.71e-06 ***
alpha1.past 5.354 1.921 2.787 0.005314 **
freq1:freq2 -29.185 9.237 -3.159 0.001581 **
freq1:alpha1.past -12.565 5.917 -2.123 0.033724 *
freq2:alpha1.past -12.954 3.003 -4.314 1.61e-05 ***
freq1:freq2:alpha1.past 25.787 8.520 3.027 0.002473 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 107.94 on 53 degrees of freedom
Residual deviance: 58.57 on 46 degrees of freedom
AIC: 155.46
Number of Fisher Scoring iterations: 6pred.inv.pastdiv.v02 = predict.glm(glm.inv.pastdiv.lin12,
newdata = subset(newdat3 %>% rename(alpha1.past = alpha1),
freq1 %in% c(0.125, 0.25, 1) & freq2 %in% c(0.125, 0.25, 1) & invaded==T),
type="link", se.fit=T) %>%
magrittr::extract(c("fit", "se.fit")) %>%
as.data.frame() %>%
cbind(subset(newdat3 %>% rename(alpha1.past = alpha1),
freq1 %in% c(0.125, 0.25, 1) & freq2 %in% c(0.125, 0.25, 1) & invaded==T), .) %>%
transform(upr=fit+se.fit, lwr=fit-se.fit) %>%
transform(r.fit=exp(fit), r.upr=exp(upr), r.lwr=exp(lwr))
fig.inv.pastdiv -
geom_line(data=pred.inv.pastdiv.v02, aes(y=r.fit),
size=0.75, linetype=1) -
geom_ribbon(data=pred.inv.pastdiv.v02,
aes(y=NULL, ymin=r.lwr, ymax=r.upr, fill = as.factor(freq1)),
alpha=0.2, color=NA) +
guides(col = guide_legend(reverse = TRUE))
5.3 Effects of regime changes alone on invasibility.
Here, we test whether regime change alone can help explain the differences in invasion outcomes.
# This model focuses on the effects of regime change without the diversity covariate.
# We are interested in understanding the effects of specific regime pairs, so the
# regimes are coded as factors.
glm.inv.lin12.factor = update(glm.inv.lin12,
. ~ as.factor(freq1) * as.factor(freq2) + offset(log(total.dilution)))
# Calculate the marginal means of the regime changes.
emmeans(glm.inv.lin12.factor,
specs = ~ as.factor(freq2) | as.factor(freq1),
type = "response") %>%
pairs(adjust = "none") %>%
as.data.frame() %>%
slice(-c(3, 5, 7)) %>%
mutate(p.adj = p.adjust(p.value, method = "holm"))# Assess the effects of dissimilarity index on invasibility.
glm.inv.beta = update(glm.inv.lin12,
. ~ beta.past + offset(log(total.dilution)))
summary(glm.inv.beta)
Call:
glm(formula = inv.raw ~ beta.past + offset(log(total.dilution)),
family = "poisson", data = data.legacy)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 15.0805 0.1219 123.746 < 2e-16 ***
beta.past 3.5346 0.9697 3.645 0.000267 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 107.939 on 53 degrees of freedom
Residual deviance: 97.665 on 52 degrees of freedom
AIC: 182.56
Number of Fisher Scoring iterations: 6LS0tCnRpdGxlOiAiQ29kZXMgZm9yOiBEaXN0dXJiYW5jZSByZWdpbWUgY2hhbmdlcyBsZWF2ZSBsb25nLWxhc3RpbmcgbGVnYWNpZXMgb24gYSBtaWNyb2JpYWwgY29tbXVuaXR54oCZcyBjb21wb3NpdGlvbiBhbmQKZnVuY3Rpb24uIgphdXRob3I6IEguIEluYW1pbmUsIEwuIExlYXIsIEEuIE1pbGxlciwgUy4gUm94YnVyZ2gsIEEuIEJ1Y2tsaW5nLCBLLiBTaGVhCmRhdGU6ICJVcGRhdGVkOiBgciBmb3JtYXQoU3lzLnRpbWUoKSwgJyVCICVkLCAlWScpYCIKb3V0cHV0OgogIGh0bWxfbm90ZWJvb2s6CiAgICBudW1iZXJfc2VjdGlvbnM6IHRydWUKICAgIHRvYzogdHJ1ZQogICAgdG9jX2Zsb2F0OiB0cnVlCiAgICB0aGVtZTogZmxhdGx5Ci0tLQoKYGBge3IgaW5jbHVkZSA9IEZBTFNFLCBjYWNoZSA9IEZBTFNFfQprbml0cjo6cmVhZF9jaHVuaygnLi9BbmFseXNpcy5SJykKYGBgCgojIFByZWFtYmxlCgpUaGUgY29kZXMgYmVsb3cgaW1wb3J0IHBhY2thZ2VzLCBkZWZpbmUgZnVuY3Rpb25zLCBhbmQgZm9ybWF0IHRoZSBkYXRhLiAKCmBgYHtyIFByZWFtYmxlfQpgYGAKCiMgRWZmZWN0cyBvZiBkaXN0dXJiYW5jZSByZWdpbWVzIG9uIGRpdmVyc2l0eQoKIyMgRGl2ZXJzaXR5IGJlZm9yZSB0aGUgcmVnaW1lIGNoYW5nZQoKSGVyZSwgd2UgYW5hbHl6ZSBob3cgZGl2ZXJzaXR5IGNoYW5nZXMgd2l0aCBkaXN0dXJiYW5jZSBmcmVxdWVuY3kgPGI+YmVmb3JlPC9iPiB0aGUgcmVnaW1lIGNoYW5nZS4KCjxiPkh5cG90aGVzaXM6PC9iPiAKRm9sbG93aW5nIHByZXZpb3VzIHN0dWRpZXMsIGRpdmVyc2l0eSBzaG91bGQgc2hvdyBhIHVuaW1vZGFsIHBhdHRlcm4gd2l0aCBpbmNyZWFzaW5nIGRpc3R1cmJhbmNlIGZyZXF1ZW5jeS4KCjxiPlJlc3VsdHM6PC9iPiAKRGlzdHVyYmFuY2UtZGl2ZXJzaXR5IHJlbGF0aW9uc2hpcCBiZWZvcmUgdGhlIHJlZ2ltZSBjaGFuZ2UgaGFzIGEgdW5pbW9kYWwgcmVsYXRpb25zaGlwOlwKPGNvZGU+YWxwaGExIH4gKGZyZXExICsgSShmcmVxMV4yKSk8L2NvZGU+XApXZSBoYXZlIHBvc2l0aXZlIDxjb2RlPmZyZXExPC9jb2RlPiBhbmQgbmVnYXRpdmUgPGNvZGU+SShmcmVxMV4yKTwvY29kZT4gY29lZmZpY2llbnRzLCBhbmQgbm8gY292YXJpYXRlcyBzaG91bGQgYmUgZHJvcHBlZC4gCgpgYGB7ciBQYXN0LmRpdmVyc2l0eX0KYGBgCgoKIyMgRGl2ZXJzaXR5IGFmdGVyIHRoZSByZWdpbWUgY2hhbmdlCgpIZXJlLCB3ZSBhbmFseXplIGhvdyBkaXZlcnNpdHkgY2hhbmdlcyB3aXRoIGRpc3R1cmJhbmNlIGZyZXF1ZW5jeSA8Yj5hZnRlcjwvYj4gdGhlIHJlZ2ltZSBjaGFuZ2UuCgo8Yj5IeXBvdGhlc2lzOjwvYj4gClNpbWlsYXIgdG8gdGhlIHJlc3VsdHMgYWJvdmUsIGRpc3R1cmJhbmNlLWRpdmVyc2l0eSByZWxhdGlvbnNoaXAgYWZ0ZXIgdGhlIG5vdmVsIHJlZ2ltZXMgd291bGQgYmUgdW5pbW9kYWwuIApGdXJ0aGVybW9yZSwgdGhlIGVmZmVjdHMgb2YgdGhlIHBhc3QgcmVnaW1lcyB3b3VsZCBoYXZlIGRpc2FwcGVhcmVkIHVuZGVyIHRoZSBub3ZlbCByZWdpbWVzLgoKCjxiPlJlc3VsdHM6PC9iPiAKVGhlIGJlc3QgbW9kZWwgaXM6XAo8Y29kZT5hbHBoYTEgfiBmcmVxMSArIGZyZXEyICsgSShmcmVxMl4yKSArIGZyZXExOmZyZXEyPC9jb2RlPi5cCldlIGhhdmUgcG9zaXRpdmUgPGNvZGU+ZnJlcTI8L2NvZGU+LCBuZWdhdGl2ZSA8Y29kZT5JKGZyZXEyXjIpPC9jb2RlPiwgYW5kIG5lZ2F0aXZlIDxjb2RlPmZyZXExOmZyZXEyPC9jb2RlPi4KPGNvZGU+ZnJlcTE8L2NvZGU+IGlzIG5vbi1zaWduaWZpY2FudC4KCjEuIFN0ZXAtd2lzZSBtb2RlbCBzZWxlY3Rpb24gb24gdGhlIGZ1bGwgbW9kZWwgc2hvd3MgdGhhdDoKICAgIC0gSW52YXNpb25zIGhhdmUgbm8gc2lnbmlmaWNhbnQgZWZmZWN0IG9uIHRoZSByZXNpZGVudCBkaXZlcnNpdHkuIAogICAgLSBQYXN0IHJlZ2ltZXMgYWZmZWN0IGRpdmVyc2l0eSB1bmRlciB0aGUgbm92ZWwgcmVnaW1lcywgYnV0IG9ubHkgdGhyb3VnaCB0aGUgaW50ZXJhY3Rpb24gdGVybSA8Y29kZT5mcmVxMjpmcmVxMTwvY29kZT4uCgoxLiA8Y29kZT5hZGQxPC9jb2RlPiBvbiBhIG1vZGVsIHdpdGhvdXQgYW55IDxjb2RlPmZyZXExPC9jb2RlPiB0ZXJtIGNvbmZpcm1zIHRoYXQgPGNvZGU+ZnJlcTE8L2NvZGU+IHNob3VsZCBiZSBpbmNsdWRlZCBpbiB0aGUgbW9kZWwuCgoxLiBBTk9WQSBvbiBhIG1vZGVsIHdoZXJlIHRoZSBjb3ZhcmlhdGVzIGFyZSB0cmVhdGVkIGFzIGZhY3RvcnMgcmVzdWx0cyBpbiBxdWFsaXRhdGl2ZWx5IHNpbWlsYXIgcmVzdWx0OiA8Y29kZT5mcmVxMTwvY29kZT4gaXMgaW1wb3J0YW50LCBidXQgPGNvZGU+aW52YXNpb248L2NvZGU+IGlzIG5vdC4KCgpgYGB7ciBOb3ZlbC5kaXZlcnNpdHksIGZpZy53aWR0aCA9IDh9CmBgYAoKCiMjIENvbW11bml0eSBkaXNzaW1pbGFyaXR5IGJlZm9yZSBhbmQgYWZ0ZXIgdGhlIHJlZ2ltZSBjaGFuZ2UuCgpIZXJlLCB3ZSBhbmFseXplIGhvdyBjb21tdW5pdHkgZGlzc2ltaWxhcml0eSBpbmRleCBjaGFuZ2VzIHdpdGggZGlzdHVyYmFuY2UgcmVnaW1lcy4KCjxiPlJlc3VsdHM6PC9iPiAKVGhlIGJlc3QgbW9kZWwgaXM6XAo8Y29kZT5iZXRhLnBhc3QgfiBmcmVxMiArIEkoZnJlcTJeMikgKyBmcmVxMSArIEkoZnJlcTFeMikgKyBmcmVxMjpmcmVxMSArIGZyZXEyOkkoZnJlcTFeMikgKyBJKGZyZXEyXjIpOmZyZXExICsgSShmcmVxMl4yKTpJKGZyZXExXjIpPC9jb2RlPi5cCgpgYGB7ciBCZXRhLmRpdmVyc2l0eSwgZmlnLndpZHRoID0gOH0KYGBgCgoKCiMgRWZmZWN0cyBvZiBkaXN0dXJiYW5jZSByZWdpbWVzIGFuZCBkaXZlcnNpdHkgb24gaW52YXNpYmlsaXR5CgpIZXJlLCB3ZSB0ZXN0IHRoZSBkaXZlcnNpdHktaW52YXNhYmlsaXR5IHJlbGF0aW9uc2hpcC4KCjxiPkh5cG90aGVzaXM6PC9iPiBQYXN0IHJlZ2ltZXMgYWZmZWN0IGludmFzaW9uIHN1Y2Nlc3MuIFNwZWNpZmljYWxseSwgbW9yZSBkaXZlcnNlIGNvbW11bml0eSB3b3VsZCBzdXBwcmVzcyBpbnZhZGVycywgc28gaW52YXNpb24gc3VjY2VzcyB3b3VsZCBmb2xsb3cgc2ltaWxhciBwYXR0ZXJuIGFzIHRoZSBkaXN0dXJiYW5jZS1kaXZlcnNpdHkgcmVsYXRpb25zaGlwLiAKCjxiPlJlc3VsdHM6PC9iPiBUaGUgYmVzdCBtb2RlbCBpczpcCjxjb2RlPmludi5yYXcgfiBmcmVxMSArIEkoZnJlcTFeMikgKyBmcmVxMiArIGFscGhhMS5wYXN0ICsgSShmcmVxMV4yKTpmcmVxMiArIGFscGhhMS5wYXN0ICogKEkoZnJlcTFeMikgKyBmcmVxMiArIEkoZnJlcTFeMik6ZnJlcTIpICsgb2Zmc2V0KGxvZyh0b3RhbC5kaWx1dGlvbikpPC9jb2RlPi5cCgoKMS4gRmlyc3QsIHdlIHRlc3RlZCB3aGV0aGVyIGJvdGggZGlzdHVyYmFuY2UgcmVnaW1lcyBhcmUgbmVjZXNzYXJ5IHRvIHVuZGVyc3RhbmQgaW52YXNpb24gb3V0Y29tZXMuIFRoZXkgYXJlOiAKPGNvZGU+YWRkMTwvY29kZT4gc2hvd3MgdGhhdCBib3RoIDxjb2RlPmZyZXExPC9jb2RlPiBhbmQgPGNvZGU+ZnJlcTI8L2NvZGU+IHNob3VsZCBiZSBpbmNsdWRlZCwgYnV0IG5vdCBib3RoIDxjb2RlPkkoZnJlcTFeMik8L2NvZGU+IGFuZCA8Y29kZT5JKGZyZXEyXjIpPC9jb2RlPi4KCjEuIFNlY29uZCwgd2UgdGVzdGVkIHdoZXRoZXIgaW5jbHVkaW5nIGRpdmVyc2l0eSBoZWxwcyBleHBsYWluIGludmFzaW9ucyBvdXRjb21lcy4gSXQgZG9lcy4gQUlDcyBzaG93IHRoYXQgaW5jbHVkaW5nIGRpdmVyc2l0eSBhbHdheXMgaW1wcm92ZXMgdGhlIG1vZGVsLgoKMS4gVGhpcmQsIHdlIHBlcmZvcm1lZCBhIHN0ZXB3aXNlIGJhY2t3YXJkIG1vZGVsIHNlbGVjdGlvbiB0byBhc3Nlc3MgaG93IGRpc3R1cmJhbmNlIHJlZ2ltZXMgYW5kIGRpdmVyc2l0eSBqb2ludGx5IGFmZmVjdCBpbnZhc2lvbiBvdXRjb21lcy4KT3VyIGFuYWx5c2VzIHNob3cgdGhhdCBib3RoIHBhc3QgYW5kIG5vdmVsIHJlZ2ltZXMgcXVhbGl0YXRpdmVseSBhbHRlciBkaXZlcnNpdHktaW52YXNpb24gcmVsYXRpb25zaGlwLgogCgoKYGBge3IgRGl2ZXJzaXR5LmludmFzaWJpbGl0eX0KYGBgCgoKCgoKIyBDb3VudGVyZmFjdHVhbHMKV2hhdCBoYXBwZW5zIGlmIHdlIGRvbid0IGNvbnNpZGVyIGRpc3R1cmJhbmNlIHJlZ2ltZSBjaGFuZ2U/CgpgYGB7ciBDb3VudGVyZmFjdHVhbHN9CmBgYAoKCgoKCgojIFN1cHBvcnRpbmcgYW5hbHlzZXMgCgojIyBBTk9WQSBvbiByZXNpZGVudCBkaXZlcnNpdHkgd2l0aCB1bm9yZGVyZWQgZmFjdG9ycy4KSGVyZSwgd2UgZnVydGhlciB0ZXN0IHdoZXRoZXIgMS4gaW52YXNpb24gZG9lc24ndCBhZmZlY3QgZGl2ZXJzaXR5LCBhbmQgMi4gcGFzdCByZWdpbWVzIGltcHJvdmUgbW9kZWwgZml0LiAKYGBge3IgTm92ZWwuZGl2ZXJzaXR5LnN1cHB9CmBgYAoKIyMgU2Vuc2l0aXZpdHkgb2YgdGhlIGRpdmVyc2l0eS1pbnZhc2liaWxpdHkgcmVsYXRpb25zaGlwcyB0byBtb2RlbCBzdHJ1Y3R1cmUuCkhlcmUsIHdlIGFzc2VzcyBob3cgb3VyIHJlc3VsdHMgb24gZGl2ZXJzaXR5LWludmFzaWJpbGl0eSByZWxhdGlvbnNoaXBzIGNoYW5nZSB3aGVuIHdlIGFuYWx5emUgYSBzaW1wbGVyIG1vZGVsLgpgYGB7ciBJbnZhZGVycy5zdXBwfQpgYGAKCiMjIEVmZmVjdHMgb2YgcmVnaW1lIGNoYW5nZXMgYWxvbmUgb24gaW52YXNpYmlsaXR5LgpIZXJlLCB3ZSB0ZXN0IHdoZXRoZXIgcmVnaW1lIGNoYW5nZSBhbG9uZSBjYW4gaGVscCBleHBsYWluIHRoZSBkaWZmZXJlbmNlcyBpbiBpbnZhc2lvbiBvdXRjb21lcy4KYGBge3IgSW52YWRlcnMucmN9CmBgYAo=