Rensch’s Rule in Costa Rican hummingbirds
University of Costa Rica

Author

Marcelo Araya-Salas

Published

February 13, 2023

Source code, data and version control of the project found at https://github.com/maRce10/renschs_rule_costa_rican_hummingbirds

 

Purpose

  • Explore sexual size dimorphism in hummingbird species from Costa Rica

  • Evaluate Rensch’s Rule

 

Report overview

 

1 Explore data

Code 
# read data
dat <- as.data.frame(read_excel("./data/raw/supplementary data.xlsx"))

dat$Clade[dat$Clade == "Billiants"] <- "Brilliants"

# read trees
trees <- read.tree("./data/raw/Posterior trees fixed names 339 species 1 swift.trees")

dat$`scientific name`[dat$`scientific name` == "Phaetornis guy"] <- "Phaethornis guy"
dat$`scientific name`[dat$`scientific name` == "Phaetornis longirostris"] <- "Phaethornis longirostris"
dat$`scientific name`[dat$`scientific name` == "Heliomaster longisrostris"] <- "Heliomaster longirostris"
dat$scientific_name <- gsub(" ", "_", dat$`scientific name`)

trees <- lapply(trees, function(tree) {
    tree$tip.label[tree$tip.label == "Amazilia_amabilis"] <- "Polyerata_amabilis"
    tree$tip.label[tree$tip.label == "Hylocharis_eliciae"] <- "Chlorestes_eliciae"
    tree$tip.label[tree$tip.label == "Chlorostilbon_canivetii"] <- "Cynanthus_canivetii"
    tree$tip.label[tree$tip.label == "Elvira_cupreiceps"] <- "Microchera_cupreiceps"
    tree$tip.label[tree$tip.label == "Calliphlox_bryantae"] <- "Philodice_bryantae"
    tree$tip.label[tree$tip.label == "Amazilia_edward"] <- "Saucerottia_edward"
    tree$tip.label[tree$tip.label == "Amazilia_saucerottei"] <- "Saucerottia_hoffmanni"
    tree$tip.label[tree$tip.label == "Amazilia_saucerottei"] <- "Saucerottia_hoffmanni"
    tree$tip.label[tree$tip.label == "Amazilia_candida"] <- "Chlorestes_candida"
    tree$tip.label[tree$tip.label == "Elvira_chionura"] <- "Microchera_chonura"

    subtree <- drop.tip(phy = tree, tip = setdiff(tree$tip.label,
        unique(dat$scientific_name)))

    return(tree)
})


names(dat) <- gsub(" ", ".", names(dat))

# maximum clade credibility consensus tree
mcd_tree <- phangorn::maxCladeCred(trees)

mcd_tree$tip.label <- gsub("_", " ", mcd_tree$tip.label)

ggtree(mcd_tree, layout = "circular", col = viridis(10)[7]) + geom_tiplab(size = 2) +
    theme(plot.margin = unit(c(30, 10, 30, 10), "mm"))

Code 
# add log 10 variables
dat$lg10.female.weight <- log10(dat$female.weight)
dat$lg10.male.weight <- log10(dat$male.weight)

2 Run statistical models

Code 
prior <- c(prior(normal(0, 10), "b"), prior(normal(0, 50), "Intercept"),
    prior(student_t(3, 0, 20), "sd"), prior(student_t(3, 0, 20), "sigma"))

iter <- 10000

chains <- 1

cores <- 1

2.1 Size comparison among sexes

2.1.1 log(females) ~ log(males)

(note that in the data supplied ln.SEX.weight variables are in natural log)

Code 
male_x_mod <- phylogenetic_uncertainty(ln.female.weight ~ ln.male.weight,
    data = dat, sp.id.column = "scientific_name", phylos = trees,
    family = gaussian(), prior = prior, iter = 1000, control = list(adapt_delta = 0.99,
        max_treedepth = 15), model.name = "male_vs_female_size", save.models = FALSE,
    save.combined = TRUE, path = "./data/processed/", chains = chains,
    cores = cores)
Code 
html_summary(read.file = "./data/processed/male_vs_female_size.rds",
    highlight = TRUE)

