Last updated: 2020-10-11
Checks: 7 0
Knit directory: genes-to-foodweb-stability/
This reproducible R Markdown analysis was created with workflowr (version 1.6.0). 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 version displayed above was the version of the Git repository at the time these results were generated.
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
Untracked files:
Untracked: output/community-persistence-keystone.RData
Untracked: output/critical-transitions-keystone.RData
Untracked: output/structural-stability-keystone.RData
Unstaged changes:
Modified: output/plant-growth-no-insects.RData
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 R Markdown and HTML files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view them.
| File | Version | Author | Date | Message |
|---|---|---|---|---|
| Rmd | 761af40 | mabarbour | 2020-10-11 | Refocus analysis on keystone gene result. |
# Load and manage data
df <- read_csv("data/arabidopsis_clean_df.csv") %>%
# renaming for brevity
rename(cage = Cage,
com = Composition,
week = Week,
temp = Temperature,
rich = Richness) %>%
mutate(cage = as.character(cage),
fweek = factor(ifelse(week < 10, paste("0", week, sep=""), week)),
temp = ifelse(temp=="20 C", 0, 1)) %>%
arrange(cage, week)
# create data for multi-state survival analysis
state_df <- df %>%
# counter information is not relevant (because it is the same), so we summarise across it
group_by(cage, fweek, week, temp, rich, Col, gsm1, AOP2, AOP2.gsoh, com) %>%
summarise_at(vars(BRBR_Survival, LYER_Survival, Mummy_Ptoids_Survival), list(mean)) %>%
ungroup() %>%
# create possible food-web states
mutate(BRBR = ifelse(BRBR_Survival == 1, "BRBR", ifelse(BRBR_Survival == 0, "0", NA)),
LYER = ifelse(LYER_Survival == 1, "LYER", ifelse(LYER_Survival == 0, "0", NA)),
Ptoid = ifelse(Mummy_Ptoids_Survival == 1, "Ptoid", ifelse(Mummy_Ptoids_Survival == 0, "0", NA))) %>%
mutate(state = paste(BRBR, LYER, Ptoid, sep = "-"),
cage = as.character(cage)) %>%
# these variables are no longer needed
select(-BRBR, -LYER, -Ptoid) %>%
# remove all instances where all species have been labelled extinct for more than 1 week "NA-NA-NA".
# also, only keep observations after 2 weeks when we added the parasitoid (full community)
filter(state != "NA-NA-NA", week > 2) %>%
mutate(week_since = week - 2)
# replace NA with zeros for state variable
state_df$state <- gsub("NA","0", state_df$state)
# everything appears in order
arrange(state_df, cage, week_since) %>% select(cage, week_since, BRBR_Survival, LYER_Survival, Mummy_Ptoids_Survival, state) # A tibble: 710 x 6
cage week_since BRBR_Survival LYER_Survival Mummy_Ptoids_Surviv… state
<chr> <dbl> <dbl> <dbl> <dbl> <chr>
1 1 1 1 1 1 BRBR-LYER-…
2 1 2 1 1 1 BRBR-LYER-…
3 1 3 1 1 1 BRBR-LYER-…
4 1 4 1 1 1 BRBR-LYER-…
5 1 5 0 0 1 0-0-Ptoid
6 1 6 NA NA 0 0-0-0
7 10 1 1 1 1 BRBR-LYER-…
8 10 2 1 1 1 BRBR-LYER-…
9 10 3 1 1 1 BRBR-LYER-…
10 10 4 1 1 1 BRBR-LYER-…
# … with 700 more rows# source in useful functions for analyses
source('code/glm-ftest.R') # ANOVA GLM
delta_prob = function(x) round(exp(x) - 1,2) # change in probability of transtion from discrete-time survival analysisLet’s look at the possible state transitions for all cages and time points
# state transitions for all cages and time points
msm::statetable.msm(state = state, subject = cage, data = state_df) to
from 0-0-0 0-0-Ptoid 0-LYER-0 0-LYER-Ptoid BRBR-LYER-Ptoid
0-0-Ptoid 22 7 0 0 0
0-LYER-0 1 0 131 0 0
0-LYER-Ptoid 5 12 26 242 0
BRBR-LYER-Ptoid 0 10 0 50 144The most common transition from the full community BRBR-LYER-Ptoid is to lose BRBR 0-LYER-Ptoid and less common to collapse by losing both aphids 0-0-Ptoid. The 0-0-Ptoid state never lasts long, since the parasitoid no longer has any resources, and collapses to 0-0-0.
From the 0-LYER-Ptoid community, losing the parasitoid 0-LYER-0 is the most common transition, but losing LYER 0-0-Ptoid is also common. Note that transitions to complete collapse 0-0-0 were more than likely to occur by losing LYER then the parasitoid, rather than losing the parasitoid and then LYER. This is because collapsing from the 0-LYER-0 state is rare (1 / (131 + 1) = 0.7575758% chance per week).
These transition patterns suggest the following:
0-0-Ptoid state always collapses to 0-0-0, I’m going to convert the 0-0-Ptoid state into 0-0-0. This will reduce the number of state transitions I need to model, but it won’t sacrifice modeling the biological process.0-LYER-0 to 0-0-0 transition is so rare, it does not make sense to model this transition with any covariates.Below, we reproduce the analysis of deviance for discrete-time survival models presented in Tables S2–S3 in the Supplementary Material.
You’ll notice that we code the effect of the null AOP2\(-\) allele as I(Col + gsm1), as it represents the average effect of genotypes with this allele. Similarly, we code the effect of the functional AOP2\(+\) allele as I(AOP2 + AOP2.gsoh).
# filter data
full_transit_df <- state_df %>%
# removes all transitions after the first one, this focuses
# on transition from the full community
na.omit() # check for weekly variation in critical transitions.
full_transit_df %>%
mutate(CT_from_full = ifelse(state %in% c("0-LYER-Ptoid","0-0-Ptoid"), 1, 0)) %>%
group_by(fweek) %>%
summarise_at("CT_from_full", mean)# A tibble: 7 x 2
fweek CT_from_full
<fct> <dbl>
1 03 0
2 04 0.0833
3 05 0.127
4 06 0.333
5 07 0.75
6 08 0.875
7 09 1 We only observe variation between weeks 4 and 8, so we restrict the data to those weeks for the analysis. It does not make sense to include weeks where there is no variation. If we did, then the factor fweek would cause predict all of the variation in the week with no variation.
# filter data again
full_transit_df <- filter(full_transit_df, week %in% c(4:8))# adjust rich and temp so coefficients = +1 genotype and +1 C, respectively
full_transit_df <- full_transit_df %>%
mutate(rich = rich - 1,
temp = temp * 3)Now we’ve put the coefficients on a similar and more intuitive scale with respect to each other.
# AOP2- allele test
glm.ftest.v2(
model = glm(data = full_transit_df,
family = quasibinomial(link = "cloglog"),
formula = terms(state %in% c("0-LYER-Ptoid","0-0-Ptoid") ~
fweek + temp + I(AOP2 + AOP2.gsoh) + I(Col + gsm1) + com + temp:com,
keep.order = T)),
test.formula = list(
c("fweek","Residuals"),
c("temp","temp:com"),
c("I(Col + gsm1)","com"))
)[[3]] %>%
select(Source = treatment,
`df (Source)` = num_df,
`df (Error)` = den_df,
Deviance = deviance,
`Mean Deviance` = mean_deviance,
F = F, P = P, Error = error) %>%
kable(., caption = "Analysis of Deviance for any critical transition from the initial food web. Appropriate test for AOP2- allele", booktabs = T) %>%
kable_styling(latex_options = c("striped", "hold_position"))Warning: glm.fit: algorithm did not converge| Source | df (Source) | df (Error) | Deviance | Mean Deviance | F | P | Error |
|---|---|---|---|---|---|---|---|
| fweek | 4 | 177 | 65.23 | 16.31 | 20.854 | <0.001 | Residuals |
| temp | 1 | 10 | 14.81 | 14.81 | 11.720 | 0.007 | temp:com |
| I(Col + gsm1) | 1 | 8 | 3.79 | 3.79 | 3.106 | 0.116 | com |
# AOP2+ allele test
glm.ftest.v2(
model = glm(data = full_transit_df,
family = quasibinomial(link = "cloglog"),
formula = terms(state %in% c("0-LYER-Ptoid","0-0-Ptoid") ~
fweek + temp + I(Col + gsm1) + I(AOP2 + AOP2.gsoh) + com + temp:com,
keep.order = T)),
test.formula = list(
c("fweek","Residuals"),
c("temp","temp:com"),
c("I(AOP2 + AOP2.gsoh)","com"))
)[[3]] %>%
select(Source = treatment,
`df (Source)` = num_df,
`df (Error)` = den_df,
Deviance = deviance,
`Mean Deviance` = mean_deviance,
F = F, P = P, Error = error) %>%
kable(., caption = "Analysis of Deviance for any critical transition from the initial food web. Appropriate test for AOP2+ allele", booktabs = T) %>%
kable_styling(latex_options = c("striped", "hold_position"))Warning: glm.fit: algorithm did not converge| Source | df (Source) | df (Error) | Deviance | Mean Deviance | F | P | Error |
|---|---|---|---|---|---|---|---|
| fweek | 4 | 177 | 65.23 | 16.31 | 20.854 | <0.001 | Residuals |
| temp | 1 | 10 | 14.81 | 14.81 | 11.720 | 0.007 | temp:com |
| I(AOP2 + AOP2.gsoh) | 1 | 8 | 0.09 | 0.09 | 0.070 | 0.798 | com |
The model didn’t converge…Let’s try the standard link = "logit".
# AOP2- allele test
glm.ftest.v2(
model = glm(data = full_transit_df,
family = quasibinomial(link = "logit"),
formula = terms(state %in% c("0-LYER-Ptoid","0-0-Ptoid") ~
fweek + temp + I(AOP2 + AOP2.gsoh) + I(Col + gsm1) + com + temp:com,
keep.order = T)),
test.formula = list(
c("fweek","Residuals"),
c("temp","temp:com"),
c("I(Col + gsm1)","com"))
)[[3]] %>%
select(Source = treatment,
`df (Source)` = num_df,
`df (Error)` = den_df,
Deviance = deviance,
`Mean Deviance` = mean_deviance,
F = F, P = P, Error = error) Source df (Source) df (Error) Deviance Mean Deviance F P
1 fweek 4 177 65.23 16.31 21.197 <0.001
2 temp 1 10 21.87 21.87 19.673 0.001
3 I(Col + gsm1) 1 8 1.78 1.78 1.680 0.231
Error
1 Residuals
2 temp:com
3 com# AOP2+ allele test
glm.ftest.v2(
model = glm(data = full_transit_df,
family = quasibinomial(link = "logit"),
formula = terms(state %in% c("0-LYER-Ptoid","0-0-Ptoid") ~
fweek + temp + I(Col + gsm1) + I(AOP2 + AOP2.gsoh) + com + temp:com,
keep.order = T)),
test.formula = list(
c("fweek","Residuals"),
c("temp","temp:com"),
c("I(AOP2 + AOP2.gsoh)","com"))
)[[3]] %>%
select(Source = treatment,
`df (Source)` = num_df,
`df (Error)` = den_df,
Deviance = deviance,
`Mean Deviance` = mean_deviance,
F = F, P = P, Error = error) Source df (Source) df (Error) Deviance Mean Deviance F
1 fweek 4 177 65.23 16.31 21.197
2 temp 1 10 21.87 21.87 19.673
3 I(AOP2 + AOP2.gsoh) 1 8 0.01 0.01 0.008
P Error
1 <0.001 Residuals
2 0.001 temp:com
3 0.931 comEverything converges now and the results are qualitatively the same as when link = "cloglog". To maintain consistency with the other models, I report the results from link = "cloglog".
# check for weekly variation in critical transitions.
full_transit_df %>%
mutate(CT_to_0LYERPtoid = ifelse(state %in% "0-LYER-Ptoid", 1, 0)) %>%
group_by(week) %>%
summarise_at("CT_to_0LYERPtoid", mean)# A tibble: 5 x 2
week CT_to_0LYERPtoid
<dbl> <dbl>
1 4 0.0833
2 5 0.127
3 6 0.333
4 7 0.5
5 8 0.75 We keep the full_transit_df since we observe variation in all weeks.
# AOP2- allele
glm.ftest.v2(
model = glm(data = full_transit_df,
family = quasibinomial(link = "cloglog"),
formula = terms(state %in% c("0-LYER-Ptoid") ~
fweek + temp + I(AOP2 + AOP2.gsoh) + I(Col + gsm1) + com + temp:com,
keep.order = T)),
test.formula = list(
c("fweek","Residuals"),
c("temp","temp:com"),
#c("I(AOP2 + AOP2.gsoh)","com"),
c("I(Col + gsm1)","com"))
)[[3]] %>%
select(Source = treatment,
`df (Source)` = num_df,
`df (Error)` = den_df,
Deviance = deviance,
`Mean Deviance` = mean_deviance,
F = F, P = P, Error = error) %>%
kable(., caption = "Analysis of Deviance for critical transition from initial food web to food chain with dominant aphid and parasitoid. Appropriate test for AOP2- allele.", booktabs = T) %>%
kable_styling(latex_options = c("striped", "hold_position"))| Source | df (Source) | df (Error) | Deviance | Mean Deviance | F | P | Error |
|---|---|---|---|---|---|---|---|
| fweek | 4 | 177 | 35.83 | 8.96 | 10.916 | <0.001 | Residuals |
| temp | 1 | 10 | 20.19 | 20.19 | 12.163 | 0.006 | temp:com |
| I(Col + gsm1) | 1 | 8 | 0.41 | 0.41 | 0.394 | 0.548 | com |
# AOP2+ allele
glm.ftest.v2(
model = glm(data = full_transit_df,
family = quasibinomial(link = "cloglog"),
formula = terms(state %in% c("0-LYER-Ptoid") ~
fweek + temp + I(Col + gsm1) + I(AOP2 + AOP2.gsoh) + com + temp:com,
keep.order = T)),
test.formula = list(
c("fweek","Residuals"),
c("temp","temp:com"),
c("I(AOP2 + AOP2.gsoh)","com"))
)[[3]] %>%
select(Source = treatment,
`df (Source)` = num_df,
`df (Error)` = den_df,
Deviance = deviance,
`Mean Deviance` = mean_deviance,
F = F, P = P, Error = error) %>%
kable(., caption = "Analysis of Deviance for critical transition from initial food web to food chain with dominant aphid and parasitoid. Appropriate test for AOP2+ allele.", booktabs = T) %>%
kable_styling(latex_options = c("striped", "hold_position"))| Source | df (Source) | df (Error) | Deviance | Mean Deviance | F | P | Error |
|---|---|---|---|---|---|---|---|
| fweek | 4 | 177 | 35.83 | 8.96 | 10.916 | <0.001 | Residuals |
| temp | 1 | 10 | 20.19 | 20.19 | 12.163 | 0.006 | temp:com |
| I(AOP2 + AOP2.gsoh) | 1 | 8 | 0.07 | 0.07 | 0.070 | 0.799 | com |
# calculate effects sizes for AOP2- and AOP2+ alleles, Fig. 3 of main text
conf_int(
glm(data = full_transit_df,
family = quasibinomial(link = "cloglog"),
formula = state %in% c("0-LYER-Ptoid") ~
fweek + temp + I(AOP2 + AOP2.gsoh) + I(Col + gsm1)),
vcov = "CR2",
test = "naive-t",
cluster = full_transit_df$com,
coefs = c("I(AOP2 + AOP2.gsoh)","I(Col + gsm1)")
) %>%
data.frame() %>%
rownames_to_column(var = "term") %>%
group_by(term) %>%
summarise_at(vars(beta), list(delta_prob = delta_prob))# A tibble: 2 x 2
term delta_prob
<chr> <dbl>
1 I(AOP2 + AOP2.gsoh) -0.06
2 I(Col + gsm1) -0.13# check for weekly variation in critical transitions.
full_transit_df %>%
mutate(CT_to_00Ptoid = ifelse(state %in% "0-0-Ptoid", 1, 0)) %>%
group_by(week) %>%
summarise_at("CT_to_00Ptoid", mean)# A tibble: 5 x 2
week CT_to_00Ptoid
<dbl> <dbl>
1 4 0
2 5 0
3 6 0
4 7 0.25
5 8 0.125We only see variation in weeks 7 and 8, so we restrict the data to those weeks.
# check frequency of critical transitions
with(filter(full_transit_df, week %in% 7:8), table(state %in% c("0-0-Ptoid")))
FALSE TRUE
31 9 And since there were only 9 total transitions to this state, we don’t attempt to fit the interaction term temp:rich.
# AOP2- allele
glm.ftest.v2(
model = glm(data = filter(full_transit_df, week %in% 7:8),
family = quasibinomial(link = "cloglog"),
formula = terms(state %in% c("0-0-Ptoid") ~
fweek + temp + I(AOP2 + AOP2.gsoh) + I(Col + gsm1) + com + temp:com,
keep.order = T)),
test.formula = list(
c("fweek","Residuals"),
c("temp","temp:com"),
#c("I(AOP2 + AOP2.gsoh)","com"),
c("I(Col + gsm1)","com"))
)[[3]] %>%
select(Source = treatment,
`df (Source)` = num_df,
`df (Error)` = den_df,
Deviance = deviance,
`Mean Deviance` = mean_deviance,
F = F, P = P, Error = error) %>%
kable(., caption = "Analysis of deviance for critical transition from initial food web to extinction of all species. Appropriate test for AOP2- allele.", booktabs = T) %>%
kable_styling(latex_options = c("striped", "hold_position"))| Source | df (Source) | df (Error) | Deviance | Mean Deviance | F | P | Error |
|---|---|---|---|---|---|---|---|
| fweek | 1 | 21 | 0.64 | 0.64 | 0.631 | 0.436 | Residuals |
| temp | 1 | 6 | 0.93 | 0.93 | 1.056 | 0.344 | temp:com |
| I(Col + gsm1) | 1 | 8 | 1.83 | 1.83 | 1.159 | 0.313 | com |
# AOP2+ allele
glm.ftest.v2(
model = glm(data = filter(full_transit_df, week %in% 7:8),
family = quasibinomial(link = "cloglog"),
formula = terms(state %in% c("0-0-Ptoid") ~
fweek + temp + I(Col + gsm1) + I(AOP2 + AOP2.gsoh) + com + temp:com,
keep.order = T)),
test.formula = list(
c("fweek","Residuals"),
c("temp","temp:com"),
c("I(AOP2 + AOP2.gsoh)","com"))
)[[3]] %>%
select(Source = treatment,
`df (Source)` = num_df,
`df (Error)` = den_df,
Deviance = deviance,
`Mean Deviance` = mean_deviance,
F = F, P = P, Error = error) %>%
kable(., caption = "Analysis of deviance for critical transition from initial food web to extinction of all species. Appropriate test for AOP2+ allele.", booktabs = T) %>%
kable_styling(latex_options = c("striped", "hold_position"))| Source | df (Source) | df (Error) | Deviance | Mean Deviance | F | P | Error |
|---|---|---|---|---|---|---|---|
| fweek | 1 | 21 | 0.64 | 0.64 | 0.631 | 0.436 | Residuals |
| temp | 1 | 6 | 0.93 | 0.93 | 1.056 | 0.344 | temp:com |
| I(AOP2 + AOP2.gsoh) | 1 | 8 | 0.08 | 0.08 | 0.048 | 0.831 | com |
# calculate effects sizes for AOP2- and AOP2, Fig. 3 of main text
conf_int(
glm(data = filter(full_transit_df, week %in% 7:8),
family = quasibinomial(link = "cloglog"),
formula = state %in% c("0-0-Ptoid") ~
fweek + temp + I(AOP2 + AOP2.gsoh) + I(Col + gsm1)),
vcov = "CR2",
test = "naive-t",
cluster = filter(full_transit_df, week %in% 7:8)$com,
coefs = c("I(AOP2 + AOP2.gsoh)","I(Col + gsm1)")
) %>%
data.frame() %>%
rownames_to_column(var = "term") %>%
group_by(term) %>%
summarise_at(vars(beta), list(delta_prob = delta_prob))# A tibble: 2 x 2
term delta_prob
<chr> <dbl>
1 I(AOP2 + AOP2.gsoh) 0.15
2 I(Col + gsm1) -0.5 ## Organize data for analysis
# LYER-Ptoid cages at least at one time point
LP_cages <- unique(filter(state_df, state == "0-LYER-Ptoid")$cage)
# should be 50 cages, and it is
length(LP_cages)[1] 50# filter and manage data
LP_transit_df <- state_df %>%
filter(cage %in% LP_cages, state != "BRBR-LYER-Ptoid") %>%
# omit BRBR from consideration
select(-BRBR_Survival) %>%
# omit rows where we already know either LYER or Ptoid went extinct
na.omit() %>%
mutate(# assume these two states are the same, i.e. to get to 0-0-0, had to go through 0-0-Ptoid
state_adj = ifelse(state %in% c("0-0-Ptoid","0-0-0"), "0-0-Ptoid", state))
# confirm levels
unique(LP_transit_df$state_adj) #fstate)[1] "0-LYER-Ptoid" "0-LYER-0" "0-0-Ptoid" # check for weekly variation in critical transitions
LP_transit_df %>%
mutate(CT_from_FoodChain = ifelse(state_adj %in% c("0-0-Ptoid","0-LYER-0"), 1, 0)) %>%
group_by(fweek) %>%
summarise_at("CT_from_FoodChain", mean)# A tibble: 14 x 2
fweek CT_from_FoodChain
<fct> <dbl>
1 04 0
2 05 0
3 06 0
4 07 0.136
5 08 0.205
6 09 0.143
7 10 0.1
8 11 0.0370
9 12 0.0769
10 13 0.0833
11 14 0.136
12 15 0.368
13 16 0.417
14 17 0 We only observe variation between weeks 7 and 16, so we restrict the data to those weeks for the analysis.
LP_transit_df <- filter(LP_transit_df, week %in% c(7:16))# adjust rich and temp coefficients, so +1 genotype corresponds to +1 C
LP_transit_df <- LP_transit_df %>%
mutate(rich = rich-1,
temp = temp*3)with(LP_transit_df, table(state_adj))state_adj
0-0-Ptoid 0-LYER-0 0-LYER-Ptoid
17 26 240 # AOP2- allele test
glm.ftest.v2(
model = glm(data = LP_transit_df,
family = quasibinomial(link = "cloglog"),
formula = terms(state_adj %in% c("0-LYER-0","0-0-Ptoid") ~
fweek + temp + I(AOP2 + AOP2.gsoh) + I(Col + gsm1) + com + temp:com,
keep.order = T)),
test.formula = list(
c("fweek","Residuals"),
c("temp","temp:com"),
#c("I(AOP2 + AOP2.gsoh)","com"),
c("I(Col + gsm1)","com"))
)[[3]] %>%
select(Source = treatment,
`df (Source)` = num_df,
`df (Error)` = den_df,
Deviance = deviance,
`Mean Deviance` = mean_deviance,
F = F, P = P, Error = error) %>%
kable(., caption = "Analysis of deviance for any critical transition from food chain. Appropriate test for AOP2- allele.", booktabs = T) %>%
kable_styling(latex_options = c("striped", "hold_position"))| Source | df (Source) | df (Error) | Deviance | Mean Deviance | F | P | Error |
|---|---|---|---|---|---|---|---|
| fweek | 9 | 252 | 18.05 | 2.01 | 2.544 | 0.008 | Residuals |
| temp | 1 | 10 | 0.01 | 0.01 | 0.005 | 0.945 | temp:com |
| I(Col + gsm1) | 1 | 8 | 5.98 | 5.98 | 11.636 | 0.009 | com |
# AOP2+ allele test
glm.ftest.v2(
model = glm(data = LP_transit_df,
family = quasibinomial(link = "cloglog"),
formula = terms(state_adj %in% c("0-LYER-0","0-0-Ptoid") ~
fweek + temp + I(Col + gsm1) + I(AOP2 + AOP2.gsoh) + com + temp:com,
keep.order = T)),
test.formula = list(
c("fweek","Residuals"),
c("temp","temp:com"),
c("I(AOP2 + AOP2.gsoh)","com"))
)[[3]] %>%
select(Source = treatment,
`df (Source)` = num_df,
`df (Error)` = den_df,
Deviance = deviance,
`Mean Deviance` = mean_deviance,
F = F, P = P, Error = error) %>%
kable(., caption = "Analysis of deviance for critical transition from food chain. Appropriate test for AOP2+ allele.", booktabs = T) %>%
kable_styling(latex_options = c("striped", "hold_position"))| Source | df (Source) | df (Error) | Deviance | Mean Deviance | F | P | Error |
|---|---|---|---|---|---|---|---|
| fweek | 9 | 252 | 18.05 | 2.01 | 2.544 | 0.008 | Residuals |
| temp | 1 | 10 | 0.01 | 0.01 | 0.005 | 0.945 | temp:com |
| I(AOP2 + AOP2.gsoh) | 1 | 8 | 0.99 | 0.99 | 1.927 | 0.203 | com |
# check for weekly variation in critical transitions
LP_transit_df %>%
mutate(CT_FoodChain_to_00Ptoid = ifelse(state_adj == c("0-0-Ptoid"), 1, 0)) %>%
group_by(week) %>%
summarise_at("CT_FoodChain_to_00Ptoid", mean)# A tibble: 10 x 2
week CT_FoodChain_to_00Ptoid
<dbl> <dbl>
1 7 0.114
2 8 0.136
3 9 0
4 10 0
5 11 0
6 12 0
7 13 0
8 14 0
9 15 0.211
10 16 0.167Only observe variation in weeks 7 and 8 and 15 and 16.
# AOP2- allele
glm.ftest.v2(
model = glm(data = filter(LP_transit_df, week %in% c(7,8,15,16)),
family = quasibinomial(link = "cloglog"),
formula = terms(state_adj %in% c("0-0-Ptoid") ~
fweek + temp + I(AOP2 + AOP2.gsoh) + I(Col + gsm1) + com + temp:com,
keep.order = T)),
test.formula = list(
c("fweek","Residuals"),
c("temp","temp:com"),
#c("I(AOP2 + AOP2.gsoh)","com"),
c("I(Col + gsm1)","com"))
)[[3]] %>%
select(Source = treatment,
`df (Source)` = num_df,
`df (Error)` = den_df,
Deviance = deviance,
`Mean Deviance` = mean_deviance,
F = F, P = P, Error = error) %>%
kable(., caption = "Analysis of deviance for critical transition from food chain to a complete collapse (loss of both aphid and parasitoid). Appropriate test of AOP2- allele.", booktabs = T) %>%
kable_styling(latex_options = c("striped", "hold_position"))| Source | df (Source) | df (Error) | Deviance | Mean Deviance | F | P | Error |
|---|---|---|---|---|---|---|---|
| fweek | 3 | 94 | 1.03 | 0.34 | 0.469 | 0.705 | Residuals |
| temp | 1 | 10 | 0.51 | 0.51 | 0.391 | 0.546 | temp:com |
| I(Col + gsm1) | 1 | 8 | 6.85 | 6.85 | 7.604 | 0.025 | com |
# AOP2+ allele
glm.ftest.v2(
model = glm(data = filter(LP_transit_df, week %in% c(7,8,15,16)),
family = quasibinomial(link = "cloglog"),
formula = terms(state_adj %in% c("0-0-Ptoid") ~
fweek + temp + I(Col + gsm1) + I(AOP2 + AOP2.gsoh) + com + temp:com,
keep.order = T)),
test.formula = list(
c("fweek","Residuals"),
c("temp","temp:com"),
c("I(AOP2 + AOP2.gsoh)","com"))
)[[3]] %>%
select(Source = treatment,
`df (Source)` = num_df,
`df (Error)` = den_df,
Deviance = deviance,
`Mean Deviance` = mean_deviance,
F = F, P = P, Error = error) %>%
kable(., caption = "Analysis of deviance for critical transition from food chain to a complete collapse (loss of both aphid and parasitoid). Appropriate test of AOP2+ allele.", booktabs = T) %>%
kable_styling(latex_options = c("striped", "hold_position"))| Source | df (Source) | df (Error) | Deviance | Mean Deviance | F | P | Error |
|---|---|---|---|---|---|---|---|
| fweek | 3 | 94 | 1.03 | 0.34 | 0.469 | 0.705 | Residuals |
| temp | 1 | 10 | 0.51 | 0.51 | 0.391 | 0.546 | temp:com |
| I(AOP2 + AOP2.gsoh) | 1 | 8 | 1.00 | 1.00 | 1.115 | 0.322 | com |
# keystone allele effects sizes in Fig. 3 of main text
conf_int(
glm(data = filter(LP_transit_df, week %in% c(7,8,15,16)),
family = quasibinomial(link = "cloglog"),
formula = state_adj %in% c("0-0-Ptoid") ~
fweek + temp + I(AOP2 + AOP2.gsoh) + I(Col + gsm1)),
vcov = "CR2",
test = "naive-t",
cluster = filter(LP_transit_df, week %in% c(7,8,15,16))$com,
coefs = c("I(AOP2 + AOP2.gsoh)","I(Col + gsm1)")
) %>%
data.frame() %>%
rownames_to_column(var = "term") %>%
group_by(term) %>%
summarise_at(vars(beta), list(delta_prob = delta_prob))# A tibble: 2 x 2
term delta_prob
<chr> <dbl>
1 I(AOP2 + AOP2.gsoh) -0.35
2 I(Col + gsm1) -0.66LP_transit_df %>%
mutate(CT_FoodChain_to_0LYER0 = ifelse(state_adj == c("0-LYER-0"), 1, 0)) %>%
group_by(week) %>%
summarise_at("CT_FoodChain_to_0LYER0", mean)# A tibble: 10 x 2
week CT_FoodChain_to_0LYER0
<dbl> <dbl>
1 7 0.0227
2 8 0.0682
3 9 0.143
4 10 0.1
5 11 0.0370
6 12 0.0769
7 13 0.0833
8 14 0.136
9 15 0.158
10 16 0.25 We observe variation in all of the same weeks as LP_transit_df, so we continue with it for the analysis.
# AOP2- allele
glm.ftest.v2(
model = glm(data = LP_transit_df,
family = quasibinomial(link = "cloglog"),
formula = terms(state_adj %in% c("0-LYER-0") ~
fweek + temp + I(AOP2 + AOP2.gsoh) + I(Col + gsm1) + com + temp:com,
keep.order = T)),
test.formula = list(
c("fweek","Residuals"),
c("temp","temp:com"),
#c("I(AOP2 + AOP2.gsoh)","com"),
c("I(Col + gsm1)","com"))
)[[3]] %>%
select(Source = treatment,
`df (Source)` = num_df,
`df (Error)` = den_df,
Deviance = deviance,
`Mean Deviance` = mean_deviance,
F = F, P = P, Error = error) %>%
# table of results
kable(., caption = "Analysis of deviance for critical transition from food chain to an aphid only food web. Appropriate test of AOP2- allele.", booktabs = T) %>%
kable_styling(latex_options = c("striped", "hold_position"))| Source | df (Source) | df (Error) | Deviance | Mean Deviance | F | P | Error |
|---|---|---|---|---|---|---|---|
| fweek | 9 | 252 | 9.99 | 1.11 | 2.094 | 0.031 | Residuals |
| temp | 1 | 10 | 0.36 | 0.36 | 0.245 | 0.631 | temp:com |
| I(Col + gsm1) | 1 | 8 | 0.66 | 0.66 | 0.367 | 0.561 | com |
# AOP2+ allele
glm.ftest.v2(
model = glm(data = LP_transit_df,
family = quasibinomial(link = "cloglog"),
formula = terms(state_adj %in% c("0-LYER-0") ~
fweek + temp + I(Col + gsm1) + I(AOP2 + AOP2.gsoh) + com + temp:com,
keep.order = T)),
test.formula = list(
c("fweek","Residuals"),
c("temp","temp:com"),
c("I(AOP2 + AOP2.gsoh)","com"))
)[[3]] %>%
select(Source = treatment,
`df (Source)` = num_df,
`df (Error)` = den_df,
Deviance = deviance,
`Mean Deviance` = mean_deviance,
F = F, P = P, Error = error) %>%
# table of results
kable(., caption = "Analysis of deviance for critical transition from food chain to an aphid only food web. Appropriate test of AOP2+ allele.", booktabs = T) %>%
kable_styling(latex_options = c("striped", "hold_position"))| Source | df (Source) | df (Error) | Deviance | Mean Deviance | F | P | Error |
|---|---|---|---|---|---|---|---|
| fweek | 9 | 252 | 9.99 | 1.11 | 2.094 | 0.031 | Residuals |
| temp | 1 | 10 | 0.36 | 0.36 | 0.245 | 0.631 | temp:com |
| I(AOP2 + AOP2.gsoh) | 1 | 8 | 0.25 | 0.25 | 0.141 | 0.717 | com |
# keystone allele effect sizes in Fig. 3
conf_int(
glm(data = LP_transit_df,
family = quasibinomial(link = "cloglog"),
formula = state_adj %in% c("0-LYER-0") ~
fweek + temp + I(AOP2 + AOP2.gsoh) + I(Col + gsm1)),
vcov = "CR2",
test = "naive-t",
cluster = LP_transit_df$com,
coefs = c("I(AOP2 + AOP2.gsoh)","I(Col + gsm1)")
) %>%
data.frame() %>%
rownames_to_column(var = "term") %>%
group_by(term) %>%
summarise_at(vars(beta), list(delta_prob = delta_prob))# A tibble: 2 x 2
term delta_prob
<chr> <dbl>
1 I(AOP2 + AOP2.gsoh) -0.15
2 I(Col + gsm1) -0.22Write out an .RData file to use for creating the Supplementary Material Results.
# save.image(file = "output/critical-transitions-keystone.RData")
sessionInfo()R version 3.6.3 (2020-02-29)
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] kableExtra_1.1.0 clubSandwich_0.3.5 msm_1.6.8 cowplot_1.0.0
[5] forcats_0.4.0 stringr_1.4.0 dplyr_0.8.3 purrr_0.3.3
[9] readr_1.3.1 tidyr_1.0.2 tibble_2.1.3 ggplot2_3.2.1
[13] tidyverse_1.3.0
loaded via a namespace (and not attached):
[1] Rcpp_1.0.2 lubridate_1.7.4 mvtnorm_1.0-10 lattice_0.20-38
[5] zoo_1.8-6 utf8_1.1.4 assertthat_0.2.1 rprojroot_1.3-2
[9] digest_0.6.20 R6_2.4.0 cellranger_1.1.0 backports_1.1.4
[13] reprex_0.3.0 evaluate_0.14 highr_0.8 httr_1.4.1
[17] pillar_1.4.2 rlang_0.4.4 lazyeval_0.2.2 readxl_1.3.1
[21] rstudioapi_0.10 whisker_0.3-2 Matrix_1.2-17 rmarkdown_2.0
[25] splines_3.6.3 webshot_0.5.1 munsell_0.5.0 broom_0.5.2
[29] compiler_3.6.3 httpuv_1.5.1 modelr_0.1.5 xfun_0.9
[33] pkgconfig_2.0.2 htmltools_0.3.6 tidyselect_0.2.5 expm_0.999-4
[37] workflowr_1.6.0 fansi_0.4.0 viridisLite_0.3.0 crayon_1.3.4
[41] dbplyr_1.4.2 withr_2.1.2 later_1.0.0 grid_3.6.3
[45] nlme_3.1-140 jsonlite_1.6 gtable_0.3.0 lifecycle_0.1.0
[49] DBI_1.0.0 git2r_0.26.1 magrittr_1.5 scales_1.0.0
[53] cli_1.1.0 stringi_1.4.3 fs_1.3.1 promises_1.0.1
[57] xml2_1.2.2 generics_0.0.2 vctrs_0.2.2 sandwich_2.5-1
[61] tools_3.6.3 glue_1.3.1 hms_0.5.3 survival_3.1-8
[65] yaml_2.2.0 colorspace_1.4-1 rvest_0.3.5 knitr_1.26
[69] haven_2.2.0