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
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/mark-recapture.Rmd) and HTML (docs/mark-recapture.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 |
|---|---|---|---|---|
| Rmd | 18d1722 | mabarbour | 2022-01-24 | Add mark-recapture and omnibus tests, fix data point error, update analyses, and edit readme |
# load and manage data
df <- read_csv("data/InsectAbundanceSurvival.csv") %>%
# renaming for brevity
rename(cage = Cage,
com = Composition,
week = Week,
temp = Temperature,
rich = Richness) %>%
mutate(cage = as.character(cage),
temp = ifelse(temp=="20 C", 0, 3), # put on scale such that coef corresponds to 1 C increase
temp_com = paste(temp, com)) %>% # create indicator for genetic composition nested within temperature for random effect
arrange(cage, week)
# create data for survival analysis
survival_df <- df %>%
# counter information is not relevant for survival (because it is the same), so we summarise across it
group_by(cage, week, temp, rich, Col, gsm1, AOP2, AOP2.gsoh, com, temp_com) %>%
summarise_at(vars(BRBR, LYER, Mummy_Ptoids, BRBR_Survival, LYER_Survival, Mummy_Ptoids_Survival), list(mean)) %>%
ungroup() %>%
mutate(week_since = week - 2) %>% # week since the full community was added (at week 3)
filter(week_since > 0) %>%
mutate(interval_start = week_since - 1,
interval_stop = week_since) %>%
# contrasts for allelic effects at each gene
mutate(aop2_vs_AOP2 = Col + gsm1 - AOP2 - AOP2.gsoh,
mam1_vs_MAM1 = gsm1 - Col,
gsoh_vs_GSOH = AOP2.gsoh - AOP2)
mark_recap_wide_df <- survival_df %>%
# include apparent extinctions. Note that at week 3, parasitoids were still being added, which is
mutate(BRBR = ifelse(is.na(BRBR_Survival) == T, 0,
ifelse(BRBR_Survival == 1 & BRBR == 0, 0, BRBR_Survival)),
LYER = ifelse(is.na(LYER_Survival) == T, 0,
ifelse(LYER_Survival == 1 & week < 17 & LYER == 0, 0, LYER_Survival)),
Mummy_Ptoids = ifelse(is.na(Mummy_Ptoids_Survival) == T, 0,
ifelse(Mummy_Ptoids_Survival == 1 & week > 3 & week < 17 & Mummy_Ptoids == 0, 0, Mummy_Ptoids_Survival))) %>%
select(cage, week, temp:temp_com, aop2_vs_AOP2:gsoh_vs_GSOH, BRBR, LYER, Mummy_Ptoids) %>%
pivot_longer(cols = c(BRBR, LYER, Mummy_Ptoids), names_to = "species", values_to = "survival") %>%
pivot_wider(names_from = "week", values_from = "survival")
# prepare for marked package
marked_df <- mark_recap_wide_df %>%
unite(`3`:`17`, col = "ch", sep = "", remove = F) %>% # encounter history
mutate(species = factor(species),
temp = factor(temp)) %>%
select(ch, species, temp, rich, aop2_vs_AOP2) %>%
as.data.frame()# process data (and set grouping variables)
marked.proc <- process.data(marked_df,
model = "CJS",
group = c("species","temp"))
# make design data (from processed data)
marked.ddl <- make.design.data(marked.proc)
# outline formulas for each parameter
Phi.form <- list(formula = ~species + temp + time + rich + aop2_vs_AOP2) # allow survival probability (Phi) to vary by species, temperature treatment, and time, as well as genetic diversity (rich) and AOP2 gene.
p.form <- list(formula=~species) # allow capture probability (p) to vary among species, which makes sense since some species are harder to detect than others (e.g. LYER prefers to hide underneath basal rosette compared to BRBR).
# make model
cjs.model <- crm(marked.proc,
marked.ddl,
model.parameters = list(Phi = Phi.form,
p = p.form),
accumulate = FALSE)
Number of evaluations: 100 -2lnl: 1584.559464
Number of evaluations: 200 -2lnl: 1545.61225
Number of evaluations: 300 -2lnl: 1545.1513
Number of evaluations: 400 -2lnl: 1544.552532
Number of evaluations: 500 -2lnl: 1544.263366
Number of evaluations: 600 -2lnl: 1543.941381
Number of evaluations: 700 -2lnl: 1543.877974
Number of evaluations: 800 -2lnl: 1543.737575
Number of evaluations: 900 -2lnl: 1543.702892
Number of evaluations: 1000 -2lnl: 1543.66882
Number of evaluations: 1100 -2lnl: 1543.614795
Number of evaluations: 1200 -2lnl: 1543.603781
Number of evaluations: 1300 -2lnl: 1543.599523
Number of evaluations: 1400 -2lnl: 1543.599059
Number of evaluations: 1500 -2lnl: 1543.598867
Number of evaluations: 1600 -2lnl: 1543.599425
Number of evaluations: 1700 -2lnl: 1543.59882
Number of evaluations: 1800 -2lnl: 1543.851487
Number of evaluations: 1900 -2lnl: 1544.902988
Number of evaluations: 2000 -2lnl: 1543.747498
Number of evaluations: 2100 -2lnl: 1546.793019
Number of evaluations: 2200 -2lnl: 1543.93277
Number of evaluations: 2300 -2lnl: 1543.726229
Number of evaluations: 2400 -2lnl: 1543.608845
Number of evaluations: 2500 -2lnl: 1544.570045
Number of evaluations: 2600 -2lnl: 1544.099186
Number of evaluations: 2700 -2lnl: 1546.276039
Number of evaluations: 2800 -2lnl: 1543.60815
Number of evaluations: 2900 -2lnl: 1546.353122
Number of evaluations: 3000 -2lnl: 1543.602188
Number of evaluations: 3100 -2lnl: 1543.644212
Number of evaluations: 3200 -2lnl: 1544.812579
Number of evaluations: 3300 -2lnl: 1548.032423
Number of evaluations: 3400 -2lnl: 1543.885738
Number of evaluations: 3500 -2lnl: 1546.171215
Number of evaluations: 3600 -2lnl: 1543.598834
Number of evaluations: 3700 -2lnl: 1543.598829cjs.model
crm Model Summary
Npar : 22
-2lnL: 1543.599
AIC : 1587.599
Beta
Estimate
Phi.(Intercept) 2.1075944
Phi.speciesLYER 4.6259506
Phi.speciesMummy_Ptoids 3.6983825
Phi.temp3 -0.1453095
Phi.time2 -0.4730241
Phi.time3 -1.6930867
Phi.time4 -3.8969044
Phi.time5 -5.2260443
Phi.time6 -4.9839936
Phi.time7 -4.3399889
Phi.time8 -3.2078034
Phi.time9 -3.2543029
Phi.time10 -3.2087875
Phi.time11 -3.4619101
Phi.time12 -4.8579032
Phi.time13 -5.2365724
Phi.time14 7.9476065
Phi.rich 0.2408103
Phi.aop2_vs_AOP2 0.2159383
p.(Intercept) 3.7487139
p.speciesLYER -2.3084112
p.speciesMummy_Ptoids -1.9112810# calculate confidence intervals of parameters
cjs.CIs <- cjs.hessian(cjs.model)
Number of evaluations: 100 -2lnl: 1543.745759
Number of evaluations: 200 -2lnl: 1544.275828
Number of evaluations: 300 -2lnl: 1551.334614
Number of evaluations: 400 -2lnl: 1543.990456
Number of evaluations: 500 -2lnl: 1544.674277
Number of evaluations: 600 -2lnl: 1543.849803
Number of evaluations: 700 -2lnl: 1543.986051
Number of evaluations: 800 -2lnl: 1543.619508
Number of evaluations: 900 -2lnl: 1543.967571
Number of evaluations: 1000 -2lnl: 1543.681394
Number of evaluations: 1100 -2lnl: 1545.519206
Number of evaluations: 1200 -2lnl: 1543.839998
Number of evaluations: 1300 -2lnl: 1543.872825
Number of evaluations: 1400 -2lnl: 1543.613415
Number of evaluations: 1500 -2lnl: 1543.790969
Number of evaluations: 1600 -2lnl: 1545.001866
Number of evaluations: 1700 -2lnl: 1552.77495
Number of evaluations: 1800 -2lnl: 1543.901213
Number of evaluations: 1900 -2lnl: 1545.847455
Number of evaluations: 2000 -2lnl: 1543.732566# we're primarily interested in the robustness of the aop2 vs AOP2 effect so we focus on it here
cjs.CIs$results$beta$Phi[c("aop2_vs_AOP2")] # logit scaleaop2_vs_AOP2
0.2159383 exp(cjs.CIs$results$beta$Phi[c("rich","aop2_vs_AOP2")]) # aop2 increases odds of survival by 24% relative to the average allele. rich aop2_vs_AOP2
1.272280 1.241026 exp(0.4161482) # upper CI (up to 52%)[1] 1.516111exp(0.0157780) # lower CI (as low as 1.5%)[1] 1.015903
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] parallel stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] survival_3.2-13 marked_1.2.6 lme4_1.1-27.1 Matrix_1.4-0
[5] forcats_0.5.1 stringr_1.4.0 dplyr_1.0.7 purrr_0.3.4
[9] readr_2.1.1 tidyr_1.1.4 tibble_3.1.6 ggplot2_3.3.5
[13] tidyverse_1.3.1 workflowr_1.6.2
loaded via a namespace (and not attached):
[1] nlme_3.1-152 fs_1.5.2 bit64_4.0.5
[4] lubridate_1.8.0 httr_1.4.2 rprojroot_2.0.2
[7] numDeriv_2016.8-1.1 tools_4.1.2 TMB_1.7.22
[10] backports_1.4.1 bslib_0.3.1 utf8_1.2.2
[13] R6_2.5.1 DBI_1.1.2 colorspace_2.0-2
[16] withr_2.4.3 tidyselect_1.1.1 bit_4.0.4
[19] compiler_4.1.2 git2r_0.28.0 cli_3.1.0
[22] rvest_1.0.2 expm_0.999-6 xml2_1.3.3
[25] sass_0.4.0 scales_1.1.1 digest_0.6.29
[28] minqa_1.2.4 rmarkdown_2.11 pkgconfig_2.0.3
[31] htmltools_0.5.2 dbplyr_2.1.1 fastmap_1.1.0
[34] rlang_0.4.12 readxl_1.3.1 rstudioapi_0.13
[37] jquerylib_0.1.4 R2admb_0.7.16.2 generics_0.1.1
[40] jsonlite_1.7.2 vroom_1.5.7 magrittr_2.0.1
[43] Rcpp_1.0.7 munsell_0.5.0 fansi_1.0.0
[46] lifecycle_1.0.1 stringi_1.7.3 whisker_0.4
[49] yaml_2.2.1 MASS_7.3-54 grid_4.1.2
[52] promises_1.2.0.1 crayon_1.4.2 lattice_0.20-45
[55] haven_2.4.3 splines_4.1.2 hms_1.1.1
[58] knitr_1.37 pillar_1.6.4 boot_1.3-28
[61] reprex_2.0.1 glue_1.6.0 evaluate_0.14
[64] data.table_1.14.2 modelr_0.1.8 vctrs_0.3.8
[67] nloptr_1.2.2.3 tzdb_0.2.0 httpuv_1.6.5
[70] cellranger_1.1.0 gtable_0.3.0 assertthat_0.2.1
[73] xfun_0.29 broom_0.7.11 coda_0.19-4
[76] later_1.3.0 truncnorm_1.0-8 optimx_2021-10.12
[79] ellipsis_0.3.2