2.2 male_vs_female_size

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20) ln.female.weight ~ ln.male.weight + (1 | gr(scientific_name, cov = vcv.phylo)) 1000 100 1 500 1 0 25587.67 25143.11 1984988351
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_Intercept 0.198 0.038 0.369 1.003 25587.67 25143.11
b_ln.male.weight 0.838 0.780 0.896 1.001 39303.85 36999.33

2.2.1 log10(females) ~ log10(males)

Code 
log10_male_x_mod <- phylogenetic_uncertainty(lg10.female.weight ~
    lg10.male.weight, data = dat, sp.id.column = "scientific_name",
    phylos = trees, family = gaussian(), prior = prior, iter = 1000,
    control = list(adapt_delta = 0.99, max_treedepth = 15), model.name = "log10_male_vs_female_size",
    save.models = FALSE, save.combined = TRUE, path = "./data/processed/",
    chains = chains, cores = cores)
Code 
html_summary(read.file = "./data/processed/log10_male_vs_female_size.rds",
    highlight = TRUE)

2.3 log10_male_vs_female_size

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20) lg10.female.weight ~ lg10.male.weight + (1 | gr(scientific_name, cov = vcv.phylo)) 1000 100 1 500 0 0 27000.11 27321.99 91716340
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_Intercept 0.086 0.016 0.160 1.003 27000.11 27321.99
b_lg10.male.weight 0.838 0.780 0.896 1.001 43152.52 37327.16

2.3.1 log(males) ~ log(females)

Code 
phylogenetic_uncertainty(ln.male.weight ~ ln.female.weight, data = dat,
    sp.id.column = "scientific_name", phylos = trees, family = gaussian(),
    prior = prior, iter = 1000, control = list(adapt_delta = 0.99,
        max_treedepth = 15), model.name = "female_vs_male_size", save.models = FALSE,
    save.combined = TRUE, path = "./data/processed/", chains = chains,
    cores = cores)
Code 
html_summary(read.file = "./data/processed/female_vs_male_size.rds",
    highlight = TRUE)

2.4 female_vs_male_size

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20) ln.male.weight ~ ln.female.weight + (1 | gr(scientific_name, cov = vcv.phylo)) 1000 100 1 500 0 0 28352.56 29483.41 1107860310
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_Intercept -0.148 -0.347 0.044 1.002 28352.56 29483.41
b_ln.female.weight 1.139 1.060 1.216 1.001 44292.65 37977.85

2.4.1 log10(males) ~ log10(females)

Code 
phylogenetic_uncertainty(lg10.male.weight ~ lg10.female.weight, data = dat,
    sp.id.column = "scientific_name", phylos = trees, family = gaussian(),
    prior = prior, iter = 1000, control = list(adapt_delta = 0.99,
        max_treedepth = 15), model.name = "log10_female_vs_male_size",
    save.models = FALSE, save.combined = TRUE, path = "./data/processed/",
    chains = chains, cores = cores)
Code 
html_summary(read.file = "./data/processed/log10_female_vs_male_size.rds",
    highlight = TRUE)

2.5 log10_female_vs_male_size

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20) lg10.male.weight ~ lg10.female.weight + (1 | gr(scientific_name, cov = vcv.phylo)) 1000 100 1 500 1 0 26779.79 29359.06 1451555959
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_Intercept -0.064 -0.151 0.020 1.002 26779.79 29359.06
b_lg10.female.weight 1.139 1.060 1.217 1.001 41417.59 38690.55

2.6 SSD (Lovich-Gibbons ratio) vs male size

  • SSD represented as aboslute Lovich-Gibbons ratio
  • sex size as log of weight

2.6.1 Absolute SSD vs log male size

Code 
phylogenetic_uncertainty(aboluteSSD ~ ln.male.weight, data = dat,
    sp.id.column = "scientific_name", phylos = trees, family = gaussian(),
    prior = prior, iter = 1000, control = list(adapt_delta = 0.99,
        max_treedepth = 15), model.name = "male_vs_abs_SSD", save.models = FALSE,
    save.combined = TRUE, path = "./data/processed/", chains = chains,
    cores = cores)
Code 
html_summary(read.file = "./data/processed/male_vs_abs_SSD.rds", highlight = TRUE)

2.7 male_vs_abs_SSD

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20) aboluteSSD ~ ln.male.weight + (1 | gr(scientific_name, cov = vcv.phylo)) 1000 100 1 500 11 0 29241.39 25006.14 664937249
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_Intercept 0.012 -0.128 0.136 1.002 29241.39 25006.14
b_ln.male.weight 0.065 0.007 0.123 1.002 52993.88 37500.93

2.7.1 Absolute SSD vs log10 male size

Code 
phylogenetic_uncertainty(aboluteSSD ~ lg10.male.weight, data = dat,
    sp.id.column = "scientific_name", phylos = trees, family = gaussian(),
    prior = prior, iter = 1000, control = list(adapt_delta = 0.99,
        max_treedepth = 15), model.name = "log10_male_vs_abs_SSD",
    save.models = FALSE, save.combined = TRUE, path = "./data/processed/",
    chains = chains, cores = cores)
Code 
html_summary(read.file = "./data/processed/log10_male_vs_abs_SSD.rds",
    highlight = TRUE)

2.8 log10_male_vs_abs_SSD

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20) aboluteSSD ~ lg10.male.weight + (1 | gr(scientific_name, cov = vcv.phylo)) 1000 100 1 500 5 0 29625.52 28350.92 1057511118
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_Intercept 0.011 -0.128 0.136 1.002 29625.52 28350.92
b_lg10.male.weight 0.150 0.017 0.285 1.001 50057.46 37967.99

2.8.1 absolute SSD vs log female size

Code 
phylogenetic_uncertainty(aboluteSSD ~ ln.female.weight, data = dat,
    sp.id.column = "scientific_name", phylos = trees, family = gaussian(),
    prior = prior, iter = 1000, control = list(adapt_delta = 0.99,
        max_treedepth = 15), model.name = "female_vs_abs_SSD", save.models = FALSE,
    save.combined = TRUE, path = "./data/processed/", chains = chains,
    cores = cores)
Code 
html_summary(read.file = "./data/processed/female_vs_abs_SSD.rds",
    highlight = TRUE)

2.9 female_vs_abs_SSD

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20) aboluteSSD ~ ln.female.weight + (1 | gr(scientific_name, cov = vcv.phylo)) 1000 100 1 500 4 0 35372.25 30277.62 1400504111
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_Intercept 0.053 -0.093 0.189 1.002 35372.25 30277.62
b_ln.female.weight 0.042 -0.027 0.112 1.002 60413.93 38792.66

2.9.1 absolute SSD vs log10 female size

Code 
phylogenetic_uncertainty(aboluteSSD ~ lg10.female.weight, data = dat,
    sp.id.column = "scientific_name", phylos = trees, family = gaussian(),
    prior = prior, iter = 1000, control = list(adapt_delta = 0.99,
        max_treedepth = 15), model.name = "log10_female_vs_abs_SSD",
    save.models = FALSE, save.combined = TRUE, path = "./data/processed/",
    chains = chains, cores = cores)
Code 
html_summary(read.file = "./data/processed/log10_female_vs_abs_SSD.rds",
    highlight = TRUE)

2.10 log10_female_vs_abs_SSD

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20) aboluteSSD ~ lg10.female.weight + (1 | gr(scientific_name, cov = vcv.phylo)) 1000 100 1 500 2 0 33658.82 28921.81 1528172800
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_Intercept 0.053 -0.095 0.189 1.001 33658.82 28921.81
b_lg10.female.weight 0.096 -0.064 0.259 1.001 57926.09 37884.08

3 Graphs

3.1 Male vs female body mass

Code 
fit <- readRDS("./data/processed/log10_male_vs_female_size.rds")

gg_fit <- conditional_effects(fit)

pred_dat <- gg_fit$lg10.male.weight
pred_dat$Clade <- "Hermits"

fit_summ <- draw_summary(posterior::as_draws_array(fit), variables = c("b_Intercept",
    "b_lg10.male.weight"), robust = TRUE, probs = c(0.025, 0.975))

ggplot(data = dat, mapping = aes(x = lg10.male.weight, y = lg10.female.weight,
    fill = Clade, shape = Clade)) + scale_fill_viridis_d(alpha = 0.5) +
    geom_segment(data = pred_dat, mapping = aes(x = min(pred_dat$lg10.male.weight),
        y = min(pred_dat$estimate__), xend = max(pred_dat$lg10.male.weight),
        yend = max(pred_dat$estimate__))) + geom_ribbon(data = pred_dat,
    aes(ymin = lower__, ymax = upper__), fill = "gray", alpha = 0.3) +
    scale_shape_manual(values = c(21, 22, 23, 21, 22, 23, 21, 22)) +
    geom_point(size = 3, color = "transparent") + theme_classic(base_size = 20) +
    labs(x = "log10(male body mass)", y = "log10(female body mass)")

Code 
# ggsave(filename =
# './output/log10_male_vs_absolute_SSD_70dpi.tiff', dpi = 70)
# ggsave(filename =
# './output/log10_male_vs_absolute_SSD_300dpi.tiff', dpi = 300)

3.2 Male body mass vs absolute SSD

Code 
fit <- readRDS("./data/processed/log10_male_vs_abs_SSD.rds")

gg_fit <- conditional_effects(fit)

pred_dat <- gg_fit$lg10.male.weight
pred_dat$Clade <- "Hermits"

fit_summ <- draw_summary(posterior::as_draws_array(fit), variables = c("b_Intercept",
    "b_lg10.male.weight"), robust = TRUE, probs = c(0.025, 0.975))

ggplot(data = dat, mapping = aes(x = lg10.male.weight, y = aboluteSSD,
    fill = Clade, shape = Clade)) + scale_fill_viridis_d(alpha = 0.5) +
    geom_segment(data = pred_dat, mapping = aes(x = min(pred_dat$lg10.male.weight),
        y = min(pred_dat$estimate__), xend = max(pred_dat$lg10.male.weight),
        yend = max(pred_dat$estimate__))) + geom_ribbon(data = pred_dat,
    aes(ymin = lower__, ymax = upper__), fill = "gray", alpha = 0.3) +
    scale_shape_manual(values = c(21, 22, 23, 21, 22, 23, 21, 22)) +
    geom_point(color = viridis(10, alpha = 0.6)[7], size = 3) + theme_classic(base_size = 20) +
    labs(x = "log10(male body mass)", y = "Dimorphism (absolute SSD)")

Code 
# ggsave(filename =
# './output/log10_female_vs_absolute_SSD_300dpi.tiff', dpi =
# 300) ggsave(filename =
# './output/log10_female_vs_absolute_SSD_70dpi.tiff', dpi = 70)

3.3 Female body mass vs absolute SSD

Code 
fit <- readRDS("./data/processed/log10_female_vs_abs_SSD.rds")

gg_fit <- conditional_effects(fit)

pred_dat <- gg_fit$lg10.female.weight
pred_dat$Clade <- "Hermits"

fit_summ <- draw_summary(posterior::as_draws_array(fit), variables = c("b_Intercept",
    "b_lg10.female.weight"), robust = TRUE, probs = c(0.025, 0.975))

ggplot(data = dat, mapping = aes(x = lg10.female.weight, y = aboluteSSD,
    fill = Clade, shape = Clade)) + scale_fill_viridis_d(alpha = 0.5) +
    # geom_segment(data = pred_dat, mapping = aes(x =
    # min(pred_dat$lg10.female.weight), y =
    # min(pred_dat$estimate__), xend =
    # max(pred_dat$lg10.female.weight), yend =
    # max(pred_dat$estimate__))) + geom_ribbon(data = pred_dat,
    # aes(ymin = lower__, ymax = upper__), fill = 'gray', alpha
    # = 0.3) +
scale_shape_manual(values = c(21, 22, 23, 21, 22, 23, 21, 22)) + geom_point(color = viridis(10,
    alpha = 0.6)[7], size = 3) + theme_classic(base_size = 20) + labs(x = "log10(female body mass)",
    y = "Dimorphism (absolute SSD)")

Code 
# ggsave(filename =
# './output/male_vs_female_size_10log_300dpi.tiff', dpi = 300)
# ggsave(filename =
# './output/male_vs_female_size_10log_70dpi.tiff', dpi = 70)
Code 
fill_color <- viridis(10)[7]

dat$Clade <- factor(dat$Clade, levels = c("Bees", "Coquettes", "Hermits", "Brilliants", "Emeralds", "Mangoes", "Mountain Gems", "Topazes"))

agg_dat <- aggregate(`scientific_name` ~ Clade, data = dat, length)
agg_dat$`Lovich-Gibbons.ratio` <- NA
agg_dat$n.labels <- paste("n =", agg_dat$`scientific_name`)

# set.seed(12)
trans <-  function(height = 0, ...) ggplot2::position_jitter(height = height, ...)

# composed box plot
ggplot(dat, aes(y = `Lovich-Gibbons.ratio`, x = Clade)) +
    geom_hline(yintercept = 0, lty = 2) +     
## add half-violin from {ggdist} package
  ggdist::stat_halfeye(
    fill = fill_color,
    alpha = 0.5,
    ## custom bandwidth
    adjust = .5,
    ## adjust height
    width = .6,
    .width = 0,
    ## move geom to the cright
    justification = -.2,
    point_colour = NA
  ) +
  geom_boxplot(fill = fill_color,
    width = .15,
    ## remove outliers
    outlier.shape = NA ## `outlier.shape = NA` works as well
  ) +
  ## add justified jitter from the {gghalves} package
  gghalves::geom_half_point(
    color = fill_color,
    ## draw jitter on the left
    side = "l",
    ## control range of jitter
    range_scale = .4,
    ## add some transparency
    alpha = .5,
    transformation = ggplot2::position_jitter(height = 0)

  ) +
    labs(x="Clade", y= "Lovich-Gibbons ratio"
  ) +
  ylim(c(-0.39, 0.145)) +
  geom_text(data = agg_dat, aes(y = rep(-0.387, nrow(agg_dat)), x = Clade, label = n.labels), nudge_x = 0, size = 6) +
     theme_classic(base_size = 18) +
theme(axis.text.x = element_text(angle = 30, hjust = 1))

4 Combined model diagnostics

Code 
check_rds_models(path = "./data/processed", html = TRUE)
priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20) aboluteSSD ~ ln.female.weight + (1 | gr(scientific_name, cov = vcv.phylo)) 1000 100 1 500 4 0 12530.009 19870.545 1400504111
2 b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20) ln.male.weight ~ ln.female.weight + (1 | gr(scientific_name, cov = vcv.phylo)) 1000 100 1 500 0 0 7078.550 10497.000 1107860310
3 b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20) aboluteSSD ~ lg10.female.weight + (1 | gr(scientific_name, cov = vcv.phylo)) 1000 100 1 500 2 0 11334.133 17552.584 1528172800
4 b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20) lg10.male.weight ~ lg10.female.weight + (1 | gr(scientific_name, cov = vcv.phylo)) 1000 100 1 500 1 0 6814.181 9015.423 1451555959
5 b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20) aboluteSSD ~ lg10.male.weight + (1 | gr(scientific_name, cov = vcv.phylo)) 1000 100 1 500 5 0 10660.031 17028.943 1057511118
6 b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20) lg10.female.weight ~ lg10.male.weight + (1 | gr(scientific_name, cov = vcv.phylo)) 1000 100 1 500 0 0 7012.412 8034.363 91716340
7 b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20) aboluteSSD ~ ln.male.weight + (1 | gr(scientific_name, cov = vcv.phylo)) 1000 100 1 500 11 0 11299.206 18114.612 664937249
8 b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20) ln.female.weight ~ ln.male.weight + (1 | gr(scientific_name, cov = vcv.phylo)) 1000 100 1 500 1 0 7649.259 10986.574 1984988351

Takeaways

  • Dimorpism increases with body size in males but not in females

 

 

Session information

R version 4.1.0 (2021-05-18)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.2 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/atlas/libblas.so.3.10.3
LAPACK: /usr/lib/x86_64-linux-gnu/atlas/liblapack.so.3.10.3

locale:
 [1] LC_CTYPE=pt_BR.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=es_CR.UTF-8        LC_COLLATE=pt_BR.UTF-8    
 [5] LC_MONETARY=es_CR.UTF-8    LC_MESSAGES=pt_BR.UTF-8   
 [7] LC_PAPER=es_CR.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] gghalves_0.1.3      phangorn_2.7.1      ggplot2_3.4.0      
 [4] viridis_0.6.2       viridisLite_0.4.1   brms_2.18.0        
 [7] Rcpp_1.0.9          sketchy_1.0.2       brmsish_1.0.0      
[10] ggtree_3.2.1        readxl_1.3.1        ape_5.6-2          
[13] xaringanExtra_0.7.0 rprojroot_2.0.3     formatR_1.11       
[16] knitr_1.42          kableExtra_1.3.4   

loaded via a namespace (and not attached):
  [1] uuid_1.1-0           backports_1.4.1      fastmatch_1.1-0     
  [4] systemfonts_1.0.4    plyr_1.8.7           igraph_1.3.5        
  [7] lazyeval_0.2.2       splines_4.1.0        crosstalk_1.2.0     
 [10] TH.data_1.1-0        rstantools_2.2.0     inline_0.3.19       
 [13] digest_0.6.31        yulab.utils_0.0.5    htmltools_0.5.4     
 [16] fansi_1.0.3          magrittr_2.0.3       checkmate_2.1.0     
 [19] remotes_2.4.2        RcppParallel_5.1.5   matrixStats_0.62.0  
 [22] xts_0.12.2           sandwich_3.0-1       svglite_2.1.0       
 [25] prettyunits_1.1.1    colorspace_2.0-3     rvest_1.0.3         
 [28] ggdist_3.2.0         xfun_0.36            dplyr_1.0.10        
 [31] callr_3.7.3          crayon_1.5.2         jsonlite_1.8.4      
 [34] lme4_1.1-27.1        survival_3.2-11      zoo_1.8-11          
 [37] glue_1.6.2           gtable_0.3.1         emmeans_1.8.1-1     
 [40] webshot_0.5.4        distributional_0.3.1 pkgbuild_1.4.0      
 [43] rstan_2.21.7         BiocGenerics_0.40.0  abind_1.4-5         
 [46] scales_1.2.1         mvtnorm_1.1-3        DBI_1.1.1           
 [49] miniUI_0.1.1.1       xtable_1.8-4         gridGraphics_0.5-1  
 [52] tidytree_0.4.0       stats4_4.1.0         StanHeaders_2.21.0-7
 [55] DT_0.26              htmlwidgets_1.5.4    httr_1.4.4          
 [58] threejs_0.3.3        posterior_1.3.1      ellipsis_0.3.2      
 [61] pkgconfig_2.0.3      loo_2.4.1.9000       farver_2.1.1        
 [64] utf8_1.2.2           labeling_0.4.2       ggplotify_0.1.0     
 [67] tidyselect_1.2.0     rlang_1.0.6          reshape2_1.4.4      
 [70] later_1.3.0          munsell_0.5.0        cellranger_1.1.0    
 [73] tools_4.1.0          cli_3.6.0            generics_0.1.3      
 [76] ggridges_0.5.4       evaluate_0.20        stringr_1.5.0       
 [79] fastmap_1.1.0        yaml_2.3.7           processx_3.8.0      
 [82] purrr_1.0.0          packrat_0.9.0        pbapply_1.6-0       
 [85] nlme_3.1-152         mime_0.12            projpred_2.0.2      
 [88] aplot_0.1.7          xml2_1.3.3           compiler_4.1.0      
 [91] bayesplot_1.9.0      shinythemes_1.2.0    rstudioapi_0.14     
 [94] gamm4_0.2-6          treeio_1.18.1        tibble_3.1.8        
 [97] stringi_1.7.12       ps_1.7.2             Brobdingnag_1.2-9   
[100] lattice_0.20-44      Matrix_1.5-1         nloptr_1.2.2.2      
[103] markdown_1.3         shinyjs_2.1.0        tensorA_0.36.2      
[106] vctrs_0.5.2          pillar_1.8.1         lifecycle_1.0.3     
[109] bridgesampling_1.1-2 estimability_1.4.1   cowplot_1.1.1       
[112] httpuv_1.6.6         patchwork_1.1.2      R6_2.5.1            
[115] promises_1.2.0.1     gridExtra_2.3        codetools_0.2-18    
[118] boot_1.3-28          colourpicker_1.2.0   MASS_7.3-54         
[121] gtools_3.9.3         assertthat_0.2.1     withr_2.5.0         
[124] shinystan_2.6.0      multcomp_1.4-17      mgcv_1.8-36         
[127] parallel_4.1.0       quadprog_1.5-8       grid_4.1.0          
[130] ggfun_0.0.7          minqa_1.2.4          tidyr_1.1.3         
[133] coda_0.19-4          rmarkdown_2.20       Biobase_2.54.0      
[136] shiny_1.7.3          base64enc_0.1-3      dygraphs_1.1.1.6