Do females and males differ in organ scaling?

This supplementary file contains the R code for processing and analysing the raw data, and creating figures for the manuscript - Do females and males differ in organs scaling? A comparative analysis across vertebrates

Data subset used in this study is part of a database on organ sizes, which is available here: https://felixpleiva.github.io/organ_size_DB/. Additionally, I complemented the dataset on organ size by revisiting a recent publication by Tsuboi et al., 2018, which includes, among other information, sex- and species-specific data on organ size for numerous species of mammals, birds, and fish. Together, these approaches, allowed us to compile information for 204 vertebrate species and data on a total of 9 different organs.

# Clean working environment
rm(list = ls())

Documentation

Licence

This repository is provided by the authors under the Attribution-NonCommercial-NoDerivatives 4.0 International licence

Citation

When using the data and/or code associated with this repository, they should be cited as follows (almost never works😞….but just in case):

Leiva, F. P., A. J. Hendriks. (2026). Do females and males differ in organ scaling? A comparative analysis across vertebrates. Zenodo. DOI XXXX.

Contact

This script is authored by Félix P. Leiva. For any questions related to this resource, please contact me at the email address: felixpleiva@gmail.com.

Data and phylogenetic tree

Load and tidy data

dat <- read_xlsx("../outputs/dat_for_analyses_full_merged.xlsx") %>%
  rename( 
    sex = sex_individual,
    life_stage = life_stage_individual,
    original_category = trait_size_category,
    body_size = body_size_mean,
    organ_size = organ_size_mean,
    species = species_reported
  ) %>%
  mutate(organ_grouped = recode(organ_grouped,
                                pituitary_gland = "pituitary_glands",
                                lung = "lungs",
                                kidney = "kidneys")) %>%
  mutate(class = recode(class,
                        Actinopterygii = "Teleostei"))


# Calculate mean per species and sex. Most species were represented primarily by adults. Life stages that were not explicitly reported were assumed to correspond to adult individuals.

dat_mean <- dat %>%
  group_by(
    phylum, class, order, family, genus,
    species, sex, organ_grouped
  ) %>%
  summarise(
    body_size = mean(body_size, na.rm = TRUE),
    organ_size = mean(organ_size, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  mutate(
    log10_body_size    = log10(body_size),
    log10_organ_size = log10(organ_size),
    phylo = gsub(" ", "_", species)
  )

Next, from the preceding code, I will extract the list of species for which I want to retrieve the phylogenetic tree and export it to my working directory.

list_of_species <- dat_mean %>%
  distinct(species) %>%
  pull(species) %>%
  as.character()

writeLines(list_of_species, "../outputs/list_of_species.txt")

Now I have exported the list of species, I will visit the TimeTree of Life website and, in the lower section labelled Load a list of species, upload the file list_of_species.csv. Subsequently, I will export it in Newick format. Note that not all species from our list are included in the resulting tree, although most are represented.

Check data structure and transform variables if necessary

glimpse(dat_mean)

Rows: 676
Columns: 13
$ phylum           <chr> "Chordata", "Chordata", "Chordata", "Chordata", "Chor…
$ class            <chr> "Aves", "Aves", "Aves", "Aves", "Aves", "Aves", "Aves…
$ order            <chr> "Anseriformes", "Anseriformes", "Anseriformes", "Anse…
$ family           <chr> "Anatidae", "Anatidae", "Anatidae", "Anatidae", "Anat…
$ genus            <chr> "Anas", "Anas", "Anas", "Anas", "Anas", "Anas", "Anas…
$ species          <chr> "Anas acuta", "Anas acuta", "Anas crecca", "Anas crec…
$ sex              <chr> "female", "male", "female", "male", "female", "female…
$ organ_grouped    <chr> "brain", "brain", "brain", "brain", "adrenal_glands",…
$ body_size        <dbl> 721, 1176, 279, 347, 2921, 1200, 2921, 2921, 2921, 29…
$ organ_size       <dbl> 4.91, 5.28, 3.18, 2.90, 0.13, 5.44, 19.50, 17.40, 3.5…
$ log10_body_size  <dbl> 2.9, 3.1, 2.4, 2.5, 3.5, 3.1, 3.5, 3.5, 3.5, 3.5, 3.5…
$ log10_organ_size <dbl> 0.691, 0.723, 0.503, 0.463, -0.886, 0.736, 1.290, 1.2…
$ phylo            <chr> "Anas_acuta", "Anas_acuta", "Anas_crecca", "Anas_crec…

Variables to convert to factor

columns_to_factor <- c(
  "phylum", 
  "class", 
  "order", 
  "family", 
  "genus", 
  "species",
  "phylo",   
  "sex",
  "organ_grouped"
)

Variables to convert to numeric

columns_to_numeric <- c(
  "body_size", 
  "organ_size",
  "log10_organ_size",
  "log10_body_size"
)

Apply transformation

# Apply the conversions to the dataset.
dat_mean <- dat_mean %>%
  mutate(
    across(all_of(columns_to_factor), as.factor),
    across(all_of(columns_to_numeric), as.numeric)
  )

# Confirm that changes have been applied correctly
glimpse(dat_mean)

Rows: 676
Columns: 13
$ phylum           <fct> Chordata, Chordata, Chordata, Chordata, Chordata, Cho…
$ class            <fct> Aves, Aves, Aves, Aves, Aves, Aves, Aves, Aves, Aves,…
$ order            <fct> Anseriformes, Anseriformes, Anseriformes, Anseriforme…
$ family           <fct> Anatidae, Anatidae, Anatidae, Anatidae, Anatidae, Ana…
$ genus            <fct> Anas, Anas, Anas, Anas, Anas, Anas, Anas, Anas, Anas,…
$ species          <fct> Anas acuta, Anas acuta, Anas crecca, Anas crecca, Ana…
$ sex              <fct> female, male, female, male, female, female, female, f…
$ organ_grouped    <fct> brain, brain, brain, brain, adrenal_glands, brain, he…
$ body_size        <dbl> 721, 1176, 279, 347, 2921, 1200, 2921, 2921, 2921, 29…
$ organ_size       <dbl> 4.91, 5.28, 3.18, 2.90, 0.13, 5.44, 19.50, 17.40, 3.5…
$ log10_body_size  <dbl> 2.9, 3.1, 2.4, 2.5, 3.5, 3.1, 3.5, 3.5, 3.5, 3.5, 3.5…
$ log10_organ_size <dbl> 0.691, 0.723, 0.503, 0.463, -0.886, 0.736, 1.290, 1.2…
$ phylo            <fct> Anas_acuta, Anas_acuta, Anas_crecca, Anas_crecca, Ana…

Filter the dataframe to include only organs with at least five species per sex and organ. Some of them have only three species per sex. From this output, the thymus, large intestine, and small intestine are excluded.

# First identify qualifying organ-sex combinations
organs_five_species <- dat_mean %>%
  count(organ_grouped, species, sex) %>%
  group_by(organ_grouped, sex) %>% # by sex and organ
  summarise(n_species = n_distinct(species), .groups = "drop") %>%
  filter(n_species >= 5) # at least 5 spp 

# Then filter original data to these combinations only
dat_mean <- dat_mean %>%
  inner_join(organs_five_species %>% select(organ_grouped, sex), 
             by = c("organ_grouped", "sex"))

# Lets check whcih organs were retained after this filter
dat_mean %>%
  summarise(n_species = n_distinct(species), .by = organ_grouped)

# A tibble: 10 × 2
   organ_grouped    n_species
   <fct>                <int>
 1 brain                  246
 2 adrenal_glands           7
 3 heart                   12
 4 kidneys                  8
 5 liver                   14
 6 pancreas                 6
 7 spleen                  11
 8 stomach                  8
 9 pituitary_glands         7
10 lungs                    8

In addition to the five-species filter, we applied a ten-species filter when the organ in question was the brain. This was done to make our analyses more comparable to those recently published by Tsuboi et al. (2018). Moreover, this approach allowed us to assess the robustness and granularity of our results within a subset of species for which more data were available. The ten-species filter encompassed a total of two families of fish, five families of birds, and three families of mammals.

# Identify qualifying families (>10 unique species per family) for brain only
good_families_brain <- dat_mean %>%
  filter(organ_grouped == "brain") %>%
  group_by(family) %>%
  summarise(n_species = n_distinct(species), .groups = "drop") %>%
  filter(n_species > 10) %>%
  pull(family)

# Filter main dataset: keep ALL non-brain data + brain data for qualifying families only
dat_mean<- dat_mean %>%
  filter(organ_grouped != "brain" | family %in% good_families_brain) %>%
  droplevels() # here i drop the rows filtered out

# Species counts per organ, class, order, family
species_summary <- dat_mean %>%
  distinct(organ_grouped, class, order, family, species, .keep_all = TRUE) %>%
  count(organ_grouped, class, order, family, name = "n_species_unique") %>%
  arrange(organ_grouped, class, desc(n_species_unique))

One of the things I discovered recently is the library (DT), which enables the creation of these interactive tables that allow for much more dynamic filtering. Super cool!

datatable(
  species_summary,
  filter = "top", 
  options = list(
    pageLength = 15,
    autoWidth = TRUE,
    dom = 'Bfrtip', 
    buttons = c('copy', 'csv', 'excel')
  ),
  colnames = c("Organ", "Class", "Order", "Family", "Species count")
)

Load and prune phylogenetic tree

tree <- read.tree("../outputs/phylogenetic_tree.nwk")

# Number of tips (species)
length(tree$tip.label)

[1] 277

Our next step is to prune the tree to retain only those species for which organ size measurements exist for both males and females. This pruning is important because our phylogenetic informed analyses require the calculation of a variance-covariance matrix based on the shared evolutionary history of the species.

# Check for database/tree species mismatches
setdiff(tree$tip.label, dat_mean$phylo)      # Tree species absent in database

 [1] "Pseudosimochromis_babaulti"            
 [2] "Oostethus_brachyurus"                  
 [3] "Dunckerocampus_dactyliophorus"         
 [4] "Madoqua_kirkii"                        
 [5] "Philantomba_monticola"                 
 [6] "Aotus_trivirgatus"                     
 [7] "Cebuella_pygmaea"                      
 [8] "Callithrix_penicillata"                
 [9] "Callimico_goeldii"                     
[10] "Leontopithecus_chrysomelas"            
[11] "Leontopithecus_rosalia"                
[12] "Leontocebus_fuscicollis_illigeri"      
[13] "Saguinus_mystax"                       
[14] "Saguinus_labiatus_labiatus"            
[15] "Saguinus_leucopus"                     
[16] "Saguinus_midas"                        
[17] "Saguinus_niger"                        
[18] "Saguinus_oedipus"                      
[19] "Saguinus_geoffroyi"                    
[20] "Cebus_olivaceus"                       
[21] "Cebus_capucinus"                       
[22] "Sapajus_xanthosternos"                 
[23] "Cebus_nigritus_robustus"               
[24] "Saimiri_sciureus_macrodon"             
[25] "Saimiri_sciureus_sciureus"             
[26] "Saimiri_sciureus_cassiquiarensis"      
[27] "Chlorocebus_pygerythrus_pygerythrus"   
[28] "Cercopithecus_mitis_stuhlmanni"        
[29] "Cercopithecus_albogularis_kolbi"       
[30] "Cercopithecus_albogularis_erythrarchus"
[31] "Cercopithecus_mona"                    
[32] "Cercopithecus_neglectus"               
[33] "Papio_hamadryas"                       
[34] "Theropithecus_gelada"                  
[35] "Macaca_radiata_radiata"                
[36] "Macaca_nemestrina"                     
[37] "Macaca_leonina"                        
[38] "Semnopithecus_entellus"                
[39] "Semnopithecus_priam"                   
[40] "Trachypithecus_vetulus_vetulus"        
[41] "Trachypithecus_phayrei_shanicus"       
[42] "Trachypithecus_phayrei_crepuscula"     
[43] "Presbytis_thomasi"                     
[44] "Presbytis_rubicunda_rubicunda"         
[45] "Presbytis_melalophos_sumatranus"       
[46] "Presbytis_femoralis_percura"           
[47] "Nasalis_larvatus"                      
[48] "Procolobus_verus"                      
[49] "Piliocolobus_tephrosceles"             
[50] "Piliocolobus_kirkii"                   
[51] "Colobus_vellerosus"                    
[52] "Colobus_polykomos"                     
[53] "Colobus_angolensis_palliatus"          
[54] "Psephotellus_varius"                   
[55] "Orthopsittaca_manilatus"               
[56] "Lophochroa_leadbeateri"                
[57] "Zanda_funerea"                         
[58] "Crithagra_dorsostriata"                
[59] "Crithagra_burtoni"                     
[60] "Bathilda_ruficauda"                    
[61] "Granatina_ianthinogaster"              
[62] "Chloebia_gouldiae"                     
[63] "Spinus_spinus"                         
[64] "Linaria_cannabina"                     
[65] "Crithagra_mozambica"                   
[66] "Spinus_cucullatus"                     
[67] "Acanthis_flammea"                      
[68] "Chloris_spinoides"                     
[69] "Chloris_ambigua"                       
[70] "Chloris_chloris"                       
[71] "Aythya_fuligula"                       
[72] "Tetrastes_bonasia"                     
[73] "Gallus_lafayettii"                     

setdiff(dat_mean$phylo, tree$tip.label)      # Database species absent in tree

 [1] "Nyroca_fuligula"              "Ammoperdix_griseogularis"    
 [3] "Bonasa_bonasia"               "Gallus_lafayetii"            
 [5] "Lagopus_mutus"                "Erythrura_gouldiae"          
 [7] "Lonchura_grandis"             "Lonchura_molucca"            
 [9] "Lonchura_pallida"             "Lonchura_quinticolor"        
[11] "Lonchura_spectabilis"         "Lonchura_tristissima"        
[13] "Mandingoa_nitidula"           "Neochmia_ruficauda"          
[15] "Pyrenestes_ostrinus"          "Uraeginthus_ianthinogaster"  
[17] "Carduelis_ambigua"            "Carduelis_cannabina"         
[19] "Carduelis_chloris"            "Carduelis_cucullata"         
[21] "Carduelis_flammea"            "Carduelis_spinoides"         
[23] "Carduelis_spinus"             "Serinus_burtoni"             
[25] "Serinus_dorsostriatus"        "Serinus_mozambicus"          
[27] "Cacatua_leadbeateri"          "Calyptorhynchus_funereus"    
[29] "Orthopsittaca_manilata"       "Polytelis_swainsonii"        
[31] "Psephotus_varius"             "Cephalophus_monticola"       
[33] "Eudorcas_thomsoni"            "Madoqua_kirki"               
[35] "Saimiri_sciureus"             "Chlorocebus_aethiops_sabaeus"
[37] "Simochromis_babaulti"         "Doryichthys_martensi"        
[39] "Doryrhamphus_dactyliophorus"  "Microphis_brachyurus"        
[41] "Trachyhampus_serratus"       

# Exclude mismatched species for consistency
tree <- keep.tip(tree, intersect(tree$tip.label, dat_mean$phylo))
dat_mean <- dat_mean[dat_mean$phylo %in% tree$tip.label, ]

# Re-check for mismatches after pruning
setdiff(tree$tip.label, dat_mean$phylo)

character(0)

setdiff(dat_mean$phylo, tree$tip.label)

character(0)

# Is ultrametric?
is.ultrametric(tree)

[1] TRUE

# how many species?
length(tree$tip.label)

[1] 204

# Compute phylogenetic covariance matrix (Brownian motion)
A <- vcv.phylo(tree, corr = TRUE)

How many species have phylogenetic information available across all taxonomic classes for the brain?

dat_mean %>%
  count(organ_grouped, species, class) %>%
  group_by(organ_grouped) %>%
  filter(organ_grouped == "brain") %>%
  summarise(
    n_species = n_distinct(species),
    .groups = 'drop'
  ) %>%
  arrange(desc(n_species))

# A tibble: 1 × 2
  organ_grouped n_species
  <fct>             <int>
1 brain               189

dat_mean %>%
  count(organ_grouped, species) %>%
  group_by(organ_grouped) %>%
  count(organ_grouped, name = "n_species") %>%
  arrange(desc(n_species))

# A tibble: 10 × 2
# Groups:   organ_grouped [10]
   organ_grouped    n_species
   <fct>                <int>
 1 brain                  189
 2 liver                   12
 3 heart                   10
 4 spleen                   9
 5 stomach                  8
 6 kidneys                  6
 7 lungs                    6
 8 adrenal_glands           5
 9 pituitary_glands         5
10 pancreas                 4

Since our objective is to account for the evolutionary history of the species, it is essential to recognise that, in a comparative context, species cannot be treated as independent units of observation. Closely related species tend to exhibit similar values for a given trait due to shared ancestry. Although no strict rule exists for the minimum number of species required, Garland and Adolph’s paper Why Not to Do Two-Species Comparative Studies: Limitations on Inferring Adaptation’ discusses relevant limitations and recommends broader sampling. As a rule of thumb, we evaluated the effects of accounting for shared evolutionary history in organs measured across more than 5 species, and compared these results with non-phylogenetic models (onward refereed as simple models). It means we won´t include pancreas size in the analysis.

Exclude pancreas from the analysis (only four species)

dat_mean <- dat_mean %>% 
  filter(organ_grouped != "pancreas")

# Confirm pancreas is gone
unique(dat_mean$organ_grouped)

[1] brain            adrenal_glands   heart            kidneys         
[5] liver            spleen           stomach          pituitary_glands
[9] lungs           
10 Levels: adrenal_glands brain heart kidneys liver lungs ... stomach

and now re-run the previous lines para merge the tree and data

# Check for database/tree species mismatches
setdiff(tree$tip.label, dat_mean$phylo)      # Tree species absent in database

character(0)

setdiff(dat_mean$phylo, tree$tip.label)      # Database species absent in tree

character(0)

# Exclude mismatched species for consistency
tree <- keep.tip(tree, intersect(tree$tip.label, dat_mean$phylo))
dat_mean <- dat_mean[dat_mean$phylo %in% tree$tip.label, ]

# Re-check for mismatches after pruning
setdiff(tree$tip.label, dat_mean$phylo)

character(0)

setdiff(dat_mean$phylo, tree$tip.label)

character(0)

# Is ultrametric?q
is.ultrametric(tree)

[1] TRUE

# how many species?
length(tree$tip.label)

[1] 204

# Compute phylogenetic covariance matrix (Brownian motion)
A <- vcv.phylo(tree, corr = TRUE)

Lets count species again by organ

dat_mean %>%
  count(organ_grouped, species) %>%
  group_by(organ_grouped) %>%
  count(organ_grouped, name = "n_species") %>%
  arrange(desc(n_species))

# A tibble: 9 × 2
# Groups:   organ_grouped [9]
  organ_grouped    n_species
  <fct>                <int>
1 brain                  189
2 liver                   12
3 heart                   10
4 spleen                   9
5 stomach                  8
6 kidneys                  6
7 lungs                    6
8 adrenal_glands           5
9 pituitary_glands         5

and taxonomic Class

dat_mean %>%
  count(class, species) %>%
  group_by(class) %>%
  count(class, name = "n_species") %>%
  arrange(desc(n_species))

# A tibble: 4 × 2
# Groups:   class [4]
  class     n_species
  <fct>         <int>
1 Aves             98
2 Teleostei        70
3 Mammalia         35
4 Reptilia          1

Now the tree match the species in the dataframe, we will create a figure of the phylogenetic tree containing the 204 species for which we have organ size data for both males and females. We will add a plot indicating whether each organ was measured (green) or not (grey).

Important note!

Our analyses for most organs are based on a relatively small number of observations (i.e., species). This may limit the interpretability of our models. Cross-validation (via the loo_compare function) was implemented in some cases following the recommendations in the cross-validation FAQ, particularly regarding influential observations. In the preceding link Aki Vehtari states: Pareto-k is also useful as a measure of influence of an observation. Highly influential observations have high k^ values. Very high k values often indicate model misspecification, outliers or mistakes in data processing. For this reason, I decided to include additional species by searching for complementary datasets. We incorporated the dataset from Tsuboi et al., which, to our knowledge, represents the most extensive available compilation of vertebrate brain size. I filtered the entries where sex was reported for same species.

Figure 1

The figure below illustrates the expected response scenarios for the scaling relationships between organ size and body mass in males and females.

# Set colors by sex 
cols <- c(male = "#008B8B", female = "#9370DB")

# Implemented  on a previous script: https://felixpleiva.github.io/correlated_evolution_moths/

sex_symbols <- tibble(
  sex    = c("female", "male"),
  symbol = c("\u2640", "\u2642"),
  body_size  = c(0.9, 0.9),
  organ_size = c(1.10, 0.95)
)

# Define theme
base_theme <- theme_bw(base_size = 12) +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    legend.position  = "none",
    axis.text.x      = element_blank(),
    axis.text.y      = element_blank(),
    axis.ticks.x     = element_blank(),
    axis.ticks.y     = element_blank(),
    axis.title       = element_text(size = 11),
    panel.border     = element_rect(colour = "black", fill = NA, size = 0.8),
    plot.margin      = margin(3, 3, 3, 3)
  )

# helper conceptual tribble to identify data structure
lines_df <- tribble(
  ~ panel, ~ sex,      ~ intercept, ~ slope,
  "a",    "male",        0.25,      0.18*3,
  "a",    "female",      0.60,      0.06*3,
  "b",    "male",        0.25,      0.12*3,
  "b",    "female",      0.55,      0.11*3,
  "c",    "male",        0.35,      0.18*3,
  "c",    "female",      0.38,      0.06*3,
  "d",    "male",        0.35,      0.11*3,
  "d",    "female",      0.36,      0.10*3
)

# position of vertical line
x_mid <- 0.5

# Adjust only panel c so that male and female intersect at x_mid...I am doing this for illustrative purposes only.
lines_df2 <- lines_df %>%
  group_by(panel) %>%
  mutate(
    intercept = case_when(
      panel != "c" ~ intercept,
      panel == "c" & sex == "male" ~ intercept,
      panel == "c" & sex == "female" ~ (
        intercept[sex == "male"] +
          slope[sex == "male"] * x_mid -
          slope[sex == "female"] * x_mid
      )
    )
  ) %>%
  ungroup()

x_seq <- seq(0, 1, length.out = 50)

plot_df <- lines_df2 %>%
  tidyr::expand_grid(body_size = x_seq) %>%
  mutate(organ_size = intercept + slope * body_size)

# Function to add letters in the panel
add_tag_inside <- function(p, tag) {
  p +
    annotate(
      "text",
      x = -Inf, y = Inf,
      label = tag,
      hjust = -0.3, vjust = 1.3,
      size = 3.5
    )
}

# Function to make four scenarios of scaling for females and males
make_panel <- function(p, x_lab, y_lab = NULL) {
  ggplot(filter(plot_df, panel == p),
         aes(x = body_size, y = organ_size,
             colour = sex, group = sex)) +
    geom_vline(xintercept = x_mid,
               linetype = "dashed",
               colour   = "grey50",
               linewidth = 0.5) +
    geom_line(linewidth = 1.1) +
    geom_text(
      data = sex_symbols,
      aes(x = body_size, y = organ_size,
          label = symbol, colour = sex),
      inherit.aes = FALSE,
      size = 4.5,
      fontface = "bold"
    ) +
    scale_colour_manual(values = cols, name = NULL) +
    labs(x = x_lab, y = y_lab) +
    coord_cartesian(xlim = c(0, 1), ylim = c(0, 1.2), expand = FALSE) +
    base_theme
}
# Make plots 
p1 <- make_panel("a", x_lab = "Body size", y_lab = "Organ size") %>% add_tag_inside("A")
p2 <- make_panel("b", x_lab = "Body size", y_lab = NULL)         %>% add_tag_inside("B")
p3 <- make_panel("c", x_lab = "Body size", y_lab = NULL)         %>% add_tag_inside("C")
p4 <- make_panel("d", x_lab = "Body size", y_lab = NULL)         %>% add_tag_inside("D")

# Simulations of posterior distributions (again 4 illustrative only) 
set.seed(6955) 

post_df <- lines_df2 %>%
  group_by(panel, sex) %>%
  do({
    tibble(
      draw = 1:1000,
      intercept = rnorm(1000, .$intercept, 0.05),
      slope     = rnorm(1000, .$slope,     0.02)
    )
  }) %>%
  ungroup()

post_long <- post_df %>%
  pivot_longer(c(intercept, slope),
               names_to = "parameter",
               values_to = "value")

make_posterior_panel <- function(panel_id, param, x_lab, y_lab = NULL) {
  post_long %>%
    filter(panel == panel_id, parameter == param) %>%
    ggplot(aes(x = value, fill = sex, colour = sex)) +
    geom_density(alpha = 0.6, adjust = 0.8, color = "white") +
    scale_fill_manual(values = cols, name = NULL) +
    scale_colour_manual(values = cols, guide = "none") +
    labs(x = x_lab, y = y_lab) +
    base_theme
}

# Posterior intercepts (e–h)
p5  <- make_posterior_panel("a", "intercept", "Posterior intercept", y_lab = "Density") %>% add_tag_inside("E")
p6  <- make_posterior_panel("b", "intercept", "Posterior intercept", y_lab = NULL)      %>% add_tag_inside("F")
p7  <- make_posterior_panel("c", "intercept", "Posterior intercept", y_lab = NULL)      %>% add_tag_inside("G")
p8  <- make_posterior_panel("d", "intercept", "Posterior intercept", y_lab = NULL)      %>% add_tag_inside("H")

# Posterior slopes (i–l)
p9  <- make_posterior_panel("a", "slope", "Posterior slope", y_lab = "Density") %>% add_tag_inside("I")
p10 <- make_posterior_panel("b", "slope", "Posterior slope", y_lab = NULL)      %>% add_tag_inside("J")
p11 <- make_posterior_panel("c", "slope", "Posterior slope", y_lab = NULL)      %>% add_tag_inside("K")
p12 <- make_posterior_panel("d", "slope", "Posterior slope", y_lab = NULL)      %>% add_tag_inside("L")

# Combine plots
Figure_1 <- wrap_plots(
  p1, p2, p3, p4,
  p5, p6, p7, p8,
  p9, p10, p11, p12,
  nrow = 3, ncol = 4, byrow = TRUE
) +
  plot_layout(
    axis_titles = "collect_x" 
  )

Figure_1

# Export figure
ggsave("../outputs/Figure_1.png", Figure_1, width = 7, height = 6, dpi = 1200)

Figure 2

# Prepare data to merge with tree
dat_mean$organ_grouped <- gsub("_", " ", dat_mean$organ_grouped)

# I make a yes/no binary variable to indicate in the tree which organs were measured for a given species
summary_data <- dat_mean %>%
  count(phylo, organ_grouped) %>%
  mutate(presence = if_else(n > 0, "Yes", "No")) %>%
  select(-n) %>%
  pivot_wider(
    names_from = organ_grouped,
    values_from = presence,
    values_fill = list(presence = "No")
  )

summary_data <- summary_data[summary_data$phylo %in% tree$tip.label, ]

datF_organs <- summary_data %>%
  column_to_rownames("phylo") %>%
  as.matrix()

# force matc between data and tree
datF_organs <- datF_organs[tree$tip.label, , drop = FALSE]

# Create phylogeny tree
circ <- ggtree(tree, layout = "fan", open.angle = 18) 
circ <- rotate_tree(circ, 90)

Figure_2 <- gheatmap(
  circ,
  datF_organs,
  width = 0.4, 
  offset = 0, 
  colnames_offset_x = 0, 
  colnames_offset_y = 0, 
  font.size = 2,
  hjust = 0
) +
  scale_fill_manual(values = c("grey90", "grey20"), name = "Organ measured")  +
  theme(
    legend.position = c(0.57, 0.55), 
    legend.background = element_rect(fill = "transparent", colour = NA),
    legend.box.background = element_rect(fill = "transparent", colour = NA),
    legend.title = element_text(size = 8),      
    legend.text = element_text(size = 6),       
    legend.key.size = unit(0.8, "cm"),         
    legend.key.height = unit(0.4, "cm"),       
    legend.key.width = unit(0.4, "cm")         
  )

Figure_2

# Export figure
ggsave('../outputs/Figure_2.pdf', Figure_2, width = 7, height = 7)

Analysis by organ

Liver

Filter data and prune tree

# Filter organ of interest
dat_organ <- dat_mean %>%
  group_by(organ_grouped, sex) %>%
  ungroup() %>%
  filter(organ_grouped == "liver")

# Check for database/tree species mismatches
setdiff(tree$tip.label, dat_organ$phylo)      # Tree species absent in database

  [1] "Petrochromis_famula"             "Lobochilotes_labiatus"          
  [3] "Simochromis_diagramma"           "Petrochromis_orthognathus"      
  [5] "Ctenochromis_horei"              "Limnotilapia_dardennii"         
  [7] "Gnathochromis_pfefferi"          "Tropheus_moorii"                
  [9] "Tropheus_brichardi"              "Haplotaxodon_microlepis"        
 [11] "Perissodus_microlepis"           "Benthochromis_tricoti"          
 [13] "Cyprichromis_leptosoma"          "Greenwoodochromis_christyi"     
 [15] "Limnochromis_staneri"            "Gnathochromis_permaxillaris"    
 [17] "Baileychromis_centropomoides"    "Triglachromis_otostigma"        
 [19] "Cyphotilapia_frontosa"           "Ophthalmotilapia_ventralis"     
 [21] "Cyathopharynx_furcifer"          "Ophthalmotilapia_nasuta"        
 [23] "Ophthalmotilapia_boops"          "Aulonocranus_dewindti"          
 [25] "Ectodus_descampsii"              "Callochromis_melanostigma"      
 [27] "Xenotilapia_melanogenys"         "Xenotilapia_flavipinnis"        
 [29] "Eretmodus_cyanostictus"          "Tanganicodus_irsacae"           
 [31] "Spathodus_marlieri"              "Spathodus_erythrodon"           
 [33] "Neolamprologus_brichardi"        "Julidochromis_regani"           
 [35] "Julidochromis_marlieri"          "Neolamprologus_tetracanthus"    
 [37] "Telmatochromis_temporalis"       "Lepidiolamprologus_elongatus"   
 [39] "Lepidiolamprologus_profundicola" "Lepidiolamprologus_nkambae"     
 [41] "Altolamprologus_compressiceps"   "Neolamprologus_brevis"          
 [43] "Chalinochromis_brichardi"        "Julidochromis_ornatus"          
 [45] "Lamprologus_ornatipinnis"        "Neolamprologus_pulcher"         
 [47] "Variabilichromis_moorii"         "Neolamprologus_tretocephalus"   
 [49] "Neolamprologus_sexfasciatus"     "Trematocara_unimaculatum"       
 [51] "Hemibates_stenosoma"             "Bathybates_ferox"               
 [53] "Poecilia_reticulata"             "Nerophis_lumbriciformis"        
 [55] "Syngnathus_abaster"              "Syngnathus_schlegeli"           
 [57] "Hippocampus_comes"               "Hippocampus_trimaculatus"       
 [59] "Hippocampus_kuda"                "Hippocampus_spinosissimus"      
 [61] "Hippocampus_abdominalis"         "Hippichthys_cyanospilos"        
 [63] "Syngnathoides_biaculeatus"       "Corythoichthys_intestinalis"    
 [65] "Corythoichthys_haematopterus"    "Doryichthys_boaja"              
 [67] "Doryrhamphus_japonicus"          "Entelurus_aequoreus"            
 [69] "Syncerus_caffer"                 "Tragelaphus_eurycerus"          
 [71] "Tragelaphus_scriptus"            "Redunca_arundinum"              
 [73] "Antidorcas_marsupialis"          "Nanger_granti"                  
 [75] "Raphicerus_campestris"           "Cephalophus_natalensis"         
 [77] "Sylvicapra_grimmia"              "Oreotragus_oreotragus"          
 [79] "Hippotragus_equinus"             "Addax_nasomaculatus"            
 [81] "Oryx_gazella"                    "Connochaetes_taurinus"          
 [83] "Damaliscus_lunatus"              "Ovis_aries"                     
 [85] "Aepyceros_melampus"              "Pteropus_alecto"                
 [87] "Pteropus_poliocephalus"          "Pteropus_scapulatus"            
 [89] "Chlorocebus_aethiops"            "Chlorocebus_sabaeus"            
 [91] "Erythrocebus_patas"              "Papio_cynocephalus"             
 [93] "Papio_anubis"                    "Cercocebus_atys"                
 [95] "Macaca_arctoides"                "Macaca_mulatta"                 
 [97] "Macaca_maura"                    "Trachypithecus_francoisi"       
 [99] "Phrynocephalus_vlangalii"        "Platycercus_caledonicus"        
[101] "Platycercus_elegans"             "Platycercus_eximius"            
[103] "Platycercus_venustus"            "Barnardius_zonarius"            
[105] "Northiella_haematogaster"        "Lathamus_discolor"              
[107] "Neophema_chrysostoma"            "Neophema_pulchella"             
[109] "Psittacula_krameri"              "Eclectus_roratus"               
[111] "Aprosmictus_jonquillaceus"       "Alisterus_scapularis"           
[113] "Alisterus_amboinensis"           "Polytelis_alexandrae"           
[115] "Trichoglossus_haematodus"        "Agapornis_lilianae"             
[117] "Agapornis_taranta"               "Loriculus_vernalis"             
[119] "Ara_ararauna"                    "Guaruba_guarouba"               
[121] "Enicognathus_leptorhynchus"      "Pionites_melanocephalus"        
[123] "Forpus_coelestis"                "Forpus_passerinus"              
[125] "Amazona_leucocephala"            "Amazona_albifrons"              
[127] "Amazona_pretrei"                 "Amazona_vinacea"                
[129] "Amazona_finschi"                 "Amazona_amazonica"              
[131] "Amazona_ochrocephala"            "Amazona_aestiva"                
[133] "Brotogeris_versicolurus"         "Brotogeris_pyrrhoptera"         
[135] "Bolborhynchus_lineola"           "Poicephalus_gulielmi"           
[137] "Poicephalus_meyeri"              "Psittacus_erithacus"            
[139] "Serinus_mennelli"                "Padda_oryzivora"                
[141] "Lonchura_punctulata"             "Lonchura_flaviprymna"           
[143] "Lonchura_bicolor"                "Hypargos_niveoguttatus"         
[145] "Pytilia_hypogrammica"            "Uraeginthus_bengalus"           
[147] "Amandava_amandava"               "Plectrophenax_nivalis"          
[149] "Serinus_canaria"                 "Serinus_serinus"                
[151] "Loxia_curvirostra"               "Carduelis_carduelis"            
[153] "Carpodacus_roseus"               "Uragus_sibiricus"               
[155] "Pyrrhula_pyrrhula"               "Pinicola_enucleator"            
[157] "Bucanetes_githagineus"           "Coccothraustes_coccothraustes"  
[159] "Mycerobas_affinis"               "Mycerobas_carnipes"             
[161] "Fringilla_montifringilla"        "Fringilla_coelebs"              
[163] "Anas_acuta"                      "Anas_crecca"                    
[165] "Mergus_serrator"                 "Mergus_merganser"               
[167] "Bucephala_clangula"              "Melanitta_nigra"                
[169] "Somateria_mollissima"            "Callonetta_leucophrys"          
[171] "Branta_bernicla"                 "Branta_leucopsis"               
[173] "Anser_anser"                     "Cygnus_columbianus"             
[175] "Tadorna_tadorna"                 "Pucrasia_macrolopha"            
[177] "Phasianus_colchicus"             "Chrysolophus_pictus"            
[179] "Lophura_nycthemera"              "Lophura_ignita"                 
[181] "Syrmaticus_ellioti"              "Syrmaticus_reevesii"            
[183] "Perdix_perdix"                   "Tetrao_urogallus"               
[185] "Lagopus_muta"                    "Lagopus_lagopus"                
[187] "Lophophorus_impejanus"           "Tragopan_temminckii"            
[189] "Pavo_cristatus"                  "Gallus_sonneratii"              
[191] "Rollulus_rouloul"                "Arborophila_torqueola"          

setdiff(dat_organ$phylo, tree$tip.label)      # Database species absent in tree

character(0)

# Exclude mismatched species for consistency
tree_organ <- keep.tip(tree, intersect(tree$tip.label, dat_organ$phylo))
dat_organ <- dat_organ[dat_organ$phylo %in% tree_organ$tip.label, ]

# Re-check for mismatches after pruning
setdiff(tree_organ$tip.label, dat_organ$phylo)

character(0)

setdiff(dat_organ$phylo, tree_organ$tip.label)

character(0)

# Is ultrametric?
is.ultrametric(tree_organ)

[1] TRUE

# how many species?
length(tree_organ$tip.label)

[1] 12

# Compute phylogenetic covariance matrix
A_organ <- ape::vcv.phylo(tree_organ, corr = TRUE)

Model simple

model_simple_liver <- brm(
  log10_organ_size ~ 
    1 + 
    log10_body_size +
    (1 + log10_body_size | sex),
  data = dat_organ,
  family = gaussian(),
    prior = c(
    prior(normal(0.75, 0.5), class = "b", coef = "log10_body_size"),
    prior(normal(0, 2), class = "Intercept"),
    prior(student_t(3, 0, 0.5), class = "sd"),
    prior(student_t(3, 0, 1),   class = "sigma")
  ),
  seed = 6955,
  sample_prior = TRUE, 
  chains = 4, 
  cores = 4, 
  threads = threading(2),
  iter = 4000, 
  warmup = 2000,
  control = list(max_treedepth = 15,
                 adapt_delta = 0.99),
  save_pars = save_pars(all = TRUE),
  silent = TRUE, refresh = 0
)

Running MCMC with 4 parallel chains, with 2 thread(s) per chain...

Chain 2 finished in 2.8 seconds.
Chain 4 finished in 3.5 seconds.
Chain 3 finished in 3.7 seconds.
Chain 1 finished in 4.4 seconds.

All 4 chains finished successfully.
Mean chain execution time: 3.6 seconds.
Total execution time: 4.5 seconds.

Model check

pp_check(model_simple_liver)

Model summary

summary(model_simple_liver)

 Family: gaussian 
  Links: mu = identity 
Formula: log10_organ_size ~ 1 + log10_body_size + (1 + log10_body_size | sex) 
   Data: dat_organ (Number of observations: 24) 
  Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
         total post-warmup draws = 8000

Multilevel Hyperparameters:
~sex (Number of levels: 2) 
                               Estimate Est.Error l-95% CI u-95% CI Rhat
sd(Intercept)                      0.30      0.31     0.01     1.12 1.00
sd(log10_body_size)                0.18      0.22     0.00     0.81 1.00
cor(Intercept,log10_body_size)    -0.08      0.59    -0.97     0.95 1.00
                               Bulk_ESS Tail_ESS
sd(Intercept)                      2690     2972
sd(log10_body_size)                1903     3263
cor(Intercept,log10_body_size)     4124     4348

Regression Coefficients:
                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept          -1.34      0.31    -1.98    -0.75 1.00     2946     2317
log10_body_size     0.85      0.15     0.48     1.14 1.00     1706     1586

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.18      0.03     0.13     0.26 1.00     5259     4010

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Coeficients by sex

# Extract combined coefficients (fixed + random)
coef(model_simple_liver)$sex %>%
{tibble(
  sex = dimnames(.)[[1]],
  intercept    = .[, "Estimate", "Intercept"],
  intercept_low  = .[, "Q2.5", "Intercept"], 
  intercept_high = .[, "Q97.5", "Intercept"]
)} %>%
bind_cols(
  tibble(
    sex = dimnames(coef(model_simple_liver)$sex)[[1]],
    slope = coef(model_simple_liver)$sex[, "Estimate", "log10_body_size"]
  ) %>% select(slope),
  tibble(
    sex = dimnames(coef(model_simple_liver)$sex)[[1]],
    slope_low  = coef(model_simple_liver)$sex[, "Q2.5", "log10_body_size"],
    slope_high = coef(model_simple_liver)$sex[, "Q97.5", "log10_body_size"]
  ) %>% select(slope_low, slope_high)
)

# A tibble: 2 × 7
  sex    intercept intercept_low intercept_high slope slope_low slope_high
  <chr>      <dbl>         <dbl>          <dbl> <dbl>     <dbl>      <dbl>
1 female     -1.29         -1.56          -1.01 0.868     0.756      0.976
2 male       -1.36         -1.63          -1.10 0.876     0.772      0.986

Model with phylogeny

model_phylo_liver <- brm(
  log10_organ_size ~ 
    1 + 
    log10_body_size +
    (1 + log10_body_size | sex) +
    (1 | gr(phylo, cov = A_organ)),
  data = dat_organ,
  data2 = list(A_organ = A_organ),
  family = gaussian(),
    prior = c(
    prior(normal(0.75, 0.5), class = "b", coef = "log10_body_size"),
    prior(normal(0, 2), class = "Intercept"),
    prior(student_t(3, 0, 0.5), class = "sd"),
    prior(student_t(3, 0, 1),   class = "sigma")
  ),
  seed = 6955,
  sample_prior = TRUE, 
  chains = 4, 
  cores = 4, 
  threads = threading(2),
  iter = 4000, 
  warmup = 2000,
  control = list(max_treedepth = 15,
                 adapt_delta = 0.99),
  save_pars = save_pars(all = TRUE),
  silent = TRUE, refresh = 0
)

Running MCMC with 4 parallel chains, with 2 thread(s) per chain...

Chain 1 finished in 13.6 seconds.
Chain 3 finished in 13.9 seconds.
Chain 4 finished in 14.4 seconds.
Chain 2 finished in 15.2 seconds.

All 4 chains finished successfully.
Mean chain execution time: 14.3 seconds.
Total execution time: 15.3 seconds.

Model check

pp_check(model_phylo_liver)

Model summary

summary(model_phylo_liver)

 Family: gaussian 
  Links: mu = identity 
Formula: log10_organ_size ~ 1 + log10_body_size + (1 + log10_body_size | sex) + (1 | gr(phylo, cov = A_organ)) 
   Data: dat_organ (Number of observations: 24) 
  Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
         total post-warmup draws = 8000

Multilevel Hyperparameters:
~phylo (Number of levels: 12) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.18      0.06     0.10     0.32 1.00     2763     4184

~sex (Number of levels: 2) 
                               Estimate Est.Error l-95% CI u-95% CI Rhat
sd(Intercept)                      0.25      0.28     0.01     0.99 1.00
sd(log10_body_size)                0.16      0.21     0.00     0.76 1.00
cor(Intercept,log10_body_size)    -0.04      0.62    -0.97     0.96 1.00
                               Bulk_ESS Tail_ESS
sd(Intercept)                      3416     3875
sd(log10_body_size)                2082     3720
cor(Intercept,log10_body_size)     5444     5035

Regression Coefficients:
                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept          -1.48      0.27    -2.04    -0.96 1.00     3236     3945
log10_body_size     0.88      0.14     0.51     1.13 1.00     1986     1545

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.07      0.01     0.04     0.10 1.00     3602     5314

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Phylogenetic signal

hyp <- paste(
  "sd_phylo__Intercept^2 /",
  "(sd_phylo__Intercept^2 + sd_sex__Intercept^2 + sigma^2) = 0"
)
(hyp <- hypothesis(model_phylo_liver, hyp, class = NULL))

Hypothesis Tests for class :
                Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
1 (sd_phylo__Interc... = 0     0.48       0.3     0.02     0.93       0.19
  Post.Prob Star
1      0.16    *
---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.

lambda_liver <- hyp$hypothesis$Estimate

Coeficients by sex

# Extract combined coefficients (fixed + random)
coef(model_phylo_liver)$sex %>%
{tibble(
  sex = dimnames(.)[[1]],
  intercept    = .[, "Estimate", "Intercept"],
  intercept_low  = .[, "Q2.5", "Intercept"], 
  intercept_high = .[, "Q97.5", "Intercept"]
)} %>%
bind_cols(
  tibble(
    sex = dimnames(coef(model_phylo_liver)$sex)[[1]],
    slope = coef(model_phylo_liver)$sex[, "Estimate", "log10_body_size"]
  ) %>% select(slope),
  tibble(
    sex = dimnames(coef(model_phylo_liver)$sex)[[1]],
    slope_low  = coef(model_phylo_liver)$sex[, "Q2.5", "log10_body_size"],
    slope_high = coef(model_phylo_liver)$sex[, "Q97.5", "log10_body_size"]
  ) %>% select(slope_low, slope_high)
)

# A tibble: 2 × 7
  sex    intercept intercept_low intercept_high slope slope_low slope_high
  <chr>      <dbl>         <dbl>          <dbl> <dbl>     <dbl>      <dbl>
1 female     -1.45         -1.73          -1.19 0.896     0.827      0.969
2 male       -1.50         -1.78          -1.24 0.893     0.825      0.966

Model comparison

loo_simple <- loo(model_simple_liver)
loo_phylo <- loo(model_phylo_liver)
loo_compare(loo_simple, loo_phylo)

                   elpd_diff se_diff
model_phylo_liver    0.0       0.0  
model_simple_liver -19.2       3.8  

In this comparison, the model that accounts for phylogeny performs better to the model that treats species as independent observational units. The phylogenetic model (model_phylo_liver) is used as the reference with elpd_diff = 0, whereas the simpler model (model_simple_liver) has an expected log predictive density (ELPD) that is about 19.6 units lower, with an associated standard error of 3.8. In the context of leave‑one‑out cross‑validation, differences exceeding roughly four to five times their standard error are typically regarded as strong evidence (Vehtari et al., 2016; Sivula et al., 2025). Here, 19.6/3.8≈5.2, which provides strong evidence in favour of the phylogenetic model.

Heart

Filter data and prune tree

# Filter organ of interest
dat_organ <- dat_mean %>%
  group_by(organ_grouped, sex) %>%
  ungroup() %>%
  filter(organ_grouped == "heart")

# Check for database/tree species mismatches
setdiff(tree$tip.label, dat_organ$phylo)      # Tree species absent in database

  [1] "Petrochromis_famula"             "Lobochilotes_labiatus"          
  [3] "Simochromis_diagramma"           "Petrochromis_orthognathus"      
  [5] "Ctenochromis_horei"              "Limnotilapia_dardennii"         
  [7] "Gnathochromis_pfefferi"          "Tropheus_moorii"                
  [9] "Tropheus_brichardi"              "Haplotaxodon_microlepis"        
 [11] "Perissodus_microlepis"           "Benthochromis_tricoti"          
 [13] "Cyprichromis_leptosoma"          "Greenwoodochromis_christyi"     
 [15] "Limnochromis_staneri"            "Gnathochromis_permaxillaris"    
 [17] "Baileychromis_centropomoides"    "Triglachromis_otostigma"        
 [19] "Cyphotilapia_frontosa"           "Ophthalmotilapia_ventralis"     
 [21] "Cyathopharynx_furcifer"          "Ophthalmotilapia_nasuta"        
 [23] "Ophthalmotilapia_boops"          "Aulonocranus_dewindti"          
 [25] "Ectodus_descampsii"              "Callochromis_melanostigma"      
 [27] "Xenotilapia_melanogenys"         "Xenotilapia_flavipinnis"        
 [29] "Eretmodus_cyanostictus"          "Tanganicodus_irsacae"           
 [31] "Spathodus_marlieri"              "Spathodus_erythrodon"           
 [33] "Neolamprologus_brichardi"        "Julidochromis_regani"           
 [35] "Julidochromis_marlieri"          "Neolamprologus_tetracanthus"    
 [37] "Telmatochromis_temporalis"       "Lepidiolamprologus_elongatus"   
 [39] "Lepidiolamprologus_profundicola" "Lepidiolamprologus_nkambae"     
 [41] "Altolamprologus_compressiceps"   "Neolamprologus_brevis"          
 [43] "Chalinochromis_brichardi"        "Julidochromis_ornatus"          
 [45] "Lamprologus_ornatipinnis"        "Neolamprologus_pulcher"         
 [47] "Variabilichromis_moorii"         "Neolamprologus_tretocephalus"   
 [49] "Neolamprologus_sexfasciatus"     "Trematocara_unimaculatum"       
 [51] "Hemibates_stenosoma"             "Bathybates_ferox"               
 [53] "Oreochromis_niloticus"           "Nerophis_lumbriciformis"        
 [55] "Syngnathus_abaster"              "Syngnathus_schlegeli"           
 [57] "Hippocampus_comes"               "Hippocampus_trimaculatus"       
 [59] "Hippocampus_kuda"                "Hippocampus_spinosissimus"      
 [61] "Hippocampus_abdominalis"         "Hippichthys_cyanospilos"        
 [63] "Syngnathoides_biaculeatus"       "Corythoichthys_intestinalis"    
 [65] "Corythoichthys_haematopterus"    "Doryichthys_boaja"              
 [67] "Doryrhamphus_japonicus"          "Entelurus_aequoreus"            
 [69] "Anguilla_anguilla"               "Syncerus_caffer"                
 [71] "Tragelaphus_eurycerus"           "Tragelaphus_scriptus"           
 [73] "Redunca_arundinum"               "Antidorcas_marsupialis"         
 [75] "Nanger_granti"                   "Raphicerus_campestris"          
 [77] "Cephalophus_natalensis"          "Sylvicapra_grimmia"             
 [79] "Oreotragus_oreotragus"           "Hippotragus_equinus"            
 [81] "Addax_nasomaculatus"             "Oryx_gazella"                   
 [83] "Connochaetes_taurinus"           "Damaliscus_lunatus"             
 [85] "Ovis_aries"                      "Aepyceros_melampus"             
 [87] "Pteropus_alecto"                 "Pteropus_poliocephalus"         
 [89] "Pteropus_scapulatus"             "Chlorocebus_aethiops"           
 [91] "Chlorocebus_sabaeus"             "Erythrocebus_patas"             
 [93] "Papio_cynocephalus"              "Papio_anubis"                   
 [95] "Cercocebus_atys"                 "Macaca_arctoides"               
 [97] "Macaca_mulatta"                  "Macaca_maura"                   
 [99] "Trachypithecus_francoisi"        "Platycercus_caledonicus"        
[101] "Platycercus_elegans"             "Platycercus_eximius"            
[103] "Platycercus_venustus"            "Barnardius_zonarius"            
[105] "Northiella_haematogaster"        "Lathamus_discolor"              
[107] "Neophema_chrysostoma"            "Neophema_pulchella"             
[109] "Psittacula_krameri"              "Eclectus_roratus"               
[111] "Aprosmictus_jonquillaceus"       "Alisterus_scapularis"           
[113] "Alisterus_amboinensis"           "Polytelis_alexandrae"           
[115] "Trichoglossus_haematodus"        "Agapornis_lilianae"             
[117] "Agapornis_taranta"               "Loriculus_vernalis"             
[119] "Ara_ararauna"                    "Guaruba_guarouba"               
[121] "Enicognathus_leptorhynchus"      "Pionites_melanocephalus"        
[123] "Forpus_coelestis"                "Forpus_passerinus"              
[125] "Amazona_leucocephala"            "Amazona_albifrons"              
[127] "Amazona_pretrei"                 "Amazona_vinacea"                
[129] "Amazona_finschi"                 "Amazona_amazonica"              
[131] "Amazona_ochrocephala"            "Amazona_aestiva"                
[133] "Brotogeris_versicolurus"         "Brotogeris_pyrrhoptera"         
[135] "Bolborhynchus_lineola"           "Poicephalus_gulielmi"           
[137] "Poicephalus_meyeri"              "Psittacus_erithacus"            
[139] "Serinus_mennelli"                "Ficedula_hypoleuca"             
[141] "Padda_oryzivora"                 "Lonchura_punctulata"            
[143] "Lonchura_flaviprymna"            "Lonchura_bicolor"               
[145] "Hypargos_niveoguttatus"          "Pytilia_hypogrammica"           
[147] "Uraeginthus_bengalus"            "Amandava_amandava"              
[149] "Plectrophenax_nivalis"           "Serinus_canaria"                
[151] "Serinus_serinus"                 "Loxia_curvirostra"              
[153] "Carduelis_carduelis"             "Carpodacus_roseus"              
[155] "Uragus_sibiricus"                "Pyrrhula_pyrrhula"              
[157] "Pinicola_enucleator"             "Bucanetes_githagineus"          
[159] "Coccothraustes_coccothraustes"   "Mycerobas_affinis"              
[161] "Mycerobas_carnipes"              "Fringilla_montifringilla"       
[163] "Fringilla_coelebs"               "Anas_acuta"                     
[165] "Anas_crecca"                     "Mergus_serrator"                
[167] "Mergus_merganser"                "Bucephala_clangula"             
[169] "Melanitta_nigra"                 "Somateria_mollissima"           
[171] "Callonetta_leucophrys"           "Branta_bernicla"                
[173] "Branta_leucopsis"                "Anser_anser"                    
[175] "Cygnus_columbianus"              "Tadorna_tadorna"                
[177] "Pucrasia_macrolopha"             "Phasianus_colchicus"            
[179] "Chrysolophus_pictus"             "Lophura_nycthemera"             
[181] "Lophura_ignita"                  "Syrmaticus_ellioti"             
[183] "Syrmaticus_reevesii"             "Perdix_perdix"                  
[185] "Tetrao_urogallus"                "Lagopus_muta"                   
[187] "Lagopus_lagopus"                 "Lophophorus_impejanus"          
[189] "Tragopan_temminckii"             "Pavo_cristatus"                 
[191] "Gallus_gallus"                   "Gallus_sonneratii"              
[193] "Rollulus_rouloul"                "Arborophila_torqueola"          

setdiff(dat_organ$phylo, tree$tip.label)      # Database species absent in tree

character(0)

# Exclude mismatched species for consistency
tree_organ <- keep.tip(tree, intersect(tree$tip.label, dat_organ$phylo))
dat_organ <- dat_organ[dat_organ$phylo %in% tree_organ$tip.label, ]

# Re-check for mismatches after pruning
setdiff(tree_organ$tip.label, dat_organ$phylo)

character(0)

setdiff(dat_organ$phylo, tree_organ$tip.label)

character(0)

# Is ultrametric?
is.ultrametric(tree_organ)

[1] TRUE

# how many species?
length(tree_organ$tip.label)

[1] 10

# Compute phylogenetic covariance matrix
A_organ <- ape::vcv.phylo(tree_organ, corr = TRUE)

Model simple

model_simple_heart <- brm(
  log10_organ_size ~ 
    1 + 
    log10_body_size +
    (1 + log10_body_size | sex),
  data = dat_organ,
  family = gaussian(),
    prior = c(
    prior(normal(0.75, 0.5), class = "b", coef = "log10_body_size"),
    prior(normal(0, 2), class = "Intercept"),
    prior(student_t(3, 0, 0.5), class = "sd"),
    prior(student_t(3, 0, 1),   class = "sigma")
  ),
  seed = 6955,
  sample_prior = TRUE, 
  chains = 4, 
  cores = 4, 
  threads = threading(2),
  iter = 4000, 
  warmup = 2000,
  control = list(max_treedepth = 15,
                 adapt_delta = 0.99),
  save_pars = save_pars(all = TRUE),
  silent = TRUE, refresh = 0
)

Running MCMC with 4 parallel chains, with 2 thread(s) per chain...

Chain 2 finished in 2.0 seconds.
Chain 4 finished in 2.6 seconds.
Chain 1 finished in 3.2 seconds.
Chain 3 finished in 3.2 seconds.

All 4 chains finished successfully.
Mean chain execution time: 2.7 seconds.
Total execution time: 3.3 seconds.

Model check

pp_check(model_simple_heart)

Model summary

summary(model_simple_heart)

 Family: gaussian 
  Links: mu = identity 
Formula: log10_organ_size ~ 1 + log10_body_size + (1 + log10_body_size | sex) 
   Data: dat_organ (Number of observations: 20) 
  Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
         total post-warmup draws = 8000

Multilevel Hyperparameters:
~sex (Number of levels: 2) 
                               Estimate Est.Error l-95% CI u-95% CI Rhat
sd(Intercept)                      0.27      0.30     0.01     1.04 1.00
sd(log10_body_size)                0.19      0.22     0.00     0.80 1.00
cor(Intercept,log10_body_size)    -0.07      0.61    -0.97     0.95 1.00
                               Bulk_ESS Tail_ESS
sd(Intercept)                      3635     3605
sd(log10_body_size)                2149     3479
cor(Intercept,log10_body_size)     4205     4565

Regression Coefficients:
                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept          -2.92      0.29    -3.49    -2.35 1.00     3737     3845
log10_body_size     1.20      0.16     0.79     1.49 1.00     2413     2315

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.25      0.05     0.18     0.36 1.00     5488     3595

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Coefficients by sex

# Extract combined coefficients (fixed + random)
coef(model_simple_heart)$sex %>%
{tibble(
  sex = dimnames(.)[[1]],
  intercept    = .[, "Estimate", "Intercept"],
  intercept_low  = .[, "Q2.5", "Intercept"], 
  intercept_high = .[, "Q97.5", "Intercept"]
)} %>%
bind_cols(
  tibble(
    sex = dimnames(coef(model_simple_heart)$sex)[[1]],
    slope = coef(model_simple_heart)$sex[, "Estimate", "log10_body_size"]
  ) %>% select(slope),
  tibble(
    sex = dimnames(coef(model_simple_heart)$sex)[[1]],
    slope_low  = coef(model_simple_heart)$sex[, "Q2.5", "log10_body_size"],
    slope_high = coef(model_simple_heart)$sex[, "Q97.5", "log10_body_size"]
  ) %>% select(slope_low, slope_high)
)

# A tibble: 2 × 7
  sex    intercept intercept_low intercept_high slope slope_low slope_high
  <chr>      <dbl>         <dbl>          <dbl> <dbl>     <dbl>      <dbl>
1 female     -2.94         -3.21          -2.67  1.24      1.13       1.35
2 male       -2.91         -3.16          -2.65  1.23      1.13       1.34

Model with phylogeny

model_phylo_heart <- brm(
  log10_organ_size ~ 
    1 + 
    log10_body_size +
    (1 + log10_body_size | sex) +
    (1 | gr(phylo, cov = A_organ)),
  data = dat_organ,
  data2 = list(A_organ = A_organ),
  family = gaussian(),
    prior = c(
    prior(normal(0.75, 0.5), class = "b", coef = "log10_body_size"),
    prior(normal(0, 2), class = "Intercept"),
    prior(student_t(3, 0, 0.5), class = "sd"),
    prior(student_t(3, 0, 1),   class = "sigma")
  ),
  seed = 6955,
  sample_prior = TRUE, 
  chains = 4, 
  cores = 4, 
  threads = threading(2),
  iter = 4000, 
  warmup = 2000,
  control = list(max_treedepth = 15,
                 adapt_delta = 0.99),
  save_pars = save_pars(all = TRUE),
  silent = TRUE, refresh = 0
)

Running MCMC with 4 parallel chains, with 2 thread(s) per chain...

Chain 1 finished in 14.4 seconds.
Chain 2 finished in 14.6 seconds.
Chain 3 finished in 14.7 seconds.
Chain 4 finished in 17.3 seconds.

All 4 chains finished successfully.
Mean chain execution time: 15.2 seconds.
Total execution time: 17.4 seconds.

Model check

pp_check(model_phylo_heart)

Model summary

summary(model_phylo_heart)

 Family: gaussian 
  Links: mu = identity 
Formula: log10_organ_size ~ 1 + log10_body_size + (1 + log10_body_size | sex) + (1 | gr(phylo, cov = A_organ)) 
   Data: dat_organ (Number of observations: 20) 
  Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
         total post-warmup draws = 8000

Multilevel Hyperparameters:
~phylo (Number of levels: 10) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.53      0.16     0.31     0.91 1.00     1853     3477

~sex (Number of levels: 2) 
                               Estimate Est.Error l-95% CI u-95% CI Rhat
sd(Intercept)                      0.20      0.25     0.00     0.92 1.00
sd(log10_body_size)                0.13      0.18     0.00     0.65 1.00
cor(Intercept,log10_body_size)    -0.05      0.63    -0.98     0.97 1.00
                               Bulk_ESS Tail_ESS
sd(Intercept)                      3157     4410
sd(log10_body_size)                2478     4175
cor(Intercept,log10_body_size)     5539     5172

Regression Coefficients:
                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept          -2.63      0.39    -3.40    -1.84 1.00     2924     3641
log10_body_size     0.94      0.14     0.63     1.20 1.00     2455     2954

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.05      0.02     0.03     0.10 1.00     1641     2427

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Phylogenetic signal

hyp <- paste(
  "sd_phylo__Intercept^2 /",
  "(sd_phylo__Intercept^2 + sd_sex__Intercept^2 + sigma^2) = 0"
)
(hyp <- hypothesis(model_phylo_heart, hyp, class = NULL))

Hypothesis Tests for class :
                Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
1 (sd_phylo__Interc... = 0     0.84      0.22     0.21        1          0
  Post.Prob Star
1         0    *
---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.

lambda_heart <- hyp$hypothesis$Estimate

Coefficients by sex

# Extract combined coefficients (fixed + random)
coef(model_phylo_heart)$sex %>%
{tibble(
  sex = dimnames(.)[[1]],
  intercept    = .[, "Estimate", "Intercept"],
  intercept_low  = .[, "Q2.5", "Intercept"], 
  intercept_high = .[, "Q97.5", "Intercept"]
)} %>%
bind_cols(
  tibble(
    sex = dimnames(coef(model_phylo_heart)$sex)[[1]],
    slope = coef(model_phylo_heart)$sex[, "Estimate", "log10_body_size"]
  ) %>% select(slope),
  tibble(
    sex = dimnames(coef(model_phylo_heart)$sex)[[1]],
    slope_low  = coef(model_phylo_heart)$sex[, "Q2.5", "log10_body_size"],
    slope_high = coef(model_phylo_heart)$sex[, "Q97.5", "log10_body_size"]
  ) %>% select(slope_low, slope_high)
)

# A tibble: 2 × 7
  sex    intercept intercept_low intercept_high slope slope_low slope_high
  <chr>      <dbl>         <dbl>          <dbl> <dbl>     <dbl>      <dbl>
1 female     -2.64         -3.29          -1.94 0.948     0.764       1.12
2 male       -2.63         -3.28          -1.94 0.954     0.784       1.11

Model comparison

loo_simple <- loo(model_simple_heart)
loo_phylo <- loo(model_phylo_heart)
loo_compare(loo_simple, loo_phylo)

                   elpd_diff se_diff
model_phylo_heart    0.0       0.0  
model_simple_heart -25.4       3.5  

In this comparison, the model that accounts for phylogeny is clearly preferable to the model that treats species as independent observational units. The phylogenetic model (model_phylo_heart) is used as the reference with elpd_diff = 0, whereas the simpler model (model_simple_heart) has an expected log predictive density (ELPD) that is about 25.4 units lower, with an associated standard error of 3.5. In the context of leave‑one‑out cross‑validation, differences exceeding roughly four times their standard error are typically regarded as strong evidence. Here, 25.4/3.5≈7.3, indicating a very large difference relative to its uncertainty, and thus providing strong support for the phylogenetic model.

Pituitary glands

Filter data and prune tree

# Filter organ of interest
dat_organ <- dat_mean %>%
  group_by(organ_grouped, sex) %>%
  ungroup() %>%
  filter(organ_grouped == "pituitary glands")

# Check for database/tree species mismatches
setdiff(tree$tip.label, dat_organ$phylo)      # Tree species absent in database

  [1] "Petrochromis_famula"             "Lobochilotes_labiatus"          
  [3] "Simochromis_diagramma"           "Petrochromis_orthognathus"      
  [5] "Ctenochromis_horei"              "Limnotilapia_dardennii"         
  [7] "Gnathochromis_pfefferi"          "Tropheus_moorii"                
  [9] "Tropheus_brichardi"              "Haplotaxodon_microlepis"        
 [11] "Perissodus_microlepis"           "Benthochromis_tricoti"          
 [13] "Cyprichromis_leptosoma"          "Greenwoodochromis_christyi"     
 [15] "Limnochromis_staneri"            "Gnathochromis_permaxillaris"    
 [17] "Baileychromis_centropomoides"    "Triglachromis_otostigma"        
 [19] "Cyphotilapia_frontosa"           "Ophthalmotilapia_ventralis"     
 [21] "Cyathopharynx_furcifer"          "Ophthalmotilapia_nasuta"        
 [23] "Ophthalmotilapia_boops"          "Aulonocranus_dewindti"          
 [25] "Ectodus_descampsii"              "Callochromis_melanostigma"      
 [27] "Xenotilapia_melanogenys"         "Xenotilapia_flavipinnis"        
 [29] "Eretmodus_cyanostictus"          "Tanganicodus_irsacae"           
 [31] "Spathodus_marlieri"              "Spathodus_erythrodon"           
 [33] "Neolamprologus_brichardi"        "Julidochromis_regani"           
 [35] "Julidochromis_marlieri"          "Neolamprologus_tetracanthus"    
 [37] "Telmatochromis_temporalis"       "Lepidiolamprologus_elongatus"   
 [39] "Lepidiolamprologus_profundicola" "Lepidiolamprologus_nkambae"     
 [41] "Altolamprologus_compressiceps"   "Neolamprologus_brevis"          
 [43] "Chalinochromis_brichardi"        "Julidochromis_ornatus"          
 [45] "Lamprologus_ornatipinnis"        "Neolamprologus_pulcher"         
 [47] "Variabilichromis_moorii"         "Neolamprologus_tretocephalus"   
 [49] "Neolamprologus_sexfasciatus"     "Trematocara_unimaculatum"       
 [51] "Hemibates_stenosoma"             "Bathybates_ferox"               
 [53] "Oreochromis_niloticus"           "Poecilia_reticulata"            
 [55] "Nerophis_lumbriciformis"         "Syngnathus_abaster"             
 [57] "Syngnathus_schlegeli"            "Hippocampus_comes"              
 [59] "Hippocampus_trimaculatus"        "Hippocampus_kuda"               
 [61] "Hippocampus_spinosissimus"       "Hippocampus_abdominalis"        
 [63] "Hippichthys_cyanospilos"         "Syngnathoides_biaculeatus"      
 [65] "Corythoichthys_intestinalis"     "Corythoichthys_haematopterus"   
 [67] "Doryichthys_boaja"               "Doryrhamphus_japonicus"         
 [69] "Entelurus_aequoreus"             "Anguilla_anguilla"              
 [71] "Syncerus_caffer"                 "Tragelaphus_eurycerus"          
 [73] "Tragelaphus_scriptus"            "Redunca_arundinum"              
 [75] "Antidorcas_marsupialis"          "Nanger_granti"                  
 [77] "Raphicerus_campestris"           "Cephalophus_natalensis"         
 [79] "Sylvicapra_grimmia"              "Oreotragus_oreotragus"          
 [81] "Hippotragus_equinus"             "Addax_nasomaculatus"            
 [83] "Oryx_gazella"                    "Connochaetes_taurinus"          
 [85] "Damaliscus_lunatus"              "Ovis_aries"                     
 [87] "Aepyceros_melampus"              "Lasiopodomys_brandtii"          
 [89] "Mus_musculus"                    "Meriones_unguiculatus"          
 [91] "Chlorocebus_aethiops"            "Chlorocebus_sabaeus"            
 [93] "Erythrocebus_patas"              "Papio_cynocephalus"             
 [95] "Papio_anubis"                    "Cercocebus_atys"                
 [97] "Macaca_arctoides"                "Macaca_mulatta"                 
 [99] "Macaca_maura"                    "Trachypithecus_francoisi"       
[101] "Phrynocephalus_vlangalii"        "Platycercus_caledonicus"        
[103] "Platycercus_elegans"             "Platycercus_eximius"            
[105] "Platycercus_venustus"            "Barnardius_zonarius"            
[107] "Northiella_haematogaster"        "Lathamus_discolor"              
[109] "Neophema_chrysostoma"            "Neophema_pulchella"             
[111] "Psittacula_krameri"              "Eclectus_roratus"               
[113] "Aprosmictus_jonquillaceus"       "Alisterus_scapularis"           
[115] "Alisterus_amboinensis"           "Polytelis_alexandrae"           
[117] "Trichoglossus_haematodus"        "Agapornis_lilianae"             
[119] "Agapornis_taranta"               "Loriculus_vernalis"             
[121] "Ara_ararauna"                    "Guaruba_guarouba"               
[123] "Enicognathus_leptorhynchus"      "Pionites_melanocephalus"        
[125] "Forpus_coelestis"                "Forpus_passerinus"              
[127] "Amazona_leucocephala"            "Amazona_albifrons"              
[129] "Amazona_pretrei"                 "Amazona_vinacea"                
[131] "Amazona_finschi"                 "Amazona_amazonica"              
[133] "Amazona_ochrocephala"            "Amazona_aestiva"                
[135] "Brotogeris_versicolurus"         "Brotogeris_pyrrhoptera"         
[137] "Bolborhynchus_lineola"           "Poicephalus_gulielmi"           
[139] "Poicephalus_meyeri"              "Psittacus_erithacus"            
[141] "Serinus_mennelli"                "Ficedula_hypoleuca"             
[143] "Padda_oryzivora"                 "Lonchura_punctulata"            
[145] "Lonchura_flaviprymna"            "Lonchura_bicolor"               
[147] "Hypargos_niveoguttatus"          "Pytilia_hypogrammica"           
[149] "Uraeginthus_bengalus"            "Amandava_amandava"              
[151] "Plectrophenax_nivalis"           "Serinus_canaria"                
[153] "Serinus_serinus"                 "Loxia_curvirostra"              
[155] "Carduelis_carduelis"             "Carpodacus_roseus"              
[157] "Uragus_sibiricus"                "Pyrrhula_pyrrhula"              
[159] "Pinicola_enucleator"             "Bucanetes_githagineus"          
[161] "Coccothraustes_coccothraustes"   "Mycerobas_affinis"              
[163] "Mycerobas_carnipes"              "Fringilla_montifringilla"       
[165] "Fringilla_coelebs"               "Anas_platyrhynchos"             
[167] "Anas_acuta"                      "Anas_crecca"                    
[169] "Mergus_serrator"                 "Mergus_merganser"               
[171] "Bucephala_clangula"              "Melanitta_nigra"                
[173] "Somateria_mollissima"            "Callonetta_leucophrys"          
[175] "Branta_bernicla"                 "Branta_leucopsis"               
[177] "Anser_anser"                     "Cygnus_columbianus"             
[179] "Tadorna_tadorna"                 "Coturnix_japonica"              
[181] "Pucrasia_macrolopha"             "Phasianus_colchicus"            
[183] "Chrysolophus_pictus"             "Lophura_nycthemera"             
[185] "Lophura_ignita"                  "Syrmaticus_ellioti"             
[187] "Syrmaticus_reevesii"             "Perdix_perdix"                  
[189] "Tetrao_urogallus"                "Lagopus_muta"                   
[191] "Lagopus_lagopus"                 "Lophophorus_impejanus"          
[193] "Tragopan_temminckii"             "Pavo_cristatus"                 
[195] "Gallus_gallus"                   "Gallus_sonneratii"              
[197] "Rollulus_rouloul"                "Arborophila_torqueola"          
[199] "Numida_meleagris"               

setdiff(dat_organ$phylo, tree$tip.label)      # Database species absent in tree

character(0)

# Exclude mismatched species for consistency
tree_organ <- keep.tip(tree, intersect(tree$tip.label, dat_organ$phylo))
dat_organ <- dat_organ[dat_organ$phylo %in% tree_organ$tip.label, ]

# Re-check for mismatches after pruning
setdiff(tree_organ$tip.label, dat_organ$phylo)

character(0)

setdiff(dat_organ$phylo, tree_organ$tip.label)

character(0)

# Is ultrametric?
is.ultrametric(tree_organ)

[1] TRUE

# how many species?
length(tree_organ$tip.label)

[1] 5

# Compute phylogenetic covariance matrix
A_organ <- ape::vcv.phylo(tree_organ, corr = TRUE)

Model simple

model_simple_pituitary_glands <- brm(
  log10_organ_size ~ 
    1 + 
    log10_body_size +
    (1 + log10_body_size | sex),
  data = dat_organ,
  family = gaussian(),
    prior = c(
    prior(normal(0.75, 0.5), class = "b", coef = "log10_body_size"),
    prior(normal(0, 2), class = "Intercept"),
    prior(student_t(3, 0, 0.5), class = "sd"),
    prior(student_t(3, 0, 1),   class = "sigma")
  ),
  seed = 6955,
  sample_prior = TRUE, 
  chains = 4, 
  cores = 4, 
  threads = threading(2),
  iter = 4000, 
  warmup = 2000,
  control = list(max_treedepth = 15,
                 adapt_delta = 0.99),
  save_pars = save_pars(all = TRUE),
  silent = TRUE, refresh = 0
)

Running MCMC with 4 parallel chains, with 2 thread(s) per chain...

Chain 4 finished in 3.9 seconds.
Chain 3 finished in 4.1 seconds.
Chain 2 finished in 4.6 seconds.
Chain 1 finished in 4.7 seconds.

All 4 chains finished successfully.
Mean chain execution time: 4.3 seconds.
Total execution time: 4.8 seconds.

Model check

pp_check(model_simple_pituitary_glands)

Model summary

summary(model_simple_pituitary_glands)

 Family: gaussian 
  Links: mu = identity 
Formula: log10_organ_size ~ 1 + log10_body_size + (1 + log10_body_size | sex) 
   Data: dat_organ (Number of observations: 10) 
  Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
         total post-warmup draws = 8000

Multilevel Hyperparameters:
~sex (Number of levels: 2) 
                               Estimate Est.Error l-95% CI u-95% CI Rhat
sd(Intercept)                      0.31      0.32     0.01     1.19 1.00
sd(log10_body_size)                0.18      0.20     0.00     0.75 1.00
cor(Intercept,log10_body_size)    -0.07      0.60    -0.97     0.95 1.00
                               Bulk_ESS Tail_ESS
sd(Intercept)                      4018     3598
sd(log10_body_size)                1788     3113
cor(Intercept,log10_body_size)     4215     4178

Regression Coefficients:
                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept          -3.74      0.34    -4.36    -2.98 1.00     4347     4113
log10_body_size     0.75      0.16     0.46     1.12 1.00     2277     1588

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.10      0.03     0.05     0.18 1.00     4141     4111

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Coeficients by sex

# Extract combined coefficients (fixed + random)
coef(model_simple_pituitary_glands)$sex %>%
{tibble(
  sex = dimnames(.)[[1]],
  intercept    = .[, "Estimate", "Intercept"],
  intercept_low  = .[, "Q2.5", "Intercept"], 
  intercept_high = .[, "Q97.5", "Intercept"]
)} %>%
bind_cols(
  tibble(
    sex = dimnames(coef(model_simple_pituitary_glands)$sex)[[1]],
    slope = coef(model_simple_pituitary_glands)$sex[, "Estimate", "log10_body_size"]
  ) %>% select(slope),
  tibble(
    sex = dimnames(coef(model_simple_pituitary_glands)$sex)[[1]],
    slope_low  = coef(model_simple_pituitary_glands)$sex[, "Q2.5", "log10_body_size"],
    slope_high = coef(model_simple_pituitary_glands)$sex[, "Q97.5", "log10_body_size"]
  ) %>% select(slope_low, slope_high)
)

# A tibble: 2 × 7
  sex    intercept intercept_low intercept_high slope slope_low slope_high
  <chr>      <dbl>         <dbl>          <dbl> <dbl>     <dbl>      <dbl>
1 female     -3.75         -4.20          -3.31 0.731     0.571      0.892
2 male       -3.79         -4.21          -3.36 0.714     0.567      0.860

Model with phylogeny

model_phylo_pituitary_glands <- brm(
  log10_organ_size ~ 
    1 + 
    log10_body_size +
    (1 + log10_body_size | sex) +
    (1 | gr(phylo, cov = A_organ)),
  data = dat_organ,
  data2 = list(A_organ = A_organ),
  family = gaussian(),
    prior = c(
    prior(normal(0.75, 0.5), class = "b", coef = "log10_body_size"),
    prior(normal(0, 2), class = "Intercept"),
    prior(student_t(3, 0, 0.5), class = "sd"),
    prior(student_t(3, 0, 1),   class = "sigma")
  ),
  seed = 6955,
  sample_prior = TRUE, 
  chains = 4, 
  cores = 4, 
  threads = threading(2),
  iter = 4000, 
  warmup = 2000,
  control = list(max_treedepth = 15,
                 adapt_delta = 0.99),
  save_pars = save_pars(all = TRUE),
  silent = TRUE, refresh = 0
)

Running MCMC with 4 parallel chains, with 2 thread(s) per chain...

Chain 2 finished in 32.6 seconds.
Chain 4 finished in 42.7 seconds.
Chain 1 finished in 53.5 seconds.
Chain 3 finished in 95.4 seconds.

All 4 chains finished successfully.
Mean chain execution time: 56.0 seconds.
Total execution time: 95.5 seconds.

Model check

pp_check(model_phylo_pituitary_glands)

Model summary

summary(model_phylo_pituitary_glands)

 Family: gaussian 
  Links: mu = identity 
Formula: log10_organ_size ~ 1 + log10_body_size + (1 + log10_body_size | sex) + (1 | gr(phylo, cov = A_organ)) 
   Data: dat_organ (Number of observations: 10) 
  Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
         total post-warmup draws = 8000

Multilevel Hyperparameters:
~phylo (Number of levels: 5) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.16      0.10     0.05     0.40 1.00     2063     2185

~sex (Number of levels: 2) 
                               Estimate Est.Error l-95% CI u-95% CI Rhat
sd(Intercept)                      0.24      0.29     0.01     1.04 1.00
sd(log10_body_size)                0.18      0.20     0.00     0.73 1.00
cor(Intercept,log10_body_size)    -0.05      0.61    -0.98     0.96 1.00
                               Bulk_ESS Tail_ESS
sd(Intercept)                      4165     4489
sd(log10_body_size)                2562     3028
cor(Intercept,log10_body_size)     5442     4775

Regression Coefficients:
                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept          -3.83      0.45    -4.70    -2.92 1.00     2667     3522
log10_body_size     0.78      0.18     0.46     1.19 1.00     2415     2626

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.03      0.02     0.01     0.07 1.00      768     1777

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Phylogenetic signal

hyp <- paste(
  "sd_phylo__Intercept^2 /",
  "(sd_phylo__Intercept^2 + sd_sex__Intercept^2 + sigma^2) = 0"
)
(hyp <- hypothesis(model_phylo_pituitary_glands, hyp, class = NULL))

Hypothesis Tests for class :
                Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
1 (sd_phylo__Interc... = 0     0.49      0.35     0.01     0.99       0.32
  Post.Prob Star
1      0.24    *
---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.

lambda_pituitary_glands <- hyp$hypothesis$Estimate

Coeficients by sex

# Extract combined coefficients (fixed + random)
coef(model_phylo_pituitary_glands)$sex %>%
{tibble(
  sex = dimnames(.)[[1]],
  intercept    = .[, "Estimate", "Intercept"],
  intercept_low  = .[, "Q2.5", "Intercept"], 
  intercept_high = .[, "Q97.5", "Intercept"]
)} %>%
bind_cols(
  tibble(
    sex = dimnames(coef(model_phylo_pituitary_glands)$sex)[[1]],
    slope = coef(model_phylo_pituitary_glands)$sex[, "Estimate", "log10_body_size"]
  ) %>% select(slope),
  tibble(
    sex = dimnames(coef(model_phylo_pituitary_glands)$sex)[[1]],
    slope_low  = coef(model_phylo_pituitary_glands)$sex[, "Q2.5", "log10_body_size"],
    slope_high = coef(model_phylo_pituitary_glands)$sex[, "Q97.5", "log10_body_size"]
  ) %>% select(slope_low, slope_high)
)

# A tibble: 2 × 7
  sex    intercept intercept_low intercept_high slope slope_low slope_high
  <chr>      <dbl>         <dbl>          <dbl> <dbl>     <dbl>      <dbl>
1 female     -3.85         -4.64          -3.11 0.775     0.524      1.04 
2 male       -3.85         -4.56          -3.15 0.745     0.517      0.973

Model comparison

loo_simple <- loo(model_simple_pituitary_glands)
loo_phylo <- loo(model_phylo_pituitary_glands)
loo_compare(loo_simple, loo_phylo)

                              elpd_diff se_diff
model_phylo_pituitary_glands    0.0       0.0  
model_simple_pituitary_glands -11.8       1.6  

In this comparison, the phylogenetic model is clearly preferable to the model treating species as independent observational units. The phylogenetic model (model_phylo_pituitary_glands) serves as reference with elpd_diff = 0, while the simpler model (model_simple_pituitary_glands) has an expected log predictive density (ELPD) that is about 11.8 units lower, with an associated standard error of 1.6. In leave-one-out cross-validation, differences exceeding approximately four times their standard error indicate strong evidence. Here, 11.8/1.6≈7.4, indicating a large difference relative to its uncertainty, and thus providing evidence favouring the phylogenetic model.

Spleen

Filter data and prune tree

# Filter organ of interest
dat_organ <- dat_mean %>%
  group_by(organ_grouped, sex) %>%
  ungroup() %>%
  filter(organ_grouped == "spleen")

# Check for database/tree species mismatches
setdiff(tree$tip.label, dat_organ$phylo)      # Tree species absent in database

  [1] "Petrochromis_famula"             "Lobochilotes_labiatus"          
  [3] "Simochromis_diagramma"           "Petrochromis_orthognathus"      
  [5] "Ctenochromis_horei"              "Limnotilapia_dardennii"         
  [7] "Gnathochromis_pfefferi"          "Tropheus_moorii"                
  [9] "Tropheus_brichardi"              "Haplotaxodon_microlepis"        
 [11] "Perissodus_microlepis"           "Benthochromis_tricoti"          
 [13] "Cyprichromis_leptosoma"          "Greenwoodochromis_christyi"     
 [15] "Limnochromis_staneri"            "Gnathochromis_permaxillaris"    
 [17] "Baileychromis_centropomoides"    "Triglachromis_otostigma"        
 [19] "Cyphotilapia_frontosa"           "Ophthalmotilapia_ventralis"     
 [21] "Cyathopharynx_furcifer"          "Ophthalmotilapia_nasuta"        
 [23] "Ophthalmotilapia_boops"          "Aulonocranus_dewindti"          
 [25] "Ectodus_descampsii"              "Callochromis_melanostigma"      
 [27] "Xenotilapia_melanogenys"         "Xenotilapia_flavipinnis"        
 [29] "Eretmodus_cyanostictus"          "Tanganicodus_irsacae"           
 [31] "Spathodus_marlieri"              "Spathodus_erythrodon"           
 [33] "Neolamprologus_brichardi"        "Julidochromis_regani"           
 [35] "Julidochromis_marlieri"          "Neolamprologus_tetracanthus"    
 [37] "Telmatochromis_temporalis"       "Lepidiolamprologus_elongatus"   
 [39] "Lepidiolamprologus_profundicola" "Lepidiolamprologus_nkambae"     
 [41] "Altolamprologus_compressiceps"   "Neolamprologus_brevis"          
 [43] "Chalinochromis_brichardi"        "Julidochromis_ornatus"          
 [45] "Lamprologus_ornatipinnis"        "Neolamprologus_pulcher"         
 [47] "Variabilichromis_moorii"         "Neolamprologus_tretocephalus"   
 [49] "Neolamprologus_sexfasciatus"     "Trematocara_unimaculatum"       
 [51] "Hemibates_stenosoma"             "Bathybates_ferox"               
 [53] "Oreochromis_niloticus"           "Poecilia_reticulata"            
 [55] "Nerophis_lumbriciformis"         "Syngnathus_abaster"             
 [57] "Syngnathus_schlegeli"            "Hippocampus_comes"              
 [59] "Hippocampus_trimaculatus"        "Hippocampus_kuda"               
 [61] "Hippocampus_spinosissimus"       "Hippocampus_abdominalis"        
 [63] "Hippichthys_cyanospilos"         "Syngnathoides_biaculeatus"      
 [65] "Corythoichthys_intestinalis"     "Corythoichthys_haematopterus"   
 [67] "Doryichthys_boaja"               "Doryrhamphus_japonicus"         
 [69] "Entelurus_aequoreus"             "Anguilla_anguilla"              
 [71] "Syncerus_caffer"                 "Tragelaphus_eurycerus"          
 [73] "Tragelaphus_scriptus"            "Redunca_arundinum"              
 [75] "Antidorcas_marsupialis"          "Nanger_granti"                  
 [77] "Raphicerus_campestris"           "Cephalophus_natalensis"         
 [79] "Sylvicapra_grimmia"              "Oreotragus_oreotragus"          
 [81] "Hippotragus_equinus"             "Addax_nasomaculatus"            
 [83] "Oryx_gazella"                    "Connochaetes_taurinus"          
 [85] "Damaliscus_lunatus"              "Ovis_aries"                     
 [87] "Aepyceros_melampus"              "Pteropus_alecto"                
 [89] "Pteropus_poliocephalus"          "Pteropus_scapulatus"            
 [91] "Chlorocebus_aethiops"            "Chlorocebus_sabaeus"            
 [93] "Erythrocebus_patas"              "Papio_cynocephalus"             
 [95] "Papio_anubis"                    "Cercocebus_atys"                
 [97] "Macaca_arctoides"                "Macaca_mulatta"                 
 [99] "Macaca_maura"                    "Trachypithecus_francoisi"       
[101] "Phrynocephalus_vlangalii"        "Platycercus_caledonicus"        
[103] "Platycercus_elegans"             "Platycercus_eximius"            
[105] "Platycercus_venustus"            "Barnardius_zonarius"            
[107] "Northiella_haematogaster"        "Lathamus_discolor"              
[109] "Neophema_chrysostoma"            "Neophema_pulchella"             
[111] "Psittacula_krameri"              "Eclectus_roratus"               
[113] "Aprosmictus_jonquillaceus"       "Alisterus_scapularis"           
[115] "Alisterus_amboinensis"           "Polytelis_alexandrae"           
[117] "Trichoglossus_haematodus"        "Agapornis_lilianae"             
[119] "Agapornis_taranta"               "Loriculus_vernalis"             
[121] "Ara_ararauna"                    "Guaruba_guarouba"               
[123] "Enicognathus_leptorhynchus"      "Pionites_melanocephalus"        
[125] "Forpus_coelestis"                "Forpus_passerinus"              
[127] "Amazona_leucocephala"            "Amazona_albifrons"              
[129] "Amazona_pretrei"                 "Amazona_vinacea"                
[131] "Amazona_finschi"                 "Amazona_amazonica"              
[133] "Amazona_ochrocephala"            "Amazona_aestiva"                
[135] "Brotogeris_versicolurus"         "Brotogeris_pyrrhoptera"         
[137] "Bolborhynchus_lineola"           "Poicephalus_gulielmi"           
[139] "Poicephalus_meyeri"              "Psittacus_erithacus"            
[141] "Serinus_mennelli"                "Padda_oryzivora"                
[143] "Lonchura_punctulata"             "Lonchura_flaviprymna"           
[145] "Lonchura_bicolor"                "Hypargos_niveoguttatus"         
[147] "Pytilia_hypogrammica"            "Uraeginthus_bengalus"           
[149] "Amandava_amandava"               "Plectrophenax_nivalis"          
[151] "Serinus_canaria"                 "Serinus_serinus"                
[153] "Loxia_curvirostra"               "Carduelis_carduelis"            
[155] "Carpodacus_roseus"               "Uragus_sibiricus"               
[157] "Pyrrhula_pyrrhula"               "Pinicola_enucleator"            
[159] "Bucanetes_githagineus"           "Coccothraustes_coccothraustes"  
[161] "Mycerobas_affinis"               "Mycerobas_carnipes"             
[163] "Fringilla_montifringilla"        "Fringilla_coelebs"              
[165] "Anas_acuta"                      "Anas_crecca"                    
[167] "Mergus_serrator"                 "Mergus_merganser"               
[169] "Bucephala_clangula"              "Melanitta_nigra"                
[171] "Somateria_mollissima"            "Callonetta_leucophrys"          
[173] "Branta_bernicla"                 "Branta_leucopsis"               
[175] "Anser_anser"                     "Cygnus_columbianus"             
[177] "Tadorna_tadorna"                 "Pucrasia_macrolopha"            
[179] "Phasianus_colchicus"             "Chrysolophus_pictus"            
[181] "Lophura_nycthemera"              "Lophura_ignita"                 
[183] "Syrmaticus_ellioti"              "Syrmaticus_reevesii"            
[185] "Perdix_perdix"                   "Tetrao_urogallus"               
[187] "Lagopus_muta"                    "Lagopus_lagopus"                
[189] "Lophophorus_impejanus"           "Tragopan_temminckii"            
[191] "Pavo_cristatus"                  "Gallus_sonneratii"              
[193] "Rollulus_rouloul"                "Arborophila_torqueola"          
[195] "Numida_meleagris"               

setdiff(dat_organ$phylo, tree$tip.label)      # Database species absent in tree

character(0)

# Exclude mismatched species for consistency
tree_organ <- keep.tip(tree, intersect(tree$tip.label, dat_organ$phylo))
dat_organ <- dat_organ[dat_organ$phylo %in% tree_organ$tip.label, ]

# Re-check for mismatches after pruning
setdiff(tree_organ$tip.label, dat_organ$phylo)

character(0)

setdiff(dat_organ$phylo, tree_organ$tip.label)

character(0)

# Is ultrametric?
is.ultrametric(tree_organ)

[1] TRUE

# how many species?
length(tree_organ$tip.label)

[1] 9

# Compute phylogenetic covariance matrix
A_organ <- ape::vcv.phylo(tree_organ, corr = TRUE)

Model simple

model_simple_spleen <- brm(
  log10_organ_size ~ 
    1 + 
    log10_body_size +
    (1 + log10_body_size | sex),
  data = dat_organ,
  family = gaussian(),
    prior = c(
    prior(normal(0.75, 0.5), class = "b", coef = "log10_body_size"),
    prior(normal(0, 2), class = "Intercept"),
    prior(student_t(3, 0, 0.5), class = "sd"),
    prior(student_t(3, 0, 1),   class = "sigma")
  ),
  seed = 6955,
  sample_prior = TRUE, 
  chains = 4, 
  cores = 4, 
  threads = threading(2),
  iter = 4000, 
  warmup = 2000,
  control = list(max_treedepth = 15,
                 adapt_delta = 0.99),
  save_pars = save_pars(all = TRUE),
  silent = TRUE, refresh = 0
)

Running MCMC with 4 parallel chains, with 2 thread(s) per chain...

Chain 1 finished in 1.4 seconds.
Chain 2 finished in 2.0 seconds.
Chain 4 finished in 2.0 seconds.
Chain 3 finished in 2.4 seconds.

All 4 chains finished successfully.
Mean chain execution time: 1.9 seconds.
Total execution time: 2.5 seconds.

Model check

pp_check(model_simple_spleen)

Model summary

summary(model_simple_spleen)

 Family: gaussian 
  Links: mu = identity 
Formula: log10_organ_size ~ 1 + log10_body_size + (1 + log10_body_size | sex) 
   Data: dat_organ (Number of observations: 18) 
  Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
         total post-warmup draws = 8000

Multilevel Hyperparameters:
~sex (Number of levels: 2) 
                               Estimate Est.Error l-95% CI u-95% CI Rhat
sd(Intercept)                      0.32      0.32     0.01     1.14 1.00
sd(log10_body_size)                0.20      0.22     0.00     0.82 1.00
cor(Intercept,log10_body_size)    -0.10      0.59    -0.97     0.94 1.00
                               Bulk_ESS Tail_ESS
sd(Intercept)                      3318     3175
sd(log10_body_size)                2302     3173
cor(Intercept,log10_body_size)     4118     4275

Regression Coefficients:
                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept          -3.01      0.39    -3.77    -2.18 1.00     4839     4271
log10_body_size     0.99      0.18     0.60     1.34 1.00     2407     1968

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.38      0.08     0.27     0.56 1.00     5903     4940

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Coeficients by sex

# Extract combined coefficients (fixed + random)
coef(model_simple_spleen)$sex %>%
{tibble(
  sex = dimnames(.)[[1]],
  intercept    = .[, "Estimate", "Intercept"],
  intercept_low  = .[, "Q2.5", "Intercept"], 
  intercept_high = .[, "Q97.5", "Intercept"]
)} %>%
bind_cols(
  tibble(
    sex = dimnames(coef(model_simple_spleen)$sex)[[1]],
    slope = coef(model_simple_spleen)$sex[, "Estimate", "log10_body_size"]
  ) %>% select(slope),
  tibble(
    sex = dimnames(coef(model_simple_spleen)$sex)[[1]],
    slope_low  = coef(model_simple_spleen)$sex[, "Q2.5", "log10_body_size"],
    slope_high = coef(model_simple_spleen)$sex[, "Q97.5", "log10_body_size"]
  ) %>% select(slope_low, slope_high)
)

# A tibble: 2 × 7
  sex    intercept intercept_low intercept_high slope slope_low slope_high
  <chr>      <dbl>         <dbl>          <dbl> <dbl>     <dbl>      <dbl>
1 female     -3.01         -3.58          -2.43  1.01     0.771       1.24
2 male       -3.06         -3.63          -2.49  1.02     0.787       1.25

Model with phylogeny

model_phylo_spleen <- brm(
  log10_organ_size ~ 
    1 + 
    log10_body_size +
    (1 + log10_body_size | sex) +
    (1 | gr(phylo, cov = A_organ)),
  data = dat_organ,
  data2 = list(A_organ = A_organ),
  family = gaussian(),
    prior = c(
    prior(normal(0.75, 0.5), class = "b", coef = "log10_body_size"),
    prior(normal(0, 2), class = "Intercept"),
    prior(student_t(3, 0, 0.5), class = "sd"),
    prior(student_t(3, 0, 1),   class = "sigma")
  ),
  seed = 6955,
  sample_prior = TRUE, 
  chains = 4, 
  cores = 4, 
  threads = threading(2),
  iter = 4000, 
  warmup = 2000,
  control = list(max_treedepth = 15,
                 adapt_delta = 0.99),
  save_pars = save_pars(all = TRUE),
  silent = TRUE, refresh = 0
)

Running MCMC with 4 parallel chains, with 2 thread(s) per chain...

Chain 1 finished in 11.2 seconds.
Chain 2 finished in 12.9 seconds.
Chain 3 finished in 14.5 seconds.
Chain 4 finished in 16.5 seconds.

All 4 chains finished successfully.
Mean chain execution time: 13.8 seconds.
Total execution time: 16.6 seconds.

Model check

pp_check(model_phylo_spleen)

Model summary

summary(model_phylo_spleen)

 Family: gaussian 
  Links: mu = identity 
Formula: log10_organ_size ~ 1 + log10_body_size + (1 + log10_body_size | sex) + (1 | gr(phylo, cov = A_organ)) 
   Data: dat_organ (Number of observations: 18) 
  Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
         total post-warmup draws = 8000

Multilevel Hyperparameters:
~phylo (Number of levels: 9) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     1.20      0.38     0.53     2.07 1.00      834      407

~sex (Number of levels: 2) 
                               Estimate Est.Error l-95% CI u-95% CI Rhat
sd(Intercept)                      0.29      0.32     0.01     1.08 1.00
sd(log10_body_size)                0.18      0.22     0.00     0.81 1.00
cor(Intercept,log10_body_size)    -0.12      0.59    -0.98     0.93 1.00
                               Bulk_ESS Tail_ESS
sd(Intercept)                      2914     3220
sd(log10_body_size)                2519     3080
cor(Intercept,log10_body_size)     6313     4358

Regression Coefficients:
                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept          -2.71      0.89    -4.34    -0.84 1.00     2778     3472
log10_body_size     0.91      0.23     0.43     1.34 1.00     2185     4010

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.09      0.06     0.04     0.27 1.00      616      407

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Phylogenetic signal

hyp <- paste(
  "sd_phylo__Intercept^2 /",
  "(sd_phylo__Intercept^2 + sd_sex__Intercept^2 + sigma^2) = 0"
)
(hyp <- hypothesis(model_phylo_spleen, hyp, class = NULL))

Hypothesis Tests for class :
                Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
1 (sd_phylo__Interc... = 0     0.89      0.17     0.36        1       0.03
  Post.Prob Star
1      0.03    *
---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.

lambda_spleen <- hyp$hypothesis$Estimate

Coeficients by sex

# Extract combined coefficients (fixed + random)
coef(model_phylo_spleen)$sex %>%
{tibble(
  sex = dimnames(.)[[1]],
  intercept    = .[, "Estimate", "Intercept"],
  intercept_low  = .[, "Q2.5", "Intercept"], 
  intercept_high = .[, "Q97.5", "Intercept"]
)} %>%
bind_cols(
  tibble(
    sex = dimnames(coef(model_phylo_spleen)$sex)[[1]],
    slope = coef(model_phylo_spleen)$sex[, "Estimate", "log10_body_size"]
  ) %>% select(slope),
  tibble(
    sex = dimnames(coef(model_phylo_spleen)$sex)[[1]],
    slope_low  = coef(model_phylo_spleen)$sex[, "Q2.5", "log10_body_size"],
    slope_high = coef(model_phylo_spleen)$sex[, "Q97.5", "log10_body_size"]
  ) %>% select(slope_low, slope_high)
)

# A tibble: 2 × 7
  sex    intercept intercept_low intercept_high slope slope_low slope_high
  <chr>      <dbl>         <dbl>          <dbl> <dbl>     <dbl>      <dbl>
1 female     -2.68         -4.23         -0.839 0.897     0.469       1.26
2 male       -2.78         -4.32         -0.973 0.929     0.525       1.28

Model comparison

loo_simple <- loo(model_simple_spleen)
loo_phylo <- loo(model_phylo_spleen)
loo_compare(loo_simple, loo_phylo)

                    elpd_diff se_diff
model_phylo_spleen    0.0       0.0  
model_simple_spleen -20.8       1.3  

In this comparison, the phylogenetic model is clearly preferable to the model treating species as independent observational units. The phylogenetic model (model_phylo_spleen) serves as reference with elpd_diff = 0, while the simpler model (model_simple_spleen) has an expected log predictive density (ELPD) 21.4 units lower, with an associated standard error of 1.3. In leave-one-out cross-validation, differences exceeding approximately four times their standard error indicate strong evidence. Here, 21.4/1.3≈16, indicating a very large difference relative to its uncertainty, and thus providing very strong evidence favouring the phylogenetic model.

Stomach

Filter data and prune tree

# Filter organ of interest
dat_organ <- dat_mean %>%
  group_by(organ_grouped, sex) %>%
  ungroup() %>%
  filter(organ_grouped == "stomach")

# Check for database/tree species mismatches
setdiff(tree$tip.label, dat_organ$phylo)      # Tree species absent in database

  [1] "Petrochromis_famula"             "Lobochilotes_labiatus"          
  [3] "Simochromis_diagramma"           "Petrochromis_orthognathus"      
  [5] "Ctenochromis_horei"              "Limnotilapia_dardennii"         
  [7] "Gnathochromis_pfefferi"          "Tropheus_moorii"                
  [9] "Tropheus_brichardi"              "Haplotaxodon_microlepis"        
 [11] "Perissodus_microlepis"           "Benthochromis_tricoti"          
 [13] "Cyprichromis_leptosoma"          "Greenwoodochromis_christyi"     
 [15] "Limnochromis_staneri"            "Gnathochromis_permaxillaris"    
 [17] "Baileychromis_centropomoides"    "Triglachromis_otostigma"        
 [19] "Cyphotilapia_frontosa"           "Ophthalmotilapia_ventralis"     
 [21] "Cyathopharynx_furcifer"          "Ophthalmotilapia_nasuta"        
 [23] "Ophthalmotilapia_boops"          "Aulonocranus_dewindti"          
 [25] "Ectodus_descampsii"              "Callochromis_melanostigma"      
 [27] "Xenotilapia_melanogenys"         "Xenotilapia_flavipinnis"        
 [29] "Eretmodus_cyanostictus"          "Tanganicodus_irsacae"           
 [31] "Spathodus_marlieri"              "Spathodus_erythrodon"           
 [33] "Neolamprologus_brichardi"        "Julidochromis_regani"           
 [35] "Julidochromis_marlieri"          "Neolamprologus_tetracanthus"    
 [37] "Telmatochromis_temporalis"       "Lepidiolamprologus_elongatus"   
 [39] "Lepidiolamprologus_profundicola" "Lepidiolamprologus_nkambae"     
 [41] "Altolamprologus_compressiceps"   "Neolamprologus_brevis"          
 [43] "Chalinochromis_brichardi"        "Julidochromis_ornatus"          
 [45] "Lamprologus_ornatipinnis"        "Neolamprologus_pulcher"         
 [47] "Variabilichromis_moorii"         "Neolamprologus_tretocephalus"   
 [49] "Neolamprologus_sexfasciatus"     "Trematocara_unimaculatum"       
 [51] "Hemibates_stenosoma"             "Bathybates_ferox"               
 [53] "Oreochromis_niloticus"           "Poecilia_reticulata"            
 [55] "Nerophis_lumbriciformis"         "Syngnathus_abaster"             
 [57] "Syngnathus_schlegeli"            "Hippocampus_comes"              
 [59] "Hippocampus_trimaculatus"        "Hippocampus_kuda"               
 [61] "Hippocampus_spinosissimus"       "Hippocampus_abdominalis"        
 [63] "Hippichthys_cyanospilos"         "Syngnathoides_biaculeatus"      
 [65] "Corythoichthys_intestinalis"     "Corythoichthys_haematopterus"   
 [67] "Doryichthys_boaja"               "Doryrhamphus_japonicus"         
 [69] "Entelurus_aequoreus"             "Anguilla_anguilla"              
 [71] "Syncerus_caffer"                 "Tragelaphus_eurycerus"          
 [73] "Tragelaphus_scriptus"            "Redunca_arundinum"              
 [75] "Antidorcas_marsupialis"          "Nanger_granti"                  
 [77] "Raphicerus_campestris"           "Cephalophus_natalensis"         
 [79] "Sylvicapra_grimmia"              "Oreotragus_oreotragus"          
 [81] "Hippotragus_equinus"             "Addax_nasomaculatus"            
 [83] "Oryx_gazella"                    "Connochaetes_taurinus"          
 [85] "Damaliscus_lunatus"              "Ovis_aries"                     
 [87] "Aepyceros_melampus"              "Pteropus_alecto"                
 [89] "Pteropus_poliocephalus"          "Pteropus_scapulatus"            
 [91] "Rattus_norvegicus"               "Chlorocebus_aethiops"           
 [93] "Chlorocebus_sabaeus"             "Erythrocebus_patas"             
 [95] "Papio_cynocephalus"              "Papio_anubis"                   
 [97] "Cercocebus_atys"                 "Macaca_arctoides"               
 [99] "Macaca_mulatta"                  "Macaca_fascicularis"            
[101] "Macaca_maura"                    "Trachypithecus_francoisi"       
[103] "Platycercus_caledonicus"         "Platycercus_elegans"            
[105] "Platycercus_eximius"             "Platycercus_venustus"           
[107] "Barnardius_zonarius"             "Northiella_haematogaster"       
[109] "Lathamus_discolor"               "Neophema_chrysostoma"           
[111] "Neophema_pulchella"              "Psittacula_krameri"             
[113] "Eclectus_roratus"                "Aprosmictus_jonquillaceus"      
[115] "Alisterus_scapularis"            "Alisterus_amboinensis"          
[117] "Polytelis_alexandrae"            "Trichoglossus_haematodus"       
[119] "Agapornis_lilianae"              "Agapornis_taranta"              
[121] "Loriculus_vernalis"              "Ara_ararauna"                   
[123] "Guaruba_guarouba"                "Enicognathus_leptorhynchus"     
[125] "Pionites_melanocephalus"         "Forpus_coelestis"               
[127] "Forpus_passerinus"               "Amazona_leucocephala"           
[129] "Amazona_albifrons"               "Amazona_pretrei"                
[131] "Amazona_vinacea"                 "Amazona_finschi"                
[133] "Amazona_amazonica"               "Amazona_ochrocephala"           
[135] "Amazona_aestiva"                 "Brotogeris_versicolurus"        
[137] "Brotogeris_pyrrhoptera"          "Bolborhynchus_lineola"          
[139] "Poicephalus_gulielmi"            "Poicephalus_meyeri"             
[141] "Psittacus_erithacus"             "Serinus_mennelli"               
[143] "Ficedula_hypoleuca"              "Padda_oryzivora"                
[145] "Lonchura_punctulata"             "Lonchura_flaviprymna"           
[147] "Lonchura_bicolor"                "Hypargos_niveoguttatus"         
[149] "Pytilia_hypogrammica"            "Uraeginthus_bengalus"           
[151] "Amandava_amandava"               "Plectrophenax_nivalis"          
[153] "Serinus_canaria"                 "Serinus_serinus"                
[155] "Loxia_curvirostra"               "Carduelis_carduelis"            
[157] "Carpodacus_roseus"               "Uragus_sibiricus"               
[159] "Pyrrhula_pyrrhula"               "Pinicola_enucleator"            
[161] "Bucanetes_githagineus"           "Coccothraustes_coccothraustes"  
[163] "Mycerobas_affinis"               "Mycerobas_carnipes"             
[165] "Fringilla_montifringilla"        "Fringilla_coelebs"              
[167] "Anas_acuta"                      "Anas_crecca"                    
[169] "Mergus_serrator"                 "Mergus_merganser"               
[171] "Bucephala_clangula"              "Melanitta_nigra"                
[173] "Somateria_mollissima"            "Callonetta_leucophrys"          
[175] "Branta_bernicla"                 "Branta_leucopsis"               
[177] "Anser_anser"                     "Cygnus_columbianus"             
[179] "Tadorna_tadorna"                 "Pucrasia_macrolopha"            
[181] "Phasianus_colchicus"             "Chrysolophus_pictus"            
[183] "Lophura_nycthemera"              "Lophura_ignita"                 
[185] "Syrmaticus_ellioti"              "Syrmaticus_reevesii"            
[187] "Perdix_perdix"                   "Tetrao_urogallus"               
[189] "Lagopus_muta"                    "Lagopus_lagopus"                
[191] "Lophophorus_impejanus"           "Tragopan_temminckii"            
[193] "Pavo_cristatus"                  "Gallus_sonneratii"              
[195] "Rollulus_rouloul"                "Arborophila_torqueola"          

setdiff(dat_organ$phylo, tree$tip.label)      # Database species absent in tree

character(0)

# Exclude mismatched species for consistency
tree_organ <- keep.tip(tree, intersect(tree$tip.label, dat_organ$phylo))
dat_organ <- dat_organ[dat_organ$phylo %in% tree_organ$tip.label, ]

# Re-check for mismatches after pruning
setdiff(tree_organ$tip.label, dat_organ$phylo)

character(0)

setdiff(dat_organ$phylo, tree_organ$tip.label)

character(0)

# Is ultrametric?
is.ultrametric(tree_organ)

[1] TRUE

# how many species?
length(tree_organ$tip.label)

[1] 8

# Compute phylogenetic covariance matrix
A_organ <- ape::vcv.phylo(tree_organ, corr = TRUE)

Model simple

model_simple_stomach <- brm(
  log10_organ_size ~ 
    1 + 
    log10_body_size +
    (1 + log10_body_size | sex),
  data = dat_organ,
  family = gaussian(),
    prior = c(
    prior(normal(0.75, 0.5), class = "b", coef = "log10_body_size"),
    prior(normal(0, 2), class = "Intercept"),
    prior(student_t(3, 0, 0.5), class = "sd"),
    prior(student_t(3, 0, 1),   class = "sigma")
  ),
  seed = 6955,
  sample_prior = TRUE, 
  chains = 4, 
  cores = 4, 
  threads = threading(2),
  iter = 4000, 
  warmup = 2000,
  control = list(max_treedepth = 15,
                 adapt_delta = 0.99),
  save_pars = save_pars(all = TRUE),
  silent = TRUE, refresh = 0
)

Running MCMC with 4 parallel chains, with 2 thread(s) per chain...

Chain 1 finished in 2.4 seconds.
Chain 2 finished in 2.6 seconds.
Chain 3 finished in 2.5 seconds.
Chain 4 finished in 2.6 seconds.

All 4 chains finished successfully.
Mean chain execution time: 2.5 seconds.
Total execution time: 2.7 seconds.

Model check

pp_check(model_simple_stomach)

Model summary

summary(model_simple_stomach)

 Family: gaussian 
  Links: mu = identity 
Formula: log10_organ_size ~ 1 + log10_body_size + (1 + log10_body_size | sex) 
   Data: dat_organ (Number of observations: 16) 
  Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
         total post-warmup draws = 8000

Multilevel Hyperparameters:
~sex (Number of levels: 2) 
                               Estimate Est.Error l-95% CI u-95% CI Rhat
sd(Intercept)                      0.29      0.30     0.01     1.08 1.00
sd(log10_body_size)                0.20      0.23     0.00     0.84 1.00
cor(Intercept,log10_body_size)    -0.09      0.60    -0.98     0.95 1.00
                               Bulk_ESS Tail_ESS
sd(Intercept)                      3636     3639
sd(log10_body_size)                1837     3384
cor(Intercept,log10_body_size)     3837     4182

Regression Coefficients:
                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept          -2.63      0.33    -3.28    -1.95 1.00     4295     3726
log10_body_size     1.18      0.18     0.71     1.48 1.00     2080     1216

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.23      0.05     0.15     0.35 1.00     5037     3724

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Coeficients by sex

# Extract combined coefficients (fixed + random)
coef(model_simple_stomach)$sex %>%
{tibble(
  sex = dimnames(.)[[1]],
  intercept    = .[, "Estimate", "Intercept"],
  intercept_low  = .[, "Q2.5", "Intercept"], 
  intercept_high = .[, "Q97.5", "Intercept"]
)} %>%
bind_cols(
  tibble(
    sex = dimnames(coef(model_simple_stomach)$sex)[[1]],
    slope = coef(model_simple_stomach)$sex[, "Estimate", "log10_body_size"]
  ) %>% select(slope),
  tibble(
    sex = dimnames(coef(model_simple_stomach)$sex)[[1]],
    slope_low  = coef(model_simple_stomach)$sex[, "Q2.5", "log10_body_size"],
    slope_high = coef(model_simple_stomach)$sex[, "Q97.5", "log10_body_size"]
  ) %>% select(slope_low, slope_high)
)

# A tibble: 2 × 7
  sex    intercept intercept_low intercept_high slope slope_low slope_high
  <chr>      <dbl>         <dbl>          <dbl> <dbl>     <dbl>      <dbl>
1 female     -2.63         -3.03          -2.22  1.23      1.05       1.39
2 male       -2.64         -3.04          -2.23  1.23      1.06       1.39

Model with phylogeny

model_phylo_stomach <- brm(
  log10_organ_size ~ 
    1 + 
    log10_body_size +
    (1 + log10_body_size | sex) +
    (1 | gr(phylo, cov = A_organ)),
  data = dat_organ,
  data2 = list(A_organ = A_organ),
  family = gaussian(),
    prior = c(
    prior(normal(0.75, 0.5), class = "b", coef = "log10_body_size"),
    prior(normal(0, 2), class = "Intercept"),
    prior(student_t(3, 0, 0.5), class = "sd"),
    prior(student_t(3, 0, 1),   class = "sigma")
  ),
  seed = 6955,
  sample_prior = TRUE, 
  chains = 4, 
  cores = 4, 
  threads = threading(2),
  iter = 4000, 
  warmup = 2000,
  control = list(max_treedepth = 15,
                 adapt_delta = 0.99),
  save_pars = save_pars(all = TRUE),
  silent = TRUE, refresh = 0
)

Running MCMC with 4 parallel chains, with 2 thread(s) per chain...

Chain 2 finished in 14.3 seconds.
Chain 1 finished in 16.7 seconds.
Chain 4 finished in 16.8 seconds.
Chain 3 finished in 22.1 seconds.

All 4 chains finished successfully.
Mean chain execution time: 17.5 seconds.
Total execution time: 22.2 seconds.

Model check

pp_check(model_phylo_stomach)

Model summary

summary(model_phylo_stomach)

 Family: gaussian 
  Links: mu = identity 
Formula: log10_organ_size ~ 1 + log10_body_size + (1 + log10_body_size | sex) + (1 | gr(phylo, cov = A_organ)) 
   Data: dat_organ (Number of observations: 16) 
  Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
         total post-warmup draws = 8000

Multilevel Hyperparameters:
~phylo (Number of levels: 8) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.39      0.13     0.22     0.70 1.00     2453     4006

~sex (Number of levels: 2) 
                               Estimate Est.Error l-95% CI u-95% CI Rhat
sd(Intercept)                      0.21      0.27     0.00     0.93 1.00
sd(log10_body_size)                0.13      0.18     0.00     0.63 1.00
cor(Intercept,log10_body_size)    -0.05      0.62    -0.98     0.96 1.00
                               Bulk_ESS Tail_ESS
sd(Intercept)                      3348     3977
sd(log10_body_size)                2424     3143
cor(Intercept,log10_body_size)     5595     5115

Regression Coefficients:
                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept          -1.96      0.42    -2.74    -1.09 1.00     3111     3808
log10_body_size     0.87      0.17     0.52     1.20 1.00     2549     3676

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.04      0.02     0.02     0.08 1.00     2226     3234

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Phylogenetic signal

hyp <- paste(
  "sd_phylo__Intercept^2 /",
  "(sd_phylo__Intercept^2 + sd_sex__Intercept^2 + sigma^2) = 0"
)
(hyp <- hypothesis(model_phylo_stomach, hyp, class = NULL))

Hypothesis Tests for class :
                Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
1 (sd_phylo__Interc... = 0     0.77      0.26     0.12     0.99       0.02
  Post.Prob Star
1      0.02    *
---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.

lambda_stomach <- hyp$hypothesis$Estimate

Coeficients by sex

# Extract combined coefficients (fixed + random)
coef(model_phylo_stomach)$sex %>%
{tibble(
  sex = dimnames(.)[[1]],
  intercept    = .[, "Estimate", "Intercept"],
  intercept_low  = .[, "Q2.5", "Intercept"], 
  intercept_high = .[, "Q97.5", "Intercept"]
)} %>%
bind_cols(
  tibble(
    sex = dimnames(coef(model_phylo_stomach)$sex)[[1]],
    slope = coef(model_phylo_stomach)$sex[, "Estimate", "log10_body_size"]
  ) %>% select(slope),
  tibble(
    sex = dimnames(coef(model_phylo_stomach)$sex)[[1]],
    slope_low  = coef(model_phylo_stomach)$sex[, "Q2.5", "log10_body_size"],
    slope_high = coef(model_phylo_stomach)$sex[, "Q97.5", "log10_body_size"]
  ) %>% select(slope_low, slope_high)
)

# A tibble: 2 × 7
  sex    intercept intercept_low intercept_high slope slope_low slope_high
  <chr>      <dbl>         <dbl>          <dbl> <dbl>     <dbl>      <dbl>
1 female     -1.94         -2.63          -1.17 0.872     0.586       1.12
2 male       -1.97         -2.65          -1.21 0.878     0.594       1.13

Model comparison

loo_simple <- loo(model_simple_stomach)
loo_phylo <- loo(model_phylo_stomach)
loo_compare(loo_simple, loo_phylo)

                     elpd_diff se_diff
model_phylo_stomach    0.0       0.0  
model_simple_stomach -23.3       3.6  

In this comparison, the phylogenetic model is clearly preferable to the model treating species as independent observational units. The phylogenetic model (model_phylo_stomach) is used as the reference with elpd_diff = 0, while the simpler model (model_simple_stomach) has an expected log predictive density (ELPD) that is about 23.3 units lower, with an associated standard error of 3.6. In leave-one-out cross-validation, differences exceeding approximately four times their standard error indicate strong evidence. Here, 23.3/3.6≈6.5, providing strong evidence favouring the phylogenetic model..

Lungs

Filter data and prune tree

# Filter organ of interest
dat_organ <- dat_mean %>%
  group_by(organ_grouped, sex) %>%
  ungroup() %>%
  filter(organ_grouped == "lungs")

# Check for database/tree species mismatches
setdiff(tree$tip.label, dat_organ$phylo)      # Tree species absent in database

  [1] "Petrochromis_famula"             "Lobochilotes_labiatus"          
  [3] "Simochromis_diagramma"           "Petrochromis_orthognathus"      
  [5] "Ctenochromis_horei"              "Limnotilapia_dardennii"         
  [7] "Gnathochromis_pfefferi"          "Tropheus_moorii"                
  [9] "Tropheus_brichardi"              "Haplotaxodon_microlepis"        
 [11] "Perissodus_microlepis"           "Benthochromis_tricoti"          
 [13] "Cyprichromis_leptosoma"          "Greenwoodochromis_christyi"     
 [15] "Limnochromis_staneri"            "Gnathochromis_permaxillaris"    
 [17] "Baileychromis_centropomoides"    "Triglachromis_otostigma"        
 [19] "Cyphotilapia_frontosa"           "Ophthalmotilapia_ventralis"     
 [21] "Cyathopharynx_furcifer"          "Ophthalmotilapia_nasuta"        
 [23] "Ophthalmotilapia_boops"          "Aulonocranus_dewindti"          
 [25] "Ectodus_descampsii"              "Callochromis_melanostigma"      
 [27] "Xenotilapia_melanogenys"         "Xenotilapia_flavipinnis"        
 [29] "Eretmodus_cyanostictus"          "Tanganicodus_irsacae"           
 [31] "Spathodus_marlieri"              "Spathodus_erythrodon"           
 [33] "Neolamprologus_brichardi"        "Julidochromis_regani"           
 [35] "Julidochromis_marlieri"          "Neolamprologus_tetracanthus"    
 [37] "Telmatochromis_temporalis"       "Lepidiolamprologus_elongatus"   
 [39] "Lepidiolamprologus_profundicola" "Lepidiolamprologus_nkambae"     
 [41] "Altolamprologus_compressiceps"   "Neolamprologus_brevis"          
 [43] "Chalinochromis_brichardi"        "Julidochromis_ornatus"          
 [45] "Lamprologus_ornatipinnis"        "Neolamprologus_pulcher"         
 [47] "Variabilichromis_moorii"         "Neolamprologus_tretocephalus"   
 [49] "Neolamprologus_sexfasciatus"     "Trematocara_unimaculatum"       
 [51] "Hemibates_stenosoma"             "Bathybates_ferox"               
 [53] "Oreochromis_niloticus"           "Poecilia_reticulata"            
 [55] "Nerophis_lumbriciformis"         "Syngnathus_abaster"             
 [57] "Syngnathus_schlegeli"            "Hippocampus_comes"              
 [59] "Hippocampus_trimaculatus"        "Hippocampus_kuda"               
 [61] "Hippocampus_spinosissimus"       "Hippocampus_abdominalis"        
 [63] "Hippichthys_cyanospilos"         "Syngnathoides_biaculeatus"      
 [65] "Corythoichthys_intestinalis"     "Corythoichthys_haematopterus"   
 [67] "Doryichthys_boaja"               "Doryrhamphus_japonicus"         
 [69] "Entelurus_aequoreus"             "Anguilla_anguilla"              
 [71] "Syncerus_caffer"                 "Tragelaphus_eurycerus"          
 [73] "Tragelaphus_scriptus"            "Redunca_arundinum"              
 [75] "Antidorcas_marsupialis"          "Nanger_granti"                  
 [77] "Raphicerus_campestris"           "Cephalophus_natalensis"         
 [79] "Sylvicapra_grimmia"              "Oreotragus_oreotragus"          
 [81] "Hippotragus_equinus"             "Addax_nasomaculatus"            
 [83] "Oryx_gazella"                    "Connochaetes_taurinus"          
 [85] "Damaliscus_lunatus"              "Ovis_aries"                     
 [87] "Aepyceros_melampus"              "Pteropus_alecto"                
 [89] "Pteropus_poliocephalus"          "Pteropus_scapulatus"            
 [91] "Chlorocebus_aethiops"            "Chlorocebus_sabaeus"            
 [93] "Erythrocebus_patas"              "Papio_cynocephalus"             
 [95] "Papio_anubis"                    "Cercocebus_atys"                
 [97] "Macaca_arctoides"                "Macaca_mulatta"                 
 [99] "Macaca_maura"                    "Trachypithecus_francoisi"       
[101] "Platycercus_caledonicus"         "Platycercus_elegans"            
[103] "Platycercus_eximius"             "Platycercus_venustus"           
[105] "Barnardius_zonarius"             "Northiella_haematogaster"       
[107] "Lathamus_discolor"               "Neophema_chrysostoma"           
[109] "Neophema_pulchella"              "Psittacula_krameri"             
[111] "Eclectus_roratus"                "Aprosmictus_jonquillaceus"      
[113] "Alisterus_scapularis"            "Alisterus_amboinensis"          
[115] "Polytelis_alexandrae"            "Trichoglossus_haematodus"       
[117] "Agapornis_lilianae"              "Agapornis_taranta"              
[119] "Loriculus_vernalis"              "Ara_ararauna"                   
[121] "Guaruba_guarouba"                "Enicognathus_leptorhynchus"     
[123] "Pionites_melanocephalus"         "Forpus_coelestis"               
[125] "Forpus_passerinus"               "Amazona_leucocephala"           
[127] "Amazona_albifrons"               "Amazona_pretrei"                
[129] "Amazona_vinacea"                 "Amazona_finschi"                
[131] "Amazona_amazonica"               "Amazona_ochrocephala"           
[133] "Amazona_aestiva"                 "Brotogeris_versicolurus"        
[135] "Brotogeris_pyrrhoptera"          "Bolborhynchus_lineola"          
[137] "Poicephalus_gulielmi"            "Poicephalus_meyeri"             
[139] "Psittacus_erithacus"             "Serinus_mennelli"               
[141] "Ficedula_hypoleuca"              "Padda_oryzivora"                
[143] "Lonchura_punctulata"             "Lonchura_flaviprymna"           
[145] "Lonchura_bicolor"                "Hypargos_niveoguttatus"         
[147] "Pytilia_hypogrammica"            "Uraeginthus_bengalus"           
[149] "Amandava_amandava"               "Plectrophenax_nivalis"          
[151] "Serinus_canaria"                 "Serinus_serinus"                
[153] "Loxia_curvirostra"               "Carduelis_carduelis"            
[155] "Carpodacus_roseus"               "Uragus_sibiricus"               
[157] "Pyrrhula_pyrrhula"               "Pinicola_enucleator"            
[159] "Bucanetes_githagineus"           "Coccothraustes_coccothraustes"  
[161] "Mycerobas_affinis"               "Mycerobas_carnipes"             
[163] "Fringilla_montifringilla"        "Fringilla_coelebs"              
[165] "Anas_platyrhynchos"              "Anas_acuta"                     
[167] "Anas_crecca"                     "Mergus_serrator"                
[169] "Mergus_merganser"                "Bucephala_clangula"             
[171] "Melanitta_nigra"                 "Somateria_mollissima"           
[173] "Callonetta_leucophrys"           "Branta_bernicla"                
[175] "Branta_leucopsis"                "Anser_anser"                    
[177] "Cygnus_columbianus"              "Tadorna_tadorna"                
[179] "Coturnix_japonica"               "Pucrasia_macrolopha"            
[181] "Phasianus_colchicus"             "Chrysolophus_pictus"            
[183] "Lophura_nycthemera"              "Lophura_ignita"                 
[185] "Syrmaticus_ellioti"              "Syrmaticus_reevesii"            
[187] "Perdix_perdix"                   "Tetrao_urogallus"               
[189] "Lagopus_muta"                    "Lagopus_lagopus"                
[191] "Lophophorus_impejanus"           "Tragopan_temminckii"            
[193] "Pavo_cristatus"                  "Gallus_gallus"                  
[195] "Gallus_sonneratii"               "Rollulus_rouloul"               
[197] "Arborophila_torqueola"           "Numida_meleagris"               

setdiff(dat_organ$phylo, tree$tip.label)      # Database species absent in tree

character(0)

# Exclude mismatched species for consistency
tree_organ <- keep.tip(tree, intersect(tree$tip.label, dat_organ$phylo))
dat_organ <- dat_organ[dat_organ$phylo %in% tree_organ$tip.label, ]

# Re-check for mismatches after pruning
setdiff(tree_organ$tip.label, dat_organ$phylo)

character(0)

setdiff(dat_organ$phylo, tree_organ$tip.label)

character(0)

# Is ultrametric?
is.ultrametric(tree_organ)

[1] TRUE

# how many species?
length(tree_organ$tip.label)

[1] 6

# Compute phylogenetic covariance matrix
A_organ <- ape::vcv.phylo(tree_organ, corr = TRUE)

Model simple

model_simple_lungs <- brm(
  log10_organ_size ~ 
    1 + 
    log10_body_size +
    (1 + log10_body_size | sex),
  data = dat_organ,
  family = gaussian(),
    prior = c(
    prior(normal(0.75, 0.5), class = "b", coef = "log10_body_size"),
    prior(normal(0, 2), class = "Intercept"),
    prior(student_t(3, 0, 0.5), class = "sd"),
    prior(student_t(3, 0, 1),   class = "sigma")
  ),
  seed = 6955,
  sample_prior = TRUE, 
  chains = 4, 
  cores = 4, 
  threads = threading(2),
  iter = 4000, 
  warmup = 2000,
  control = list(max_treedepth = 15,
                 adapt_delta = 0.99),
  save_pars = save_pars(all = TRUE),
  silent = TRUE, refresh = 0
)

Running MCMC with 4 parallel chains, with 2 thread(s) per chain...

Chain 1 finished in 1.9 seconds.
Chain 4 finished in 1.7 seconds.
Chain 3 finished in 1.9 seconds.
Chain 2 finished in 2.2 seconds.

All 4 chains finished successfully.
Mean chain execution time: 1.9 seconds.
Total execution time: 2.3 seconds.

Model check

pp_check(model_simple_lungs)

Model summary

summary(model_simple_lungs)

 Family: gaussian 
  Links: mu = identity 
Formula: log10_organ_size ~ 1 + log10_body_size + (1 + log10_body_size | sex) 
   Data: dat_organ (Number of observations: 12) 
  Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
         total post-warmup draws = 8000

Multilevel Hyperparameters:
~sex (Number of levels: 2) 
                               Estimate Est.Error l-95% CI u-95% CI Rhat
sd(Intercept)                      0.31      0.33     0.01     1.16 1.00
sd(log10_body_size)                0.21      0.23     0.00     0.81 1.00
cor(Intercept,log10_body_size)    -0.08      0.60    -0.97     0.94 1.00
                               Bulk_ESS Tail_ESS
sd(Intercept)                      3402     3559
sd(log10_body_size)                2134     2245
cor(Intercept,log10_body_size)     3557     4522

Regression Coefficients:
                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept          -2.52      0.36    -3.22    -1.75 1.00     4350     4283
log10_body_size     1.05      0.19     0.60     1.38 1.00     2581     2303

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.27      0.08     0.17     0.46 1.00     4371     4081

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Coeficients by sex

# Extract combined coefficients (fixed + random)
coef(model_simple_lungs)$sex %>%
{tibble(
  sex = dimnames(.)[[1]],
  intercept    = .[, "Estimate", "Intercept"],
  intercept_low  = .[, "Q2.5", "Intercept"], 
  intercept_high = .[, "Q97.5", "Intercept"]
)} %>%
bind_cols(
  tibble(
    sex = dimnames(coef(model_simple_lungs)$sex)[[1]],
    slope = coef(model_simple_lungs)$sex[, "Estimate", "log10_body_size"]
  ) %>% select(slope),
  tibble(
    sex = dimnames(coef(model_simple_lungs)$sex)[[1]],
    slope_low  = coef(model_simple_lungs)$sex[, "Q2.5", "log10_body_size"],
    slope_high = coef(model_simple_lungs)$sex[, "Q97.5", "log10_body_size"]
  ) %>% select(slope_low, slope_high)
)

# A tibble: 2 × 7
  sex    intercept intercept_low intercept_high slope slope_low slope_high
  <chr>      <dbl>         <dbl>          <dbl> <dbl>     <dbl>      <dbl>
1 female     -2.54         -3.04          -2.03  1.10     0.858       1.33
2 male       -2.52         -3.01          -2.00  1.07     0.838       1.28

Model with phylogeny

model_phylo_lungs <- brm(
  log10_organ_size ~ 
    1 + 
    log10_body_size +
    (1 + log10_body_size | sex) +
    (1 | gr(phylo, cov = A_organ)),
  data = dat_organ,
  data2 = list(A_organ = A_organ),
  family = gaussian(),
    prior = c(
    prior(normal(0.75, 0.5), class = "b", coef = "log10_body_size"),
    prior(normal(0, 2), class = "Intercept"),
    prior(student_t(3, 0, 0.5), class = "sd"),
    prior(student_t(3, 0, 1),   class = "sigma")
  ),
  seed = 6955,
  sample_prior = TRUE, 
  chains = 4, 
  cores = 4, 
  threads = threading(2),
  iter = 4000, 
  warmup = 2000,
  control = list(max_treedepth = 15,
                 adapt_delta = 0.99),
  save_pars = save_pars(all = TRUE),
  silent = TRUE, refresh = 0
)

Running MCMC with 4 parallel chains, with 2 thread(s) per chain...

Chain 3 finished in 8.6 seconds.
Chain 4 finished in 9.7 seconds.
Chain 2 finished in 10.9 seconds.
Chain 1 finished in 12.1 seconds.

All 4 chains finished successfully.
Mean chain execution time: 10.3 seconds.
Total execution time: 12.1 seconds.

Model check

pp_check(model_phylo_lungs)

Model summary

summary(model_phylo_lungs)

 Family: gaussian 
  Links: mu = identity 
Formula: log10_organ_size ~ 1 + log10_body_size + (1 + log10_body_size | sex) + (1 | gr(phylo, cov = A_organ)) 
   Data: dat_organ (Number of observations: 12) 
  Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
         total post-warmup draws = 8000

Multilevel Hyperparameters:
~phylo (Number of levels: 6) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.48      0.24     0.18     1.05 1.00     1882     3320

~sex (Number of levels: 2) 
                               Estimate Est.Error l-95% CI u-95% CI Rhat
sd(Intercept)                      0.29      0.29     0.01     1.08 1.00
sd(log10_body_size)                0.21      0.22     0.01     0.82 1.00
cor(Intercept,log10_body_size)    -0.12      0.59    -0.98     0.94 1.00
                               Bulk_ESS Tail_ESS
sd(Intercept)                      3492     3380
sd(log10_body_size)                2708     3539
cor(Intercept,log10_body_size)     5438     4491

Regression Coefficients:
                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept          -2.37      0.52    -3.35    -1.23 1.00     4180     4323
log10_body_size     0.93      0.20     0.50     1.33 1.00     2609     3012

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.08      0.04     0.03     0.18 1.00     1269     2514

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Phylogenetic signal

hyp <- paste(
  "sd_phylo__Intercept^2 /",
  "(sd_phylo__Intercept^2 + sd_sex__Intercept^2 + sigma^2) = 0"
)
(hyp <- hypothesis(model_phylo_lungs, hyp, class = NULL))

Hypothesis Tests for class :
                Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
1 (sd_phylo__Interc... = 0     0.69      0.27     0.08     0.99       0.04
  Post.Prob Star
1      0.04    *
---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.

lambda_lungs <- hyp$hypothesis$Estimate

Coeficients by sex

# Extract combined coefficients (fixed + random)
coef(model_phylo_lungs)$sex %>%
{tibble(
  sex = dimnames(.)[[1]],
  intercept    = .[, "Estimate", "Intercept"],
  intercept_low  = .[, "Q2.5", "Intercept"], 
  intercept_high = .[, "Q97.5", "Intercept"]
)} %>%
bind_cols(
  tibble(
    sex = dimnames(coef(model_phylo_lungs)$sex)[[1]],
    slope = coef(model_phylo_lungs)$sex[, "Estimate", "log10_body_size"]
  ) %>% select(slope),
  tibble(
    sex = dimnames(coef(model_phylo_lungs)$sex)[[1]],
    slope_low  = coef(model_phylo_lungs)$sex[, "Q2.5", "log10_body_size"],
    slope_high = coef(model_phylo_lungs)$sex[, "Q97.5", "log10_body_size"]
  ) %>% select(slope_low, slope_high)
)

# A tibble: 2 × 7
  sex    intercept intercept_low intercept_high slope slope_low slope_high
  <chr>      <dbl>         <dbl>          <dbl> <dbl>     <dbl>      <dbl>
1 female     -2.43         -3.23          -1.39 0.971     0.648       1.22
2 male       -2.34         -3.11          -1.32 0.914     0.620       1.14

Model comparison

loo_simple <- loo(model_simple_lungs)
loo_phylo <- loo(model_phylo_lungs)
loo_compare(loo_simple, loo_phylo)

                   elpd_diff se_diff
model_phylo_lungs    0.0       0.0  
model_simple_lungs -13.4       1.3  

In this comparison, the model that accounts for phylogeny is clearly preferable to the model that treats species as independent observational units. The phylogenetic model (model_phylo_lungs) is used as the reference with elpd_diff = 0, while the simpler model (model_simple_lungs) has an expected log predictive density (ELPD) that is about 12.7 units lower, with an associated standard error of 1.4. In the context of leave‑one‑out cross‑validation, differences exceeding roughly four times their standard error are typically regarded as strong evidence. Here, 12.7/1.4≈9.1, indicating a very large difference relative to its uncertainty, and thus providing strong support for the phylogenetic model.

Adrenal glands

Filter data and prune tree

# Filter organ of interest
dat_organ <- dat_mean %>%
  group_by(organ_grouped, sex) %>%
  ungroup() %>%
  filter(organ_grouped == "adrenal glands")

# Check for database/tree species mismatches
setdiff(tree$tip.label, dat_organ$phylo)      # Tree species absent in database

  [1] "Petrochromis_famula"             "Lobochilotes_labiatus"          
  [3] "Simochromis_diagramma"           "Petrochromis_orthognathus"      
  [5] "Ctenochromis_horei"              "Limnotilapia_dardennii"         
  [7] "Gnathochromis_pfefferi"          "Tropheus_moorii"                
  [9] "Tropheus_brichardi"              "Haplotaxodon_microlepis"        
 [11] "Perissodus_microlepis"           "Benthochromis_tricoti"          
 [13] "Cyprichromis_leptosoma"          "Greenwoodochromis_christyi"     
 [15] "Limnochromis_staneri"            "Gnathochromis_permaxillaris"    
 [17] "Baileychromis_centropomoides"    "Triglachromis_otostigma"        
 [19] "Cyphotilapia_frontosa"           "Ophthalmotilapia_ventralis"     
 [21] "Cyathopharynx_furcifer"          "Ophthalmotilapia_nasuta"        
 [23] "Ophthalmotilapia_boops"          "Aulonocranus_dewindti"          
 [25] "Ectodus_descampsii"              "Callochromis_melanostigma"      
 [27] "Xenotilapia_melanogenys"         "Xenotilapia_flavipinnis"        
 [29] "Eretmodus_cyanostictus"          "Tanganicodus_irsacae"           
 [31] "Spathodus_marlieri"              "Spathodus_erythrodon"           
 [33] "Neolamprologus_brichardi"        "Julidochromis_regani"           
 [35] "Julidochromis_marlieri"          "Neolamprologus_tetracanthus"    
 [37] "Telmatochromis_temporalis"       "Lepidiolamprologus_elongatus"   
 [39] "Lepidiolamprologus_profundicola" "Lepidiolamprologus_nkambae"     
 [41] "Altolamprologus_compressiceps"   "Neolamprologus_brevis"          
 [43] "Chalinochromis_brichardi"        "Julidochromis_ornatus"          
 [45] "Lamprologus_ornatipinnis"        "Neolamprologus_pulcher"         
 [47] "Variabilichromis_moorii"         "Neolamprologus_tretocephalus"   
 [49] "Neolamprologus_sexfasciatus"     "Trematocara_unimaculatum"       
 [51] "Hemibates_stenosoma"             "Bathybates_ferox"               
 [53] "Oreochromis_niloticus"           "Poecilia_reticulata"            
 [55] "Nerophis_lumbriciformis"         "Syngnathus_abaster"             
 [57] "Syngnathus_schlegeli"            "Hippocampus_comes"              
 [59] "Hippocampus_trimaculatus"        "Hippocampus_kuda"               
 [61] "Hippocampus_spinosissimus"       "Hippocampus_abdominalis"        
 [63] "Hippichthys_cyanospilos"         "Syngnathoides_biaculeatus"      
 [65] "Corythoichthys_intestinalis"     "Corythoichthys_haematopterus"   
 [67] "Doryichthys_boaja"               "Doryrhamphus_japonicus"         
 [69] "Entelurus_aequoreus"             "Anguilla_anguilla"              
 [71] "Syncerus_caffer"                 "Tragelaphus_eurycerus"          
 [73] "Tragelaphus_scriptus"            "Redunca_arundinum"              
 [75] "Antidorcas_marsupialis"          "Nanger_granti"                  
 [77] "Raphicerus_campestris"           "Cephalophus_natalensis"         
 [79] "Sylvicapra_grimmia"              "Oreotragus_oreotragus"          
 [81] "Hippotragus_equinus"             "Addax_nasomaculatus"            
 [83] "Oryx_gazella"                    "Connochaetes_taurinus"          
 [85] "Damaliscus_lunatus"              "Ovis_aries"                     
 [87] "Aepyceros_melampus"              "Pteropus_alecto"                
 [89] "Pteropus_poliocephalus"          "Pteropus_scapulatus"            
 [91] "Lasiopodomys_brandtii"           "Meriones_unguiculatus"          
 [93] "Chlorocebus_aethiops"            "Chlorocebus_sabaeus"            
 [95] "Erythrocebus_patas"              "Papio_cynocephalus"             
 [97] "Papio_anubis"                    "Cercocebus_atys"                
 [99] "Macaca_arctoides"                "Macaca_mulatta"                 
[101] "Macaca_maura"                    "Trachypithecus_francoisi"       
[103] "Phrynocephalus_vlangalii"        "Platycercus_caledonicus"        
[105] "Platycercus_elegans"             "Platycercus_eximius"            
[107] "Platycercus_venustus"            "Barnardius_zonarius"            
[109] "Northiella_haematogaster"        "Lathamus_discolor"              
[111] "Neophema_chrysostoma"            "Neophema_pulchella"             
[113] "Psittacula_krameri"              "Eclectus_roratus"               
[115] "Aprosmictus_jonquillaceus"       "Alisterus_scapularis"           
[117] "Alisterus_amboinensis"           "Polytelis_alexandrae"           
[119] "Trichoglossus_haematodus"        "Agapornis_lilianae"             
[121] "Agapornis_taranta"               "Loriculus_vernalis"             
[123] "Ara_ararauna"                    "Guaruba_guarouba"               
[125] "Enicognathus_leptorhynchus"      "Pionites_melanocephalus"        
[127] "Forpus_coelestis"                "Forpus_passerinus"              
[129] "Amazona_leucocephala"            "Amazona_albifrons"              
[131] "Amazona_pretrei"                 "Amazona_vinacea"                
[133] "Amazona_finschi"                 "Amazona_amazonica"              
[135] "Amazona_ochrocephala"            "Amazona_aestiva"                
[137] "Brotogeris_versicolurus"         "Brotogeris_pyrrhoptera"         
[139] "Bolborhynchus_lineola"           "Poicephalus_gulielmi"           
[141] "Poicephalus_meyeri"              "Psittacus_erithacus"            
[143] "Serinus_mennelli"                "Ficedula_hypoleuca"             
[145] "Padda_oryzivora"                 "Lonchura_punctulata"            
[147] "Lonchura_flaviprymna"            "Lonchura_bicolor"               
[149] "Hypargos_niveoguttatus"          "Pytilia_hypogrammica"           
[151] "Uraeginthus_bengalus"            "Amandava_amandava"              
[153] "Plectrophenax_nivalis"           "Serinus_canaria"                
[155] "Serinus_serinus"                 "Loxia_curvirostra"              
[157] "Carduelis_carduelis"             "Carpodacus_roseus"              
[159] "Uragus_sibiricus"                "Pyrrhula_pyrrhula"              
[161] "Pinicola_enucleator"             "Bucanetes_githagineus"          
[163] "Coccothraustes_coccothraustes"   "Mycerobas_affinis"              
[165] "Mycerobas_carnipes"              "Fringilla_montifringilla"       
[167] "Fringilla_coelebs"               "Anas_acuta"                     
[169] "Anas_crecca"                     "Mergus_serrator"                
[171] "Mergus_merganser"                "Bucephala_clangula"             
[173] "Melanitta_nigra"                 "Somateria_mollissima"           
[175] "Callonetta_leucophrys"           "Branta_bernicla"                
[177] "Branta_leucopsis"                "Anser_anser"                    
[179] "Cygnus_columbianus"              "Tadorna_tadorna"                
[181] "Pucrasia_macrolopha"             "Phasianus_colchicus"            
[183] "Chrysolophus_pictus"             "Lophura_nycthemera"             
[185] "Lophura_ignita"                  "Syrmaticus_ellioti"             
[187] "Syrmaticus_reevesii"             "Perdix_perdix"                  
[189] "Tetrao_urogallus"                "Lagopus_muta"                   
[191] "Lagopus_lagopus"                 "Lophophorus_impejanus"          
[193] "Tragopan_temminckii"             "Pavo_cristatus"                 
[195] "Gallus_gallus"                   "Gallus_sonneratii"              
[197] "Rollulus_rouloul"                "Arborophila_torqueola"          
[199] "Numida_meleagris"               

setdiff(dat_organ$phylo, tree$tip.label)      # Database species absent in tree

character(0)

# Exclude mismatched species for consistency
tree_organ <- keep.tip(tree, intersect(tree$tip.label, dat_organ$phylo))
dat_organ <- dat_organ[dat_organ$phylo %in% tree_organ$tip.label, ]

# Re-check for mismatches after pruning
setdiff(tree_organ$tip.label, dat_organ$phylo)

character(0)

setdiff(dat_organ$phylo, tree_organ$tip.label)

character(0)

# Is ultrametric?
is.ultrametric(tree_organ)

[1] TRUE

# how many species?
length(tree_organ$tip.label)

[1] 5

# Compute phylogenetic covariance matrix
A_organ <- ape::vcv.phylo(tree_organ, corr = TRUE)

Model simple

model_simple_adrenal_glands <- brm(
  log10_organ_size ~ 
    1 + 
    log10_body_size +
    (1 + log10_body_size | sex),
  data = dat_organ,
  family = gaussian(),
    prior = c(
    prior(normal(0.75, 0.5), class = "b", coef = "log10_body_size"),
    prior(normal(0, 2), class = "Intercept"),
    prior(student_t(3, 0, 0.5), class = "sd"),
    prior(student_t(3, 0, 1),   class = "sigma")
  ),
  seed = 6955,
  sample_prior = TRUE, 
  chains = 4, 
  cores = 4, 
  threads = threading(2),
  iter = 4000, 
  warmup = 2000,
  control = list(max_treedepth = 15,
                 adapt_delta = 0.99),
  save_pars = save_pars(all = TRUE),
  silent = TRUE, refresh = 0
)

Running MCMC with 4 parallel chains, with 2 thread(s) per chain...

Chain 4 finished in 1.5 seconds.
Chain 1 finished in 1.7 seconds.
Chain 3 finished in 1.8 seconds.
Chain 2 finished in 1.9 seconds.

All 4 chains finished successfully.
Mean chain execution time: 1.7 seconds.
Total execution time: 2.0 seconds.

Model check

pp_check(model_simple_adrenal_glands)

Model summary

summary(model_simple_adrenal_glands)

 Family: gaussian 
  Links: mu = identity 
Formula: log10_organ_size ~ 1 + log10_body_size + (1 + log10_body_size | sex) 
   Data: dat_organ (Number of observations: 10) 
  Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
         total post-warmup draws = 8000

Multilevel Hyperparameters:
~sex (Number of levels: 2) 
                               Estimate Est.Error l-95% CI u-95% CI Rhat
sd(Intercept)                      0.32      0.33     0.01     1.18 1.00
sd(log10_body_size)                0.20      0.22     0.01     0.80 1.00
cor(Intercept,log10_body_size)    -0.11      0.59    -0.98     0.94 1.00
                               Bulk_ESS Tail_ESS
sd(Intercept)                      3565     2537
sd(log10_body_size)                2615     3544
cor(Intercept,log10_body_size)     3227     3452

Regression Coefficients:
                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept          -3.00      0.37    -3.72    -2.24 1.00     4586     4119
log10_body_size     0.74      0.17     0.41     1.11 1.00     2893     2630

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.26      0.09     0.15     0.48 1.00     3971     3995

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Coeficients by sex

# Extract combined coefficients (fixed + random)
coef(model_simple_adrenal_glands)$sex %>%
{tibble(
  sex = dimnames(.)[[1]],
  intercept    = .[, "Estimate", "Intercept"],
  intercept_low  = .[, "Q2.5", "Intercept"], 
  intercept_high = .[, "Q97.5", "Intercept"]
)} %>%
bind_cols(
  tibble(
    sex = dimnames(coef(model_simple_adrenal_glands)$sex)[[1]],
    slope = coef(model_simple_adrenal_glands)$sex[, "Estimate", "log10_body_size"]
  ) %>% select(slope),
  tibble(
    sex = dimnames(coef(model_simple_adrenal_glands)$sex)[[1]],
    slope_low  = coef(model_simple_adrenal_glands)$sex[, "Q2.5", "log10_body_size"],
    slope_high = coef(model_simple_adrenal_glands)$sex[, "Q97.5", "log10_body_size"]
  ) %>% select(slope_low, slope_high)
)

# A tibble: 2 × 7
  sex    intercept intercept_low intercept_high slope slope_low slope_high
  <chr>      <dbl>         <dbl>          <dbl> <dbl>     <dbl>      <dbl>
1 female     -2.99         -3.54          -2.44 0.720     0.508      0.926
2 male       -3.05         -3.59          -2.50 0.727     0.520      0.932

Model with phylogeny

model_phylo_adrenal_glands <- brm(
  log10_organ_size ~ 
    1 + 
    log10_body_size +
    (1 + log10_body_size | sex) +
    (1 | gr(phylo, cov = A_organ)),
  data = dat_organ,
  data2 = list(A_organ = A_organ),
  family = gaussian(),
    prior = c(
    prior(normal(0.75, 0.5), class = "b", coef = "log10_body_size"),
    prior(normal(0, 2), class = "Intercept"),
    prior(student_t(3, 0, 0.5), class = "sd"),
    prior(student_t(3, 0, 1),   class = "sigma")
  ),
  seed = 6955,
  sample_prior = TRUE, 
  chains = 4, 
  cores = 4, 
  threads = threading(2),
  iter = 4000, 
  warmup = 2000,
  control = list(max_treedepth = 15,
                 adapt_delta = 0.99),
  save_pars = save_pars(all = TRUE),
  silent = TRUE, refresh = 0
)

Running MCMC with 4 parallel chains, with 2 thread(s) per chain...

Chain 2 finished in 3.4 seconds.
Chain 3 finished in 3.8 seconds.
Chain 1 finished in 4.0 seconds.
Chain 4 finished in 4.0 seconds.

All 4 chains finished successfully.
Mean chain execution time: 3.8 seconds.
Total execution time: 4.2 seconds.

Model checks

pp_check(model_phylo_adrenal_glands)

Model summary

summary(model_phylo_adrenal_glands)

 Family: gaussian 
  Links: mu = identity 
Formula: log10_organ_size ~ 1 + log10_body_size + (1 + log10_body_size | sex) + (1 | gr(phylo, cov = A_organ)) 
   Data: dat_organ (Number of observations: 10) 
  Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
         total post-warmup draws = 8000

Multilevel Hyperparameters:
~phylo (Number of levels: 5) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.35      0.23     0.02     0.92 1.00     2085     2237

~sex (Number of levels: 2) 
                               Estimate Est.Error l-95% CI u-95% CI Rhat
sd(Intercept)                      0.31      0.32     0.01     1.14 1.00
sd(log10_body_size)                0.19      0.22     0.00     0.80 1.00
cor(Intercept,log10_body_size)    -0.09      0.60    -0.97     0.95 1.00
                               Bulk_ESS Tail_ESS
sd(Intercept)                      4236     3612
sd(log10_body_size)                2443     4133
cor(Intercept,log10_body_size)     5135     4652

Regression Coefficients:
                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept          -3.05      0.47    -4.01    -2.06 1.00     4480     4249
log10_body_size     0.76      0.18     0.41     1.16 1.00     2808     2514

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.18      0.09     0.08     0.39 1.00     2083     3940

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Phylogenetic signal

hyp <- paste(
  "sd_phylo__Intercept^2 /",
  "(sd_phylo__Intercept^2 + sd_sex__Intercept^2 + sigma^2) = 0"
)
(hyp <- hypothesis(model_phylo_adrenal_glands, hyp, class = NULL))

Hypothesis Tests for class :
                Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
1 (sd_phylo__Interc... = 0     0.47      0.32        0     0.96       0.32
  Post.Prob Star
1      0.24    *
---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.

lambda_adrenal_glands <- hyp$hypothesis$Estimate

Coeficients by sex

# Extract combined coefficients (fixed + random)
coef(model_phylo_adrenal_glands)$sex %>%
{tibble(
  sex = dimnames(.)[[1]],
  intercept    = .[, "Estimate", "Intercept"],
  intercept_low  = .[, "Q2.5", "Intercept"], 
  intercept_high = .[, "Q97.5", "Intercept"]
)} %>%
bind_cols(
  tibble(
    sex = dimnames(coef(model_phylo_adrenal_glands)$sex)[[1]],
    slope = coef(model_phylo_adrenal_glands)$sex[, "Estimate", "log10_body_size"]
  ) %>% select(slope),
  tibble(
    sex = dimnames(coef(model_phylo_adrenal_glands)$sex)[[1]],
    slope_low  = coef(model_phylo_adrenal_glands)$sex[, "Q2.5", "log10_body_size"],
    slope_high = coef(model_phylo_adrenal_glands)$sex[, "Q97.5", "log10_body_size"]
  ) %>% select(slope_low, slope_high)
)

# A tibble: 2 × 7
  sex    intercept intercept_low intercept_high slope slope_low slope_high
  <chr>      <dbl>         <dbl>          <dbl> <dbl>     <dbl>      <dbl>
1 female     -3.03         -3.84          -2.22 0.743     0.519      0.968
2 male       -3.11         -3.91          -2.32 0.755     0.537      0.976

Model comparison

loo_simple <- loo(model_simple_adrenal_glands)
loo_phylo <- loo(model_phylo_adrenal_glands)
loo_compare(loo_simple, loo_phylo)

                            elpd_diff se_diff
model_phylo_adrenal_glands   0.0       0.0   
model_simple_adrenal_glands -2.5       0.7   

In this comparison, the model that accounts for phylogeny is clearly preferable to the model that treats species as independent observational units. The phylogenetic model (model_phylo_adrenal glands) is used as the reference with elpd_diff = 0, whereas the simpler model (model_simple_adrenal glands) has an expected log predictive density (ELPD) that is about 2.6 units lower, with an associated standard error of 0.7 In the context of leave‑one‑out cross‑validation, differences exceeding roughly four times their standard error are typically regarded as strong evidence. Here, 2.6/0.7≈3.7, indicating both models are similarly good.

Kidneys

Filter data and prune tree

# Filter organ of interest
dat_organ <- dat_mean %>%
  group_by(organ_grouped, sex) %>%
  ungroup() %>%
  filter(organ_grouped == "kidneys") %>%
  mutate(
    log10_body_size = as.numeric
    (scale(log10_body_size, center = TRUE, scale = FALSE))
  )

# Check for database/tree species mismatches
setdiff(tree$tip.label, dat_organ$phylo)      # Tree species absent in database

  [1] "Petrochromis_famula"             "Lobochilotes_labiatus"          
  [3] "Simochromis_diagramma"           "Petrochromis_orthognathus"      
  [5] "Ctenochromis_horei"              "Limnotilapia_dardennii"         
  [7] "Gnathochromis_pfefferi"          "Tropheus_moorii"                
  [9] "Tropheus_brichardi"              "Haplotaxodon_microlepis"        
 [11] "Perissodus_microlepis"           "Benthochromis_tricoti"          
 [13] "Cyprichromis_leptosoma"          "Greenwoodochromis_christyi"     
 [15] "Limnochromis_staneri"            "Gnathochromis_permaxillaris"    
 [17] "Baileychromis_centropomoides"    "Triglachromis_otostigma"        
 [19] "Cyphotilapia_frontosa"           "Ophthalmotilapia_ventralis"     
 [21] "Cyathopharynx_furcifer"          "Ophthalmotilapia_nasuta"        
 [23] "Ophthalmotilapia_boops"          "Aulonocranus_dewindti"          
 [25] "Ectodus_descampsii"              "Callochromis_melanostigma"      
 [27] "Xenotilapia_melanogenys"         "Xenotilapia_flavipinnis"        
 [29] "Eretmodus_cyanostictus"          "Tanganicodus_irsacae"           
 [31] "Spathodus_marlieri"              "Spathodus_erythrodon"           
 [33] "Neolamprologus_brichardi"        "Julidochromis_regani"           
 [35] "Julidochromis_marlieri"          "Neolamprologus_tetracanthus"    
 [37] "Telmatochromis_temporalis"       "Lepidiolamprologus_elongatus"   
 [39] "Lepidiolamprologus_profundicola" "Lepidiolamprologus_nkambae"     
 [41] "Altolamprologus_compressiceps"   "Neolamprologus_brevis"          
 [43] "Chalinochromis_brichardi"        "Julidochromis_ornatus"          
 [45] "Lamprologus_ornatipinnis"        "Neolamprologus_pulcher"         
 [47] "Variabilichromis_moorii"         "Neolamprologus_tretocephalus"   
 [49] "Neolamprologus_sexfasciatus"     "Trematocara_unimaculatum"       
 [51] "Hemibates_stenosoma"             "Bathybates_ferox"               
 [53] "Oreochromis_niloticus"           "Poecilia_reticulata"            
 [55] "Nerophis_lumbriciformis"         "Syngnathus_abaster"             
 [57] "Syngnathus_schlegeli"            "Hippocampus_comes"              
 [59] "Hippocampus_trimaculatus"        "Hippocampus_kuda"               
 [61] "Hippocampus_spinosissimus"       "Hippocampus_abdominalis"        
 [63] "Hippichthys_cyanospilos"         "Syngnathoides_biaculeatus"      
 [65] "Corythoichthys_intestinalis"     "Corythoichthys_haematopterus"   
 [67] "Doryichthys_boaja"               "Doryrhamphus_japonicus"         
 [69] "Entelurus_aequoreus"             "Anguilla_anguilla"              
 [71] "Syncerus_caffer"                 "Tragelaphus_eurycerus"          
 [73] "Tragelaphus_scriptus"            "Redunca_arundinum"              
 [75] "Antidorcas_marsupialis"          "Nanger_granti"                  
 [77] "Raphicerus_campestris"           "Cephalophus_natalensis"         
 [79] "Sylvicapra_grimmia"              "Oreotragus_oreotragus"          
 [81] "Hippotragus_equinus"             "Addax_nasomaculatus"            
 [83] "Oryx_gazella"                    "Connochaetes_taurinus"          
 [85] "Damaliscus_lunatus"              "Ovis_aries"                     
 [87] "Aepyceros_melampus"              "Pteropus_alecto"                
 [89] "Pteropus_poliocephalus"          "Pteropus_scapulatus"            
 [91] "Chlorocebus_aethiops"            "Chlorocebus_sabaeus"            
 [93] "Erythrocebus_patas"              "Papio_cynocephalus"             
 [95] "Papio_anubis"                    "Cercocebus_atys"                
 [97] "Macaca_arctoides"                "Macaca_mulatta"                 
 [99] "Macaca_maura"                    "Trachypithecus_francoisi"       
[101] "Phrynocephalus_vlangalii"        "Platycercus_caledonicus"        
[103] "Platycercus_elegans"             "Platycercus_eximius"            
[105] "Platycercus_venustus"            "Barnardius_zonarius"            
[107] "Northiella_haematogaster"        "Lathamus_discolor"              
[109] "Neophema_chrysostoma"            "Neophema_pulchella"             
[111] "Psittacula_krameri"              "Eclectus_roratus"               
[113] "Aprosmictus_jonquillaceus"       "Alisterus_scapularis"           
[115] "Alisterus_amboinensis"           "Polytelis_alexandrae"           
[117] "Trichoglossus_haematodus"        "Agapornis_lilianae"             
[119] "Agapornis_taranta"               "Loriculus_vernalis"             
[121] "Ara_ararauna"                    "Guaruba_guarouba"               
[123] "Enicognathus_leptorhynchus"      "Pionites_melanocephalus"        
[125] "Forpus_coelestis"                "Forpus_passerinus"              
[127] "Amazona_leucocephala"            "Amazona_albifrons"              
[129] "Amazona_pretrei"                 "Amazona_vinacea"                
[131] "Amazona_finschi"                 "Amazona_amazonica"              
[133] "Amazona_ochrocephala"            "Amazona_aestiva"                
[135] "Brotogeris_versicolurus"         "Brotogeris_pyrrhoptera"         
[137] "Bolborhynchus_lineola"           "Poicephalus_gulielmi"           
[139] "Poicephalus_meyeri"              "Psittacus_erithacus"            
[141] "Serinus_mennelli"                "Ficedula_hypoleuca"             
[143] "Padda_oryzivora"                 "Lonchura_punctulata"            
[145] "Lonchura_flaviprymna"            "Lonchura_bicolor"               
[147] "Hypargos_niveoguttatus"          "Pytilia_hypogrammica"           
[149] "Uraeginthus_bengalus"            "Amandava_amandava"              
[151] "Plectrophenax_nivalis"           "Serinus_canaria"                
[153] "Serinus_serinus"                 "Loxia_curvirostra"              
[155] "Carduelis_carduelis"             "Carpodacus_roseus"              
[157] "Uragus_sibiricus"                "Pyrrhula_pyrrhula"              
[159] "Pinicola_enucleator"             "Bucanetes_githagineus"          
[161] "Coccothraustes_coccothraustes"   "Mycerobas_affinis"              
[163] "Mycerobas_carnipes"              "Fringilla_montifringilla"       
[165] "Fringilla_coelebs"               "Anas_acuta"                     
[167] "Anas_crecca"                     "Mergus_serrator"                
[169] "Mergus_merganser"                "Bucephala_clangula"             
[171] "Melanitta_nigra"                 "Somateria_mollissima"           
[173] "Callonetta_leucophrys"           "Branta_bernicla"                
[175] "Branta_leucopsis"                "Anser_anser"                    
[177] "Cygnus_columbianus"              "Tadorna_tadorna"                
[179] "Coturnix_japonica"               "Pucrasia_macrolopha"            
[181] "Phasianus_colchicus"             "Chrysolophus_pictus"            
[183] "Lophura_nycthemera"              "Lophura_ignita"                 
[185] "Syrmaticus_ellioti"              "Syrmaticus_reevesii"            
[187] "Perdix_perdix"                   "Tetrao_urogallus"               
[189] "Lagopus_muta"                    "Lagopus_lagopus"                
[191] "Lophophorus_impejanus"           "Tragopan_temminckii"            
[193] "Pavo_cristatus"                  "Gallus_gallus"                  
[195] "Gallus_sonneratii"               "Rollulus_rouloul"               
[197] "Arborophila_torqueola"           "Numida_meleagris"               

setdiff(dat_organ$phylo, tree$tip.label)      # Database species absent in tree

character(0)

# Exclude mismatched species for consistency
tree_organ <- keep.tip(tree, intersect(tree$tip.label, dat_organ$phylo))
dat_organ <- dat_organ[dat_organ$phylo %in% tree_organ$tip.label, ]

# Re-check for mismatches after pruning
setdiff(tree_organ$tip.label, dat_organ$phylo)

character(0)

setdiff(dat_organ$phylo, tree_organ$tip.label)

character(0)

# Is ultrametric?
is.ultrametric(tree_organ)

[1] TRUE

# how many species?
length(tree_organ$tip.label)

[1] 6

# Compute phylogenetic covariance matrix
A_organ <- ape::vcv.phylo(tree_organ, corr = TRUE)

Model simple

model_simple_kidneys <- brm(
  log10_organ_size ~ 
    1 + 
    log10_body_size +
    (1 + log10_body_size | sex),
  data = dat_organ,
  family = gaussian(),
    prior = c(
    prior(normal(0.75, 0.5), class = "b", coef = "log10_body_size"),
    prior(normal(0, 2), class = "Intercept"),
    prior(student_t(3, 0, 0.5), class = "sd"),
    prior(student_t(3, 0, 1),   class = "sigma")
  ),
  seed = 6955,
  sample_prior = TRUE, 
  chains = 4, 
  cores = 4, 
  threads = threading(2),
  iter = 4000, 
  warmup = 2000,
  control = list(max_treedepth = 15,
                 adapt_delta = 0.99),
  save_pars = save_pars(all = TRUE),
  silent = TRUE, refresh = 0
)

Running MCMC with 4 parallel chains, with 2 thread(s) per chain...

Chain 1 finished in 1.8 seconds.
Chain 2 finished in 1.8 seconds.
Chain 3 finished in 2.0 seconds.
Chain 4 finished in 2.7 seconds.

All 4 chains finished successfully.
Mean chain execution time: 2.1 seconds.
Total execution time: 2.9 seconds.

Model checks

pp_check(model_simple_kidneys)

Model summary

summary(model_simple_kidneys)

 Family: gaussian 
  Links: mu = identity 
Formula: log10_organ_size ~ 1 + log10_body_size + (1 + log10_body_size | sex) 
   Data: dat_organ (Number of observations: 12) 
  Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
         total post-warmup draws = 8000

Multilevel Hyperparameters:
~sex (Number of levels: 2) 
                               Estimate Est.Error l-95% CI u-95% CI Rhat
sd(Intercept)                      0.20      0.26     0.00     0.88 1.00
sd(log10_body_size)                0.18      0.21     0.00     0.78 1.00
cor(Intercept,log10_body_size)    -0.01      0.61    -0.97     0.96 1.00
                               Bulk_ESS Tail_ESS
sd(Intercept)                      2158     3146
sd(log10_body_size)                2092     2753
cor(Intercept,log10_body_size)     4904     4636

Regression Coefficients:
                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept           0.27      0.20    -0.18     0.73 1.00     1706     1554
log10_body_size     0.81      0.16     0.43     1.12 1.00     2361     1672

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.09      0.03     0.06     0.16 1.00     4386     3614

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Coeficients by sex

# Extract combined coefficients (fixed + random)
coef(model_simple_kidneys)$sex %>%
{tibble(
  sex = dimnames(.)[[1]],
  intercept    = .[, "Estimate", "Intercept"],
  intercept_low  = .[, "Q2.5", "Intercept"], 
  intercept_high = .[, "Q97.5", "Intercept"]
)} %>%
bind_cols(
  tibble(
    sex = dimnames(coef(model_simple_kidneys)$sex)[[1]],
    slope = coef(model_simple_kidneys)$sex[, "Estimate", "log10_body_size"]
  ) %>% select(slope),
  tibble(
    sex = dimnames(coef(model_simple_kidneys)$sex)[[1]],
    slope_low  = coef(model_simple_kidneys)$sex[, "Q2.5", "log10_body_size"],
    slope_high = coef(model_simple_kidneys)$sex[, "Q97.5", "log10_body_size"]
  ) %>% select(slope_low, slope_high)
)

# A tibble: 2 × 7
  sex    intercept intercept_low intercept_high slope slope_low slope_high
  <chr>      <dbl>         <dbl>          <dbl> <dbl>     <dbl>      <dbl>
1 female     0.263         0.193          0.337 0.827     0.743      0.910
2 male       0.270         0.198          0.342 0.815     0.732      0.895

Model with phylogeny

model_phylo_kidneys <- brm(
  log10_organ_size ~ 
    1 + 
    log10_body_size +
    (1 + log10_body_size | sex) +
    (1 | gr(phylo, cov = A_organ)),
  data = dat_organ,
  data2 = list(A_organ = A_organ),
  family = gaussian(),
    prior = c(
    prior(normal(0.75, 0.5), class = "b", coef = "log10_body_size"),
    prior(normal(0, 2), class = "Intercept"),
    prior(student_t(3, 0, 0.5), class = "sd"),
    prior(student_t(3, 0, 1),   class = "sigma")
  ),
  seed = 6955,
  sample_prior = TRUE, 
  chains = 4, 
  cores = 4, 
  threads = threading(2),
  iter = 4000, 
  warmup = 2000,
  control = list(max_treedepth = 15,
                 adapt_delta = 0.99),
  save_pars = save_pars(all = TRUE),
  silent = TRUE, refresh = 0
)

Running MCMC with 4 parallel chains, with 2 thread(s) per chain...

Chain 3 finished in 12.7 seconds.
Chain 2 finished in 14.2 seconds.
Chain 1 finished in 17.2 seconds.
Chain 4 finished in 18.4 seconds.

All 4 chains finished successfully.
Mean chain execution time: 15.6 seconds.
Total execution time: 18.6 seconds.

Model checks

pp_check(model_phylo_kidneys)

Model summary

summary(model_phylo_kidneys)

 Family: gaussian 
  Links: mu = identity 
Formula: log10_organ_size ~ 1 + log10_body_size + (1 + log10_body_size | sex) + (1 | gr(phylo, cov = A_organ)) 
   Data: dat_organ (Number of observations: 12) 
  Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
         total post-warmup draws = 8000

Multilevel Hyperparameters:
~phylo (Number of levels: 6) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.27      0.14     0.06     0.63 1.00      953     1349

~sex (Number of levels: 2) 
                               Estimate Est.Error l-95% CI u-95% CI Rhat
sd(Intercept)                      0.17      0.24     0.00     0.84 1.00
sd(log10_body_size)                0.15      0.20     0.00     0.70 1.00
cor(Intercept,log10_body_size)    -0.04      0.62    -0.98     0.97 1.00
                               Bulk_ESS Tail_ESS
sd(Intercept)                      2286     3600
sd(log10_body_size)                2586     3691
cor(Intercept,log10_body_size)     5434     4523

Regression Coefficients:
                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept           0.31      0.28    -0.25     0.86 1.00     2934     2758
log10_body_size     0.77      0.14     0.47     1.08 1.00     3390     3035

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.03      0.02     0.01     0.09 1.00      683     1540

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Phylogenetic signal

hyp <- paste(
  "sd_phylo__Intercept^2 /",
  "(sd_phylo__Intercept^2 + sd_sex__Intercept^2 + sigma^2) = 0"
)
(hyp <- hypothesis(model_phylo_kidneys, hyp, class = NULL))

Hypothesis Tests for class :
                Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
1 (sd_phylo__Interc... = 0     0.71      0.32     0.03        1       0.11
  Post.Prob Star
1       0.1    *
---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.

lambda_kidneys <- hyp$hypothesis$Estimate

Coeficients by sex

# Extract combined coefficients (fixed + random)
coef(model_phylo_kidneys)$sex %>%
{tibble(
  sex = dimnames(.)[[1]],
  intercept    = .[, "Estimate", "Intercept"],
  intercept_low  = .[, "Q2.5", "Intercept"], 
  intercept_high = .[, "Q97.5", "Intercept"]
)} %>%
bind_cols(
  tibble(
    sex = dimnames(coef(model_phylo_kidneys)$sex)[[1]],
    slope = coef(model_phylo_kidneys)$sex[, "Estimate", "log10_body_size"]
  ) %>% select(slope),
  tibble(
    sex = dimnames(coef(model_phylo_kidneys)$sex)[[1]],
    slope_low  = coef(model_phylo_kidneys)$sex[, "Q2.5", "log10_body_size"],
    slope_high = coef(model_phylo_kidneys)$sex[, "Q97.5", "log10_body_size"]
  ) %>% select(slope_low, slope_high)
)

# A tibble: 2 × 7
  sex    intercept intercept_low intercept_high slope slope_low slope_high
  <chr>      <dbl>         <dbl>          <dbl> <dbl>     <dbl>      <dbl>
1 female     0.308        -0.113          0.732 0.777     0.657      0.889
2 male       0.318        -0.101          0.744 0.765     0.646      0.878

Model comparison

loo_simple <- loo(model_simple_kidneys)
loo_phylo <- loo(model_phylo_kidneys)
loo_compare(loo_simple, loo_phylo)

                     elpd_diff se_diff
model_phylo_kidneys    0.0       0.0  
model_simple_kidneys -10.6       1.6  

In this comparison, the phylogenetic model is clearly preferable to the model treating species as independent observational units. The phylogenetic model (model_phylo_kidneys) serves as reference with elpd_diff = 0, while the simpler model (model_simple_kidneys) has an expected log predictive density (ELPD) that is about 10.3 units lower, with an associated standard error of 1.5. In leave-one-out cross-validation, differences exceeding approximately four times their standard error indicate strong evidence. Here, 10.3/1.5≈6.9, providing strong evidence favouring the phylogenetic model.

Analysis by Class (Brain only)

In the following analysis, I examine the data without subsetting by family to test whether the patterns observed at the level of vertebrate families are maintained. These analyses are motivated by the paper The Taxon-Level Problem in the Evolution of Mammalian Brain Size: Facts and Artifacts, which shows that allometric slopes can change with taxonomic level and sampling, and may partly reflect statistical artefacts rather than biological processes. Here, I apply the same rationale and additionally assess whether differences between males and females persist when using all available species with brain size data.

# Identify qualifying families (>10 unique species per family)
good_families_all <- dat_mean %>%
  filter(organ_grouped == "brain") %>%  # Brain organ only
  group_by(class, family) %>%
  summarise(n_species = n_distinct(species), .groups = "drop") %>%
  filter(n_species > 10) %>%  # species threshold
  pull(family)

# Filter main dataset - preserves all measurements for qualifying families
dat_brain_all <- dat_mean %>%
  filter(organ_grouped == "brain", 
         family %in% good_families_all)

# families per class, order and family
dat_brain_all %>% 
  distinct(class, order, family, species, .keep_all = TRUE) %>% 
  count(class, order, family, name = "n_species_unique") %>% 
  arrange(class, desc(n_species_unique))

# A tibble: 8 × 4
  class     order           family          n_species_unique
  <fct>     <fct>           <fct>                      <int>
1 Aves      Psittaciformes  Psittacidae                   39
2 Aves      Galliformes     Phasianidae                   17
3 Aves      Passeriformes   Fringillidae                  16
4 Aves      Anseriformes    Anatidae                      14
5 Mammalia  Cetartiodactyla Bovidae                       17
6 Mammalia  Primates        Cercopithecidae               11
7 Teleostei Perciformes     Cichlidae                     52
8 Teleostei Syngnathiformes Syngnathidae                  15

Fishes

Filter data and prune tree

# Filter full data to those families only
dat_organ_fishes <- dat_brain_all %>%
         filter(class == "Teleostei") %>%
  ungroup()

# Check for database/tree species mismatches
setdiff(tree$tip.label, dat_organ_fishes$phylo)      # Tree species absent in database

  [1] "Oreochromis_niloticus"         "Poecilia_reticulata"          
  [3] "Anguilla_anguilla"             "Syncerus_caffer"              
  [5] "Tragelaphus_eurycerus"         "Tragelaphus_scriptus"         
  [7] "Redunca_arundinum"             "Antidorcas_marsupialis"       
  [9] "Nanger_granti"                 "Raphicerus_campestris"        
 [11] "Cephalophus_natalensis"        "Sylvicapra_grimmia"           
 [13] "Oreotragus_oreotragus"         "Hippotragus_equinus"          
 [15] "Addax_nasomaculatus"           "Oryx_gazella"                 
 [17] "Connochaetes_taurinus"         "Damaliscus_lunatus"           
 [19] "Ovis_aries"                    "Aepyceros_melampus"           
 [21] "Pteropus_alecto"               "Pteropus_poliocephalus"       
 [23] "Pteropus_scapulatus"           "Lasiopodomys_brandtii"        
 [25] "Rattus_norvegicus"             "Mus_musculus"                 
 [27] "Meriones_unguiculatus"         "Chlorocebus_aethiops"         
 [29] "Chlorocebus_sabaeus"           "Erythrocebus_patas"           
 [31] "Papio_cynocephalus"            "Papio_anubis"                 
 [33] "Cercocebus_atys"               "Macaca_arctoides"             
 [35] "Macaca_mulatta"                "Macaca_fascicularis"          
 [37] "Macaca_maura"                  "Trachypithecus_francoisi"     
 [39] "Phrynocephalus_vlangalii"      "Platycercus_caledonicus"      
 [41] "Platycercus_elegans"           "Platycercus_eximius"          
 [43] "Platycercus_venustus"          "Barnardius_zonarius"          
 [45] "Northiella_haematogaster"      "Lathamus_discolor"            
 [47] "Neophema_chrysostoma"          "Neophema_pulchella"           
 [49] "Psittacula_krameri"            "Eclectus_roratus"             
 [51] "Aprosmictus_jonquillaceus"     "Alisterus_scapularis"         
 [53] "Alisterus_amboinensis"         "Polytelis_alexandrae"         
 [55] "Trichoglossus_haematodus"      "Agapornis_lilianae"           
 [57] "Agapornis_taranta"             "Loriculus_vernalis"           
 [59] "Ara_ararauna"                  "Guaruba_guarouba"             
 [61] "Enicognathus_leptorhynchus"    "Pionites_melanocephalus"      
 [63] "Forpus_coelestis"              "Forpus_passerinus"            
 [65] "Amazona_leucocephala"          "Amazona_albifrons"            
 [67] "Amazona_pretrei"               "Amazona_vinacea"              
 [69] "Amazona_finschi"               "Amazona_amazonica"            
 [71] "Amazona_ochrocephala"          "Amazona_aestiva"              
 [73] "Brotogeris_versicolurus"       "Brotogeris_pyrrhoptera"       
 [75] "Bolborhynchus_lineola"         "Poicephalus_gulielmi"         
 [77] "Poicephalus_meyeri"            "Psittacus_erithacus"          
 [79] "Serinus_mennelli"              "Ficedula_hypoleuca"           
 [81] "Padda_oryzivora"               "Lonchura_punctulata"          
 [83] "Lonchura_flaviprymna"          "Lonchura_bicolor"             
 [85] "Hypargos_niveoguttatus"        "Pytilia_hypogrammica"         
 [87] "Uraeginthus_bengalus"          "Amandava_amandava"            
 [89] "Plectrophenax_nivalis"         "Serinus_canaria"              
 [91] "Serinus_serinus"               "Loxia_curvirostra"            
 [93] "Carduelis_carduelis"           "Carpodacus_roseus"            
 [95] "Uragus_sibiricus"              "Pyrrhula_pyrrhula"            
 [97] "Pinicola_enucleator"           "Bucanetes_githagineus"        
 [99] "Coccothraustes_coccothraustes" "Mycerobas_affinis"            
[101] "Mycerobas_carnipes"            "Fringilla_montifringilla"     
[103] "Fringilla_coelebs"             "Anas_platyrhynchos"           
[105] "Anas_acuta"                    "Anas_crecca"                  
[107] "Mergus_serrator"               "Mergus_merganser"             
[109] "Bucephala_clangula"            "Melanitta_nigra"              
[111] "Somateria_mollissima"          "Callonetta_leucophrys"        
[113] "Branta_bernicla"               "Branta_leucopsis"             
[115] "Anser_anser"                   "Cygnus_columbianus"           
[117] "Tadorna_tadorna"               "Coturnix_japonica"            
[119] "Pucrasia_macrolopha"           "Phasianus_colchicus"          
[121] "Chrysolophus_pictus"           "Lophura_nycthemera"           
[123] "Lophura_ignita"                "Syrmaticus_ellioti"           
[125] "Syrmaticus_reevesii"           "Perdix_perdix"                
[127] "Tetrao_urogallus"              "Lagopus_muta"                 
[129] "Lagopus_lagopus"               "Lophophorus_impejanus"        
[131] "Tragopan_temminckii"           "Pavo_cristatus"               
[133] "Gallus_gallus"                 "Gallus_sonneratii"            
[135] "Rollulus_rouloul"              "Arborophila_torqueola"        
[137] "Numida_meleagris"             

setdiff(dat_organ_fishes$phylo, tree$tip.label)      # Database species absent in tree

character(0)

# Exclude mismatched species for consistency
tree_organ <- keep.tip(tree, intersect(tree$tip.label, dat_organ_fishes$phylo))
dat_organ_fishes <- dat_organ_fishes[dat_organ_fishes$phylo %in% tree_organ$tip.label, ]

# Re-check for mismatches after pruning
setdiff(tree_organ$tip.label, dat_organ_fishes$phylo)

character(0)

setdiff(dat_organ_fishes$phylo, tree_organ$tip.label)

character(0)

# Is ultrametric?
is.ultrametric(tree_organ)

[1] TRUE

# how many species?
length(tree_organ$tip.label)

[1] 67

# Compute phylogenetic covariance matrix
A_organ <- ape::vcv.phylo(tree_organ, corr = TRUE)

Model simple

model_simple_brain_fishes <- brm(
  log10_organ_size ~ 
    1 + 
    log10_body_size +
    (1 + log10_body_size | sex),
  data = dat_organ_fishes,
  family = gaussian(),
    prior = c(
    prior(normal(0.75, 0.5), class = "b", coef = "log10_body_size"),
    prior(normal(0, 2), class = "Intercept"),
    prior(student_t(3, 0, 0.5), class = "sd"),
    prior(student_t(3, 0, 1),   class = "sigma")
  ),
  seed = 6955,
  sample_prior = TRUE, 
  chains = 4, 
  cores = 4, 
  threads = threading(2),
  iter = 4000, 
  warmup = 2000,
  control = list(max_treedepth = 15,
                 adapt_delta = 0.99),
  save_pars = save_pars(all = TRUE),
  silent = TRUE, refresh = 0
)

Running MCMC with 4 parallel chains, with 2 thread(s) per chain...

Chain 2 finished in 5.9 seconds.
Chain 4 finished in 6.2 seconds.
Chain 3 finished in 7.8 seconds.
Chain 1 finished in 8.0 seconds.

All 4 chains finished successfully.
Mean chain execution time: 7.0 seconds.
Total execution time: 8.2 seconds.

Model check

pp_check(model_simple_brain_fishes)

Model summary

summary(model_simple_brain_fishes)

 Family: gaussian 
  Links: mu = identity 
Formula: log10_organ_size ~ 1 + log10_body_size + (1 + log10_body_size | sex) 
   Data: dat_organ_fishes (Number of observations: 134) 
  Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
         total post-warmup draws = 8000

Multilevel Hyperparameters:
~sex (Number of levels: 2) 
                               Estimate Est.Error l-95% CI u-95% CI Rhat
sd(Intercept)                      0.25      0.29     0.01     1.05 1.00
sd(log10_body_size)                0.18      0.21     0.00     0.74 1.00
cor(Intercept,log10_body_size)    -0.07      0.61    -0.97     0.95 1.00
                               Bulk_ESS Tail_ESS
sd(Intercept)                      1339      349
sd(log10_body_size)                2170     2949
cor(Intercept,log10_body_size)     4613     4495

Regression Coefficients:
                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept          -2.06      0.28    -2.54    -1.23 1.01      776      298
log10_body_size     0.81      0.15     0.49     1.14 1.00     1328     2035

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.28      0.02     0.25     0.32 1.00     6122     4482

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Coefficient by sex

# Extract combined coefficients (fixed + random)
coef(model_simple_brain_fishes)$sex %>%
{tibble(
  sex = dimnames(.)[[1]],
  intercept    = .[, "Estimate", "Intercept"],
  intercept_low  = .[, "Q2.5", "Intercept"], 
  intercept_high = .[, "Q97.5", "Intercept"]
)} %>%
bind_cols(
  tibble(
    sex = dimnames(coef(model_simple_brain_fishes)$sex)[[1]],
    slope = coef(model_simple_brain_fishes)$sex[, "Estimate", "log10_body_size"]
  ) %>% select(slope),
  tibble(
    sex = dimnames(coef(model_simple_brain_fishes)$sex)[[1]],
    slope_low  = coef(model_simple_brain_fishes)$sex[, "Q2.5", "log10_body_size"],
    slope_high = coef(model_simple_brain_fishes)$sex[, "Q97.5", "log10_body_size"]
  ) %>% select(slope_low, slope_high)
)

# A tibble: 2 × 7
  sex    intercept intercept_low intercept_high slope slope_low slope_high
  <chr>      <dbl>         <dbl>          <dbl> <dbl>     <dbl>      <dbl>
1 female     -2.07         -2.19          -1.95 0.801     0.688      0.906
2 male       -2.12         -2.25          -2.00 0.825     0.721      0.939

Model with phylogeny

model_phylo_brain_fishes <- brm(
  log10_organ_size ~ 
    1 + 
    log10_body_size +
    (1 + log10_body_size | sex) +
    (1 | gr(phylo, cov = A_organ)),
  data = dat_organ_fishes,
  data2 = list(A_organ = A_organ),
  family = gaussian(),
    prior = c(
    prior(normal(0.75, 0.5), class = "b", coef = "log10_body_size"),
    prior(normal(0, 2), class = "Intercept"),
    prior(student_t(3, 0, 0.5), class = "sd"),
    prior(student_t(3, 0, 1),   class = "sigma")
  ),
  seed = 6955,
  sample_prior = TRUE, 
  chains = 4, 
  cores = 4, 
  threads = threading(2),
  iter = 4000, 
  warmup = 2000,
  control = list(max_treedepth = 15,
                 adapt_delta = 0.99),
  save_pars = save_pars(all = TRUE),
  silent = TRUE, refresh = 0
)

Running MCMC with 4 parallel chains, with 2 thread(s) per chain...

Chain 1 finished in 79.1 seconds.
Chain 2 finished in 83.0 seconds.
Chain 4 finished in 107.1 seconds.
Chain 3 finished in 108.1 seconds.

All 4 chains finished successfully.
Mean chain execution time: 94.3 seconds.
Total execution time: 108.3 seconds.

Model check

pp_check(model_phylo_brain_fishes)

Model summary

summary(model_phylo_brain_fishes)

 Family: gaussian 
  Links: mu = identity 
Formula: log10_organ_size ~ 1 + log10_body_size + (1 + log10_body_size | sex) + (1 | gr(phylo, cov = A_organ)) 
   Data: dat_organ_fishes (Number of observations: 134) 
  Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
         total post-warmup draws = 8000

Multilevel Hyperparameters:
~phylo (Number of levels: 67) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.31      0.04     0.25     0.38 1.00     1410     3027

~sex (Number of levels: 2) 
                               Estimate Est.Error l-95% CI u-95% CI Rhat
sd(Intercept)                      0.18      0.24     0.00     0.87 1.00
sd(log10_body_size)                0.15      0.20     0.00     0.71 1.00
cor(Intercept,log10_body_size)    -0.07      0.62    -0.98     0.96 1.00
                               Bulk_ESS Tail_ESS
sd(Intercept)                      2170     3138
sd(log10_body_size)                1924     3721
cor(Intercept,log10_body_size)     5994     4839

Regression Coefficients:
                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept          -2.07      0.26    -2.56    -1.53 1.00     1544     2554
log10_body_size     0.49      0.14     0.25     0.88 1.00     2285     1980

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.04      0.00     0.03     0.05 1.00     2656     4751

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Phylogenetic signal

hyp <- paste(
  "sd_phylo__Intercept^2 /",
  "(sd_phylo__Intercept^2 + sd_sex__Intercept^2 + sigma^2) = 0"
)
(hyp <- hypothesis(model_phylo_brain_fishes, hyp, class = NULL))

Hypothesis Tests for class :
                Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
1 (sd_phylo__Interc... = 0     0.77      0.27     0.11     0.99       0.02
  Post.Prob Star
1      0.02    *
---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.

lambda_brain <- hyp$hypothesis$Estimate

Coeficients by sex

# Extract combined coefficients (fixed + random)
coef(model_phylo_brain_fishes)$sex %>%
{tibble(
  sex = dimnames(.)[[1]],
  intercept    = .[, "Estimate", "Intercept"],
  intercept_low  = .[, "Q2.5", "Intercept"], 
  intercept_high = .[, "Q97.5", "Intercept"]
)} %>%
bind_cols(
  tibble(
    sex = dimnames(coef(model_phylo_brain_fishes)$sex)[[1]],
    slope = coef(model_phylo_brain_fishes)$sex[, "Estimate", "log10_body_size"]
  ) %>% select(slope),
  tibble(
    sex = dimnames(coef(model_phylo_brain_fishes)$sex)[[1]],
    slope_low  = coef(model_phylo_brain_fishes)$sex[, "Q2.5", "log10_body_size"],
    slope_high = coef(model_phylo_brain_fishes)$sex[, "Q97.5", "log10_body_size"]
  ) %>% select(slope_low, slope_high)
)

# A tibble: 2 × 7
  sex    intercept intercept_low intercept_high slope slope_low slope_high
  <chr>      <dbl>         <dbl>          <dbl> <dbl>     <dbl>      <dbl>
1 female     -2.07         -2.42          -1.72 0.463     0.414      0.513
2 male       -2.09         -2.44          -1.74 0.475     0.426      0.526

Model comparison

loo_simple <- loo(model_simple_brain_fishes)
loo_phylo <- loo(model_phylo_brain_fishes)
loo_compare(loo_simple, loo_phylo)

                          elpd_diff se_diff
model_phylo_brain_fishes     0.0       0.0 
model_simple_brain_fishes -220.0      15.1 

Birds

Filter data and prune tree

# Filter full data to those families only
dat_organ_birds <- dat_brain_all %>%
         filter(class == "Aves") %>%
  ungroup()

# Check for database/tree species mismatches
setdiff(tree$tip.label, dat_organ_birds$phylo)      # Tree species absent in database

  [1] "Petrochromis_famula"             "Lobochilotes_labiatus"          
  [3] "Simochromis_diagramma"           "Petrochromis_orthognathus"      
  [5] "Ctenochromis_horei"              "Limnotilapia_dardennii"         
  [7] "Gnathochromis_pfefferi"          "Tropheus_moorii"                
  [9] "Tropheus_brichardi"              "Haplotaxodon_microlepis"        
 [11] "Perissodus_microlepis"           "Benthochromis_tricoti"          
 [13] "Cyprichromis_leptosoma"          "Greenwoodochromis_christyi"     
 [15] "Limnochromis_staneri"            "Gnathochromis_permaxillaris"    
 [17] "Baileychromis_centropomoides"    "Triglachromis_otostigma"        
 [19] "Cyphotilapia_frontosa"           "Ophthalmotilapia_ventralis"     
 [21] "Cyathopharynx_furcifer"          "Ophthalmotilapia_nasuta"        
 [23] "Ophthalmotilapia_boops"          "Aulonocranus_dewindti"          
 [25] "Ectodus_descampsii"              "Callochromis_melanostigma"      
 [27] "Xenotilapia_melanogenys"         "Xenotilapia_flavipinnis"        
 [29] "Eretmodus_cyanostictus"          "Tanganicodus_irsacae"           
 [31] "Spathodus_marlieri"              "Spathodus_erythrodon"           
 [33] "Neolamprologus_brichardi"        "Julidochromis_regani"           
 [35] "Julidochromis_marlieri"          "Neolamprologus_tetracanthus"    
 [37] "Telmatochromis_temporalis"       "Lepidiolamprologus_elongatus"   
 [39] "Lepidiolamprologus_profundicola" "Lepidiolamprologus_nkambae"     
 [41] "Altolamprologus_compressiceps"   "Neolamprologus_brevis"          
 [43] "Chalinochromis_brichardi"        "Julidochromis_ornatus"          
 [45] "Lamprologus_ornatipinnis"        "Neolamprologus_pulcher"         
 [47] "Variabilichromis_moorii"         "Neolamprologus_tretocephalus"   
 [49] "Neolamprologus_sexfasciatus"     "Trematocara_unimaculatum"       
 [51] "Hemibates_stenosoma"             "Bathybates_ferox"               
 [53] "Oreochromis_niloticus"           "Poecilia_reticulata"            
 [55] "Nerophis_lumbriciformis"         "Syngnathus_abaster"             
 [57] "Syngnathus_schlegeli"            "Hippocampus_comes"              
 [59] "Hippocampus_trimaculatus"        "Hippocampus_kuda"               
 [61] "Hippocampus_spinosissimus"       "Hippocampus_abdominalis"        
 [63] "Hippichthys_cyanospilos"         "Syngnathoides_biaculeatus"      
 [65] "Corythoichthys_intestinalis"     "Corythoichthys_haematopterus"   
 [67] "Doryichthys_boaja"               "Doryrhamphus_japonicus"         
 [69] "Entelurus_aequoreus"             "Anguilla_anguilla"              
 [71] "Syncerus_caffer"                 "Tragelaphus_eurycerus"          
 [73] "Tragelaphus_scriptus"            "Redunca_arundinum"              
 [75] "Antidorcas_marsupialis"          "Nanger_granti"                  
 [77] "Raphicerus_campestris"           "Cephalophus_natalensis"         
 [79] "Sylvicapra_grimmia"              "Oreotragus_oreotragus"          
 [81] "Hippotragus_equinus"             "Addax_nasomaculatus"            
 [83] "Oryx_gazella"                    "Connochaetes_taurinus"          
 [85] "Damaliscus_lunatus"              "Ovis_aries"                     
 [87] "Aepyceros_melampus"              "Pteropus_alecto"                
 [89] "Pteropus_poliocephalus"          "Pteropus_scapulatus"            
 [91] "Lasiopodomys_brandtii"           "Rattus_norvegicus"              
 [93] "Mus_musculus"                    "Meriones_unguiculatus"          
 [95] "Chlorocebus_aethiops"            "Chlorocebus_sabaeus"            
 [97] "Erythrocebus_patas"              "Papio_cynocephalus"             
 [99] "Papio_anubis"                    "Cercocebus_atys"                
[101] "Macaca_arctoides"                "Macaca_mulatta"                 
[103] "Macaca_fascicularis"             "Macaca_maura"                   
[105] "Trachypithecus_francoisi"        "Phrynocephalus_vlangalii"       
[107] "Ficedula_hypoleuca"              "Padda_oryzivora"                
[109] "Lonchura_punctulata"             "Lonchura_flaviprymna"           
[111] "Lonchura_bicolor"                "Hypargos_niveoguttatus"         
[113] "Pytilia_hypogrammica"            "Uraeginthus_bengalus"           
[115] "Amandava_amandava"               "Coturnix_japonica"              
[117] "Gallus_gallus"                   "Numida_meleagris"               

setdiff(dat_organ_birds$phylo, tree$tip.label)      # Database species absent in tree

character(0)

# Exclude mismatched species for consistency
tree_organ <- keep.tip(tree, intersect(tree$tip.label, dat_organ_birds$phylo))
dat_organ_birds <- dat_organ_birds[dat_organ_birds$phylo %in% tree_organ$tip.label, ]

# Re-check for mismatches after pruning
setdiff(tree_organ$tip.label, dat_organ_birds$phylo)

character(0)

setdiff(dat_organ_birds$phylo, tree_organ$tip.label)

character(0)

# Is ultrametric?
is.ultrametric(tree_organ)

[1] TRUE

# how many species?
length(tree_organ$tip.label)

[1] 86

# Compute phylogenetic covariance matrix
A_organ <- ape::vcv.phylo(tree_organ, corr = TRUE)

Model simple

model_simple_brain_birds <- brm(
  log10_organ_size ~ 
    1 + 
    log10_body_size +
    (1 + log10_body_size | sex),
  data = dat_organ_birds,
  family = gaussian(),
    prior = c(
    prior(normal(0.75, 0.5), class = "b", coef = "log10_body_size"),
    prior(normal(0, 2), class = "Intercept"),
    prior(student_t(3, 0, 0.5), class = "sd"),
    prior(student_t(3, 0, 1),   class = "sigma")
  ),
  seed = 6955,
  sample_prior = TRUE, 
  chains = 4, 
  cores = 4, 
  threads = threading(2),
  iter = 4000, 
  warmup = 2000,
  control = list(max_treedepth = 15,
                 adapt_delta = 0.99),
  save_pars = save_pars(all = TRUE),
  silent = TRUE, refresh = 0
)

Running MCMC with 4 parallel chains, with 2 thread(s) per chain...

Chain 4 finished in 17.1 seconds.
Chain 1 finished in 18.3 seconds.
Chain 3 finished in 18.2 seconds.
Chain 2 finished in 20.3 seconds.

All 4 chains finished successfully.
Mean chain execution time: 18.5 seconds.
Total execution time: 20.5 seconds.

Model check

pp_check(model_simple_brain_birds)

Model summary

summary(model_simple_brain_birds)

 Family: gaussian 
  Links: mu = identity 
Formula: log10_organ_size ~ 1 + log10_body_size + (1 + log10_body_size | sex) 
   Data: dat_organ_birds (Number of observations: 172) 
  Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
         total post-warmup draws = 8000

Multilevel Hyperparameters:
~sex (Number of levels: 2) 
                               Estimate Est.Error l-95% CI u-95% CI Rhat
sd(Intercept)                      0.24      0.29     0.01     1.02 1.00
sd(log10_body_size)                0.16      0.20     0.00     0.73 1.00
cor(Intercept,log10_body_size)    -0.07      0.61    -0.98     0.96 1.00
                               Bulk_ESS Tail_ESS
sd(Intercept)                      2521     3152
sd(log10_body_size)                1636     3011
cor(Intercept,log10_body_size)     3818     4309

Regression Coefficients:
                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept          -0.52      0.24    -1.01    -0.04 1.00     3017     2986
log10_body_size     0.46      0.13     0.22     0.83 1.00     1707     1106

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.20      0.01     0.18     0.23 1.00     5545     4217

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Coefficient by sex

# Extract combined coefficients (fixed + random)
coef(model_simple_brain_birds)$sex %>%
{tibble(
  sex = dimnames(.)[[1]],
  intercept    = .[, "Estimate", "Intercept"],
  intercept_low  = .[, "Q2.5", "Intercept"], 
  intercept_high = .[, "Q97.5", "Intercept"]
)} %>%
bind_cols(
  tibble(
    sex = dimnames(coef(model_simple_brain_birds)$sex)[[1]],
    slope = coef(model_simple_brain_birds)$sex[, "Estimate", "log10_body_size"]
  ) %>% select(slope),
  tibble(
    sex = dimnames(coef(model_simple_brain_birds)$sex)[[1]],
    slope_low  = coef(model_simple_brain_birds)$sex[, "Q2.5", "log10_body_size"],
    slope_high = coef(model_simple_brain_birds)$sex[, "Q97.5", "log10_body_size"]
  ) %>% select(slope_low, slope_high)
)

# A tibble: 2 × 7
  sex    intercept intercept_low intercept_high slope slope_low slope_high
  <chr>      <dbl>         <dbl>          <dbl> <dbl>     <dbl>      <dbl>
1 female    -0.528        -0.664         -0.399 0.440     0.385      0.499
2 male      -0.498        -0.625         -0.364 0.432     0.376      0.486

Model with phylogeny

model_phylo_brain_birds <- brm(
  log10_organ_size ~ 
    1 + 
    log10_body_size +
    (1 + log10_body_size | sex) +
    (1 | gr(phylo, cov = A_organ)),
  data = dat_organ_birds,
  data2 = list(A_organ = A_organ),
  family = gaussian(),
    prior = c(
    prior(normal(0.75, 0.5), class = "b", coef = "log10_body_size"),
    prior(normal(0, 2), class = "Intercept"),
    prior(student_t(3, 0, 0.5), class = "sd"),
    prior(student_t(3, 0, 1),   class = "sigma")
  ),
  seed = 6955,
  sample_prior = TRUE, 
  chains = 4, 
  cores = 4, 
  threads = threading(2),
  iter = 4000, 
  warmup = 2000,
  control = list(max_treedepth = 15,
                 adapt_delta = 0.99),
  save_pars = save_pars(all = TRUE),
  silent = TRUE, refresh = 0
)

Running MCMC with 4 parallel chains, with 2 thread(s) per chain...

Chain 4 finished in 115.6 seconds.
Chain 3 finished in 134.7 seconds.
Chain 2 finished in 142.1 seconds.
Chain 1 finished in 143.0 seconds.

All 4 chains finished successfully.
Mean chain execution time: 133.8 seconds.
Total execution time: 143.1 seconds.

Model check

pp_check(model_phylo_brain_birds)

Model summary

summary(model_phylo_brain_birds)

 Family: gaussian 
  Links: mu = identity 
Formula: log10_organ_size ~ 1 + log10_body_size + (1 + log10_body_size | sex) + (1 | gr(phylo, cov = A_organ)) 
   Data: dat_organ_birds (Number of observations: 172) 
  Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
         total post-warmup draws = 8000

Multilevel Hyperparameters:
~phylo (Number of levels: 86) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.17      0.02     0.13     0.21 1.00     1348     2941

~sex (Number of levels: 2) 
                               Estimate Est.Error l-95% CI u-95% CI Rhat
sd(Intercept)                      0.21      0.26     0.00     0.91 1.00
sd(log10_body_size)                0.13      0.18     0.00     0.62 1.00
cor(Intercept,log10_body_size)    -0.08      0.62    -0.98     0.95 1.00
                               Bulk_ESS Tail_ESS
sd(Intercept)                      2438     2361
sd(log10_body_size)                1572     2801
cor(Intercept,log10_body_size)     5408     5089

Regression Coefficients:
                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept          -0.82      0.24    -1.33    -0.36 1.00     2569     3258
log10_body_size     0.54      0.12     0.27     0.82 1.00     1339     1326

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.05      0.00     0.04     0.06 1.00     2014     4436

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Phylogenetic signal

hyp <- paste(
  "sd_phylo__Intercept^2 /",
  "(sd_phylo__Intercept^2 + sd_sex__Intercept^2 + sigma^2) = 0"
)
(hyp <- hypothesis(model_phylo_brain_birds, hyp, class = NULL))

Hypothesis Tests for class :
                Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
1 (sd_phylo__Interc... = 0     0.56      0.31     0.03     0.93       0.15
  Post.Prob Star
1      0.13    *
---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.

lambda_brain <- hyp$hypothesis$Estimate

Coeficients by sex

# Extract combined coefficients (fixed + random)
coef(model_phylo_brain_birds)$sex %>%
{tibble(
  sex = dimnames(.)[[1]],
  intercept    = .[, "Estimate", "Intercept"],
  intercept_low  = .[, "Q2.5", "Intercept"], 
  intercept_high = .[, "Q97.5", "Intercept"]
)} %>%
bind_cols(
  tibble(
    sex = dimnames(coef(model_phylo_brain_birds)$sex)[[1]],
    slope = coef(model_phylo_brain_birds)$sex[, "Estimate", "log10_body_size"]
  ) %>% select(slope),
  tibble(
    sex = dimnames(coef(model_phylo_brain_birds)$sex)[[1]],
    slope_low  = coef(model_phylo_brain_birds)$sex[, "Q2.5", "log10_body_size"],
    slope_high = coef(model_phylo_brain_birds)$sex[, "Q97.5", "log10_body_size"]
  ) %>% select(slope_low, slope_high)
)

# A tibble: 2 × 7
  sex    intercept intercept_low intercept_high slope slope_low slope_high
  <chr>      <dbl>         <dbl>          <dbl> <dbl>     <dbl>      <dbl>
1 female    -0.832         -1.04         -0.596 0.537     0.466      0.601
2 male      -0.802         -1.01         -0.564 0.527     0.458      0.590

Model comparison

loo_simple <- loo(model_simple_brain_birds)
loo_phylo <- loo(model_phylo_brain_birds)
loo_compare(loo_simple, loo_phylo)

                         elpd_diff se_diff
model_phylo_brain_birds     0.0       0.0 
model_simple_brain_birds -198.3      12.0 

Mammals

Filter data and prune tree

# Filter full data to those families only
dat_organ_mammals <- dat_brain_all %>%
         filter(class == "Mammalia") %>%
  ungroup()

# Check for database/tree species mismatches
setdiff(tree$tip.label, dat_organ_mammals$phylo)      # Tree species absent in database

  [1] "Petrochromis_famula"             "Lobochilotes_labiatus"          
  [3] "Simochromis_diagramma"           "Petrochromis_orthognathus"      
  [5] "Ctenochromis_horei"              "Limnotilapia_dardennii"         
  [7] "Gnathochromis_pfefferi"          "Tropheus_moorii"                
  [9] "Tropheus_brichardi"              "Haplotaxodon_microlepis"        
 [11] "Perissodus_microlepis"           "Benthochromis_tricoti"          
 [13] "Cyprichromis_leptosoma"          "Greenwoodochromis_christyi"     
 [15] "Limnochromis_staneri"            "Gnathochromis_permaxillaris"    
 [17] "Baileychromis_centropomoides"    "Triglachromis_otostigma"        
 [19] "Cyphotilapia_frontosa"           "Ophthalmotilapia_ventralis"     
 [21] "Cyathopharynx_furcifer"          "Ophthalmotilapia_nasuta"        
 [23] "Ophthalmotilapia_boops"          "Aulonocranus_dewindti"          
 [25] "Ectodus_descampsii"              "Callochromis_melanostigma"      
 [27] "Xenotilapia_melanogenys"         "Xenotilapia_flavipinnis"        
 [29] "Eretmodus_cyanostictus"          "Tanganicodus_irsacae"           
 [31] "Spathodus_marlieri"              "Spathodus_erythrodon"           
 [33] "Neolamprologus_brichardi"        "Julidochromis_regani"           
 [35] "Julidochromis_marlieri"          "Neolamprologus_tetracanthus"    
 [37] "Telmatochromis_temporalis"       "Lepidiolamprologus_elongatus"   
 [39] "Lepidiolamprologus_profundicola" "Lepidiolamprologus_nkambae"     
 [41] "Altolamprologus_compressiceps"   "Neolamprologus_brevis"          
 [43] "Chalinochromis_brichardi"        "Julidochromis_ornatus"          
 [45] "Lamprologus_ornatipinnis"        "Neolamprologus_pulcher"         
 [47] "Variabilichromis_moorii"         "Neolamprologus_tretocephalus"   
 [49] "Neolamprologus_sexfasciatus"     "Trematocara_unimaculatum"       
 [51] "Hemibates_stenosoma"             "Bathybates_ferox"               
 [53] "Oreochromis_niloticus"           "Poecilia_reticulata"            
 [55] "Nerophis_lumbriciformis"         "Syngnathus_abaster"             
 [57] "Syngnathus_schlegeli"            "Hippocampus_comes"              
 [59] "Hippocampus_trimaculatus"        "Hippocampus_kuda"               
 [61] "Hippocampus_spinosissimus"       "Hippocampus_abdominalis"        
 [63] "Hippichthys_cyanospilos"         "Syngnathoides_biaculeatus"      
 [65] "Corythoichthys_intestinalis"     "Corythoichthys_haematopterus"   
 [67] "Doryichthys_boaja"               "Doryrhamphus_japonicus"         
 [69] "Entelurus_aequoreus"             "Anguilla_anguilla"              
 [71] "Pteropus_alecto"                 "Pteropus_poliocephalus"         
 [73] "Pteropus_scapulatus"             "Lasiopodomys_brandtii"          
 [75] "Rattus_norvegicus"               "Mus_musculus"                   
 [77] "Meriones_unguiculatus"           "Phrynocephalus_vlangalii"       
 [79] "Platycercus_caledonicus"         "Platycercus_elegans"            
 [81] "Platycercus_eximius"             "Platycercus_venustus"           
 [83] "Barnardius_zonarius"             "Northiella_haematogaster"       
 [85] "Lathamus_discolor"               "Neophema_chrysostoma"           
 [87] "Neophema_pulchella"              "Psittacula_krameri"             
 [89] "Eclectus_roratus"                "Aprosmictus_jonquillaceus"      
 [91] "Alisterus_scapularis"            "Alisterus_amboinensis"          
 [93] "Polytelis_alexandrae"            "Trichoglossus_haematodus"       
 [95] "Agapornis_lilianae"              "Agapornis_taranta"              
 [97] "Loriculus_vernalis"              "Ara_ararauna"                   
 [99] "Guaruba_guarouba"                "Enicognathus_leptorhynchus"     
[101] "Pionites_melanocephalus"         "Forpus_coelestis"               
[103] "Forpus_passerinus"               "Amazona_leucocephala"           
[105] "Amazona_albifrons"               "Amazona_pretrei"                
[107] "Amazona_vinacea"                 "Amazona_finschi"                
[109] "Amazona_amazonica"               "Amazona_ochrocephala"           
[111] "Amazona_aestiva"                 "Brotogeris_versicolurus"        
[113] "Brotogeris_pyrrhoptera"          "Bolborhynchus_lineola"          
[115] "Poicephalus_gulielmi"            "Poicephalus_meyeri"             
[117] "Psittacus_erithacus"             "Serinus_mennelli"               
[119] "Ficedula_hypoleuca"              "Padda_oryzivora"                
[121] "Lonchura_punctulata"             "Lonchura_flaviprymna"           
[123] "Lonchura_bicolor"                "Hypargos_niveoguttatus"         
[125] "Pytilia_hypogrammica"            "Uraeginthus_bengalus"           
[127] "Amandava_amandava"               "Plectrophenax_nivalis"          
[129] "Serinus_canaria"                 "Serinus_serinus"                
[131] "Loxia_curvirostra"               "Carduelis_carduelis"            
[133] "Carpodacus_roseus"               "Uragus_sibiricus"               
[135] "Pyrrhula_pyrrhula"               "Pinicola_enucleator"            
[137] "Bucanetes_githagineus"           "Coccothraustes_coccothraustes"  
[139] "Mycerobas_affinis"               "Mycerobas_carnipes"             
[141] "Fringilla_montifringilla"        "Fringilla_coelebs"              
[143] "Anas_platyrhynchos"              "Anas_acuta"                     
[145] "Anas_crecca"                     "Mergus_serrator"                
[147] "Mergus_merganser"                "Bucephala_clangula"             
[149] "Melanitta_nigra"                 "Somateria_mollissima"           
[151] "Callonetta_leucophrys"           "Branta_bernicla"                
[153] "Branta_leucopsis"                "Anser_anser"                    
[155] "Cygnus_columbianus"              "Tadorna_tadorna"                
[157] "Coturnix_japonica"               "Pucrasia_macrolopha"            
[159] "Phasianus_colchicus"             "Chrysolophus_pictus"            
[161] "Lophura_nycthemera"              "Lophura_ignita"                 
[163] "Syrmaticus_ellioti"              "Syrmaticus_reevesii"            
[165] "Perdix_perdix"                   "Tetrao_urogallus"               
[167] "Lagopus_muta"                    "Lagopus_lagopus"                
[169] "Lophophorus_impejanus"           "Tragopan_temminckii"            
[171] "Pavo_cristatus"                  "Gallus_gallus"                  
[173] "Gallus_sonneratii"               "Rollulus_rouloul"               
[175] "Arborophila_torqueola"           "Numida_meleagris"               

setdiff(dat_organ_mammals$phylo, tree$tip.label)      # Database species absent in tree

character(0)

# Exclude mismatched species for consistency
tree_organ <- keep.tip(tree, intersect(tree$tip.label, dat_organ_mammals$phylo))
dat_organ_mammals <- dat_organ_mammals[dat_organ_mammals$phylo %in% tree_organ$tip.label, ]

# Re-check for mismatches after pruning
setdiff(tree_organ$tip.label, dat_organ_mammals$phylo)

character(0)

setdiff(dat_organ_mammals$phylo, tree_organ$tip.label)

character(0)

# Is ultrametric?
is.ultrametric(tree_organ)

[1] TRUE

# how many species?
length(tree_organ$tip.label)

[1] 28

# Compute phylogenetic covariance matrix
A_organ <- ape::vcv.phylo(tree_organ, corr = TRUE)

Model simple

model_simple_brain_mammals <- brm(
  log10_organ_size ~ 
    1 + 
    log10_body_size +
    (1 + log10_body_size | sex),
  data = dat_organ_mammals,
  family = gaussian(),
    prior = c(
    prior(normal(0.75, 0.5), class = "b", coef = "log10_body_size"),
    prior(normal(0, 2), class = "Intercept"),
    prior(student_t(3, 0, 0.5), class = "sd"),
    prior(student_t(3, 0, 1),   class = "sigma")
  ),
  seed = 6955,
  sample_prior = TRUE, 
  chains = 4, 
  cores = 4, 
  threads = threading(2),
  iter = 4000, 
  warmup = 2000,
  control = list(max_treedepth = 15,
                 adapt_delta = 0.99),
  save_pars = save_pars(all = TRUE),
  silent = TRUE, refresh = 0
)

Running MCMC with 4 parallel chains, with 2 thread(s) per chain...

Chain 3 finished in 8.5 seconds.
Chain 4 finished in 9.0 seconds.
Chain 2 finished in 10.8 seconds.
Chain 1 finished in 13.0 seconds.

All 4 chains finished successfully.
Mean chain execution time: 10.3 seconds.
Total execution time: 13.0 seconds.

Model check

pp_check(model_simple_brain_mammals)

Model summary

summary(model_simple_brain_mammals)

 Family: gaussian 
  Links: mu = identity 
Formula: log10_organ_size ~ 1 + log10_body_size + (1 + log10_body_size | sex) 
   Data: dat_organ_mammals (Number of observations: 56) 
  Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
         total post-warmup draws = 8000

Multilevel Hyperparameters:
~sex (Number of levels: 2) 
                               Estimate Est.Error l-95% CI u-95% CI Rhat
sd(Intercept)                      0.28      0.32     0.01     1.12 1.00
sd(log10_body_size)                0.14      0.17     0.00     0.62 1.00
cor(Intercept,log10_body_size)    -0.08      0.61    -0.98     0.95 1.00
                               Bulk_ESS Tail_ESS
sd(Intercept)                      2616     3166
sd(log10_body_size)                1337     3038
cor(Intercept,log10_body_size)     4740     4246

Regression Coefficients:
                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept           0.46      0.30    -0.20     1.00 1.00     4101     3938
log10_body_size     0.37      0.11     0.11     0.59 1.00     2148     1411

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.15      0.01     0.12     0.18 1.00     5994     4451

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Coefficient by sex

# Extract combined coefficients (fixed + random)
coef(model_simple_brain_mammals)$sex %>%
{tibble(
  sex = dimnames(.)[[1]],
  intercept    = .[, "Estimate", "Intercept"],
  intercept_low  = .[, "Q2.5", "Intercept"], 
  intercept_high = .[, "Q97.5", "Intercept"]
)} %>%
bind_cols(
  tibble(
    sex = dimnames(coef(model_simple_brain_mammals)$sex)[[1]],
    slope = coef(model_simple_brain_mammals)$sex[, "Estimate", "log10_body_size"]
  ) %>% select(slope),
  tibble(
    sex = dimnames(coef(model_simple_brain_mammals)$sex)[[1]],
    slope_low  = coef(model_simple_brain_mammals)$sex[, "Q2.5", "log10_body_size"],
    slope_high = coef(model_simple_brain_mammals)$sex[, "Q97.5", "log10_body_size"]
  ) %>% select(slope_low, slope_high)
)

# A tibble: 2 × 7
  sex    intercept intercept_low intercept_high slope slope_low slope_high
  <chr>      <dbl>         <dbl>          <dbl> <dbl>     <dbl>      <dbl>
1 female     0.506         0.212          0.800 0.373     0.306      0.439
2 male       0.475         0.166          0.780 0.379     0.310      0.450

Model with phylogeny

model_phylo_brain_mammals <- brm(
  log10_organ_size ~ 
    1 + 
    log10_body_size +
    (1 + log10_body_size | sex) +
    (1 | gr(phylo, cov = A_organ)),
  data = dat_organ_mammals,
  data2 = list(A_organ = A_organ),
  family = gaussian(),
    prior = c(
    prior(normal(0.75, 0.5), class = "b", coef = "log10_body_size"),
    prior(normal(0, 2), class = "Intercept"),
    prior(student_t(3, 0, 0.5), class = "sd"),
    prior(student_t(3, 0, 1),   class = "sigma")
  ),
  seed = 6955,
  sample_prior = TRUE, 
  chains = 4, 
  cores = 4, 
  threads = threading(2),
  iter = 4000, 
  warmup = 2000,
  control = list(max_treedepth = 15,
                 adapt_delta = 0.99),
  save_pars = save_pars(all = TRUE),
  silent = TRUE, refresh = 0
)

Running MCMC with 4 parallel chains, with 2 thread(s) per chain...

Chain 2 finished in 52.1 seconds.
Chain 4 finished in 53.9 seconds.
Chain 3 finished in 58.9 seconds.
Chain 1 finished in 59.4 seconds.

All 4 chains finished successfully.
Mean chain execution time: 56.1 seconds.
Total execution time: 59.5 seconds.

Model check

pp_check(model_phylo_brain_mammals)

Model summary

summary(model_phylo_brain_mammals)

 Family: gaussian 
  Links: mu = identity 
Formula: log10_organ_size ~ 1 + log10_body_size + (1 + log10_body_size | sex) + (1 | gr(phylo, cov = A_organ)) 
   Data: dat_organ_mammals (Number of observations: 56) 
  Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
         total post-warmup draws = 8000

Multilevel Hyperparameters:
~phylo (Number of levels: 28) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.26      0.05     0.18     0.38 1.00     1069     2283

~sex (Number of levels: 2) 
                               Estimate Est.Error l-95% CI u-95% CI Rhat
sd(Intercept)                      0.22      0.28     0.00     0.98 1.00
sd(log10_body_size)                0.10      0.15     0.00     0.56 1.00
cor(Intercept,log10_body_size)    -0.08      0.62    -0.98     0.96 1.00
                               Bulk_ESS Tail_ESS
sd(Intercept)                      2771     3512
sd(log10_body_size)                1637     3363
cor(Intercept,log10_body_size)     4927     4630

Regression Coefficients:
                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept           0.12      0.38    -0.62     0.87 1.00     1725     2838
log10_body_size     0.44      0.10     0.20     0.65 1.00     1508     1513

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.04      0.01     0.03     0.06 1.00      987     2392

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Phylogenetic signal

hyp <- paste(
  "sd_phylo__Intercept^2 /",
  "(sd_phylo__Intercept^2 + sd_sex__Intercept^2 + sigma^2) = 0"
)
(hyp <- hypothesis(model_phylo_brain_mammals, hyp, class = NULL))

Hypothesis Tests for class :
                Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
1 (sd_phylo__Interc... = 0     0.67       0.3     0.06     0.99       0.06
  Post.Prob Star
1      0.06    *
---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.

lambda_brain <- hyp$hypothesis$Estimate

Coeficients by sex

# Extract combined coefficients (fixed + random)
coef(model_phylo_brain_mammals)$sex %>%
{tibble(
  sex = dimnames(.)[[1]],
  intercept    = .[, "Estimate", "Intercept"],
  intercept_low  = .[, "Q2.5", "Intercept"], 
  intercept_high = .[, "Q97.5", "Intercept"]
)} %>%
bind_cols(
  tibble(
    sex = dimnames(coef(model_phylo_brain_mammals)$sex)[[1]],
    slope = coef(model_phylo_brain_mammals)$sex[, "Estimate", "log10_body_size"]
  ) %>% select(slope),
  tibble(
    sex = dimnames(coef(model_phylo_brain_mammals)$sex)[[1]],
    slope_low  = coef(model_phylo_brain_mammals)$sex[, "Q2.5", "log10_body_size"],
    slope_high = coef(model_phylo_brain_mammals)$sex[, "Q97.5", "log10_body_size"]
  ) %>% select(slope_low, slope_high)
)

# A tibble: 2 × 7
  sex    intercept intercept_low intercept_high slope slope_low slope_high
  <chr>      <dbl>         <dbl>          <dbl> <dbl>     <dbl>      <dbl>
1 female     0.159        -0.389          0.808 0.451     0.330      0.557
2 male       0.133        -0.424          0.795 0.456     0.331      0.565

Model comparison

loo_simple <- loo(model_simple_brain_mammals)
loo_phylo <- loo(model_phylo_brain_mammals)
loo_compare(loo_simple, loo_phylo)

                           elpd_diff se_diff
model_phylo_brain_mammals    0.0       0.0  
model_simple_brain_mammals -58.1       6.3  

Across all vertebrates (Brain only)

Filter data and prune tree

good_families_all <- dat_mean %>%
  filter(organ_grouped == "brain") %>%  # Brain organ only
  group_by(class, family) %>%
  summarise(n_species = n_distinct(species), .groups = "drop") %>%
  filter(n_species > 10) %>%  # species threshold
  pull(family)

# Filter main dataset - preserves all measurements for qualifying families
dat_brain_all <- dat_mean %>%
  filter(organ_grouped == "brain", 
         family %in% good_families_all)

# Check for database/tree species mismatches
setdiff(tree$tip.label, dat_brain_all$phylo)      # Tree species absent in database

 [1] "Oreochromis_niloticus"    "Poecilia_reticulata"     
 [3] "Anguilla_anguilla"        "Pteropus_alecto"         
 [5] "Pteropus_poliocephalus"   "Pteropus_scapulatus"     
 [7] "Lasiopodomys_brandtii"    "Rattus_norvegicus"       
 [9] "Mus_musculus"             "Meriones_unguiculatus"   
[11] "Phrynocephalus_vlangalii" "Ficedula_hypoleuca"      
[13] "Padda_oryzivora"          "Lonchura_punctulata"     
[15] "Lonchura_flaviprymna"     "Lonchura_bicolor"        
[17] "Hypargos_niveoguttatus"   "Pytilia_hypogrammica"    
[19] "Uraeginthus_bengalus"     "Amandava_amandava"       
[21] "Coturnix_japonica"        "Gallus_gallus"           
[23] "Numida_meleagris"        

setdiff(dat_brain_all$phylo, tree$tip.label)      # Database species absent in tree

character(0)

# Exclude mismatched species for consistency
tree_organ <- keep.tip(tree, intersect(tree$tip.label, dat_brain_all$phylo))
dat_brain_all <- dat_brain_all[dat_brain_all$phylo %in% tree_organ$tip.label, ]

# Re-check for mismatches after pruning
setdiff(tree_organ$tip.label, dat_brain_all$phylo)

character(0)

setdiff(dat_brain_all$phylo, tree_organ$tip.label)

character(0)

# Is ultrametric?
is.ultrametric(tree_organ)

[1] TRUE

# how many species?
length(tree_organ$tip.label)

[1] 181

# Compute phylogenetic covariance matrix
A_organ <- ape::vcv.phylo(tree_organ, corr = TRUE)

Model simple

model_simple_brain_all <- brm(
  log10_organ_size ~ 
    1 + 
    log10_body_size +
    (1 + log10_body_size | sex),
  data = dat_brain_all,
  family = gaussian(),
    prior = c(
    prior(normal(0.75, 0.5), class = "b", coef = "log10_body_size"),
    prior(normal(0, 2), class = "Intercept"),
    prior(student_t(3, 0, 0.5), class = "sd"),
    prior(student_t(3, 0, 1),   class = "sigma")
  ),
  seed = 6955,
  sample_prior = TRUE, 
  chains = 4, 
  cores = 4, 
  threads = threading(2),
  iter = 4000, 
  warmup = 2000,
  control = list(max_treedepth = 15,
                 adapt_delta = 0.99),
  save_pars = save_pars(all = TRUE),
  silent = TRUE, refresh = 0
)

Running MCMC with 4 parallel chains, with 2 thread(s) per chain...

Chain 4 finished in 20.7 seconds.
Chain 1 finished in 23.3 seconds.
Chain 2 finished in 23.8 seconds.
Chain 3 finished in 45.0 seconds.

All 4 chains finished successfully.
Mean chain execution time: 28.2 seconds.
Total execution time: 45.2 seconds.

Model check

pp_check(model_simple_brain_all)

Model summary

summary(model_simple_brain_all)

 Family: gaussian 
  Links: mu = identity 
Formula: log10_organ_size ~ 1 + log10_body_size + (1 + log10_body_size | sex) 
   Data: dat_brain_all (Number of observations: 362) 
  Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
         total post-warmup draws = 8000

Multilevel Hyperparameters:
~sex (Number of levels: 2) 
                               Estimate Est.Error l-95% CI u-95% CI Rhat
sd(Intercept)                      0.23      0.27     0.01     0.98 1.00
sd(log10_body_size)                0.15      0.19     0.00     0.69 1.00
cor(Intercept,log10_body_size)    -0.03      0.61    -0.97     0.96 1.00
                               Bulk_ESS Tail_ESS
sd(Intercept)                      2665     3415
sd(log10_body_size)                1756     3157
cor(Intercept,log10_body_size)     4013     4461

Regression Coefficients:
                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept          -1.88      0.23    -2.37    -1.41 1.00     3103     3586
log10_body_size     0.90      0.13     0.55     1.15 1.00     1816     1568

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.45      0.02     0.42     0.48 1.00     6418     5369

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Coefficient by sex

# Extract combined coefficients (fixed + random)
coef(model_simple_brain_all)$sex %>%
{tibble(
  sex = dimnames(.)[[1]],
  intercept    = .[, "Estimate", "Intercept"],
  intercept_low  = .[, "Q2.5", "Intercept"], 
  intercept_high = .[, "Q97.5", "Intercept"]
)} %>%
bind_cols(
  tibble(
    sex = dimnames(coef(model_simple_brain_all)$sex)[[1]],
    slope = coef(model_simple_brain_all)$sex[, "Estimate", "log10_body_size"]
  ) %>% select(slope),
  tibble(
    sex = dimnames(coef(model_simple_brain_all)$sex)[[1]],
    slope_low  = coef(model_simple_brain_all)$sex[, "Q2.5", "log10_body_size"],
    slope_high = coef(model_simple_brain_all)$sex[, "Q97.5", "log10_body_size"]
  ) %>% select(slope_low, slope_high)
)

# A tibble: 2 × 7
  sex    intercept intercept_low intercept_high slope slope_low slope_high
  <chr>      <dbl>         <dbl>          <dbl> <dbl>     <dbl>      <dbl>
1 female     -1.86         -1.97          -1.75 0.921     0.878      0.965
2 male       -1.88         -1.99          -1.78 0.918     0.876      0.962

Model with phylogeny

model_phylo_brain_all <- brm(
  log10_organ_size ~ 
    1 + 
    log10_body_size +
    (1 + log10_body_size | sex) +
    (1 | gr(phylo, cov = A_organ)),
  data = dat_brain_all,
  data2 = list(A_organ = A_organ),
  family = gaussian(),
    prior = c(
    prior(normal(0.75, 0.5), class = "b", coef = "log10_body_size"),
    prior(normal(0, 2), class = "Intercept"),
    prior(student_t(3, 0, 0.5), class = "sd"),
    prior(student_t(3, 0, 1),   class = "sigma")
  ),
  seed = 6955,
  sample_prior = TRUE, 
  chains = 4, 
  cores = 4, 
  threads = threading(2),
  iter = 4000, 
  warmup = 2000,
  control = list(max_treedepth = 15,
                 adapt_delta = 0.99),
  save_pars = save_pars(all = TRUE),
  silent = TRUE, refresh = 0
)

Running MCMC with 4 parallel chains, with 2 thread(s) per chain...

Chain 1 finished in 432.9 seconds.
Chain 2 finished in 440.8 seconds.
Chain 3 finished in 497.5 seconds.
Chain 4 finished in 528.8 seconds.

All 4 chains finished successfully.
Mean chain execution time: 475.0 seconds.
Total execution time: 529.0 seconds.

Model check

pp_check(model_phylo_brain_all)

Model summary

summary(model_phylo_brain_all)

 Family: gaussian 
  Links: mu = identity 
Formula: log10_organ_size ~ 1 + log10_body_size + (1 + log10_body_size | sex) + (1 | gr(phylo, cov = A_organ)) 
   Data: dat_brain_all (Number of observations: 362) 
  Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
         total post-warmup draws = 8000

Multilevel Hyperparameters:
~phylo (Number of levels: 181) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.49      0.03     0.42     0.56 1.00     1238     2177

~sex (Number of levels: 2) 
                               Estimate Est.Error l-95% CI u-95% CI Rhat
sd(Intercept)                      0.13      0.22     0.00     0.75 1.00
sd(log10_body_size)                0.09      0.16     0.00     0.57 1.00
cor(Intercept,log10_body_size)    -0.02      0.65    -0.99     0.98 1.00
                               Bulk_ESS Tail_ESS
sd(Intercept)                      2306     4282
sd(log10_body_size)                1474     2307
cor(Intercept,log10_body_size)     6400     4290

Regression Coefficients:
                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept          -1.06      0.33    -1.69    -0.38 1.00     1976     3000
log10_body_size     0.50      0.10     0.31     0.76 1.00     1925     1208

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.05      0.00     0.04     0.05 1.00     2150     4564

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Phylogenetic signal

hyp <- paste(
  "sd_phylo__Intercept^2 /",
  "(sd_phylo__Intercept^2 + sd_sex__Intercept^2 + sigma^2) = 0"
)
(hyp <- hypothesis(model_phylo_brain_all, hyp, class = NULL))

Hypothesis Tests for class :
                Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
1 (sd_phylo__Interc... = 0     0.89      0.18     0.29     0.99          0
  Post.Prob Star
1         0    *
---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.

lambda_brain <- hyp$hypothesis$Estimate

Coeficients by sex

# Extract combined coefficients (fixed + random)
coef(model_phylo_brain_all)$sex %>%
{tibble(
  sex = dimnames(.)[[1]],
  intercept    = .[, "Estimate", "Intercept"],
  intercept_low  = .[, "Q2.5", "Intercept"], 
  intercept_high = .[, "Q97.5", "Intercept"]
)} %>%
bind_cols(
  tibble(
    sex = dimnames(coef(model_phylo_brain_all)$sex)[[1]],
    slope = coef(model_phylo_brain_all)$sex[, "Estimate", "log10_body_size"]
  ) %>% select(slope),
  tibble(
    sex = dimnames(coef(model_phylo_brain_all)$sex)[[1]],
    slope_low  = coef(model_phylo_brain_all)$sex[, "Q2.5", "log10_body_size"],
    slope_high = coef(model_phylo_brain_all)$sex[, "Q97.5", "log10_body_size"]
  ) %>% select(slope_low, slope_high)
)

# A tibble: 2 × 7
  sex    intercept intercept_low intercept_high slope slope_low slope_high
  <chr>      <dbl>         <dbl>          <dbl> <dbl>     <dbl>      <dbl>
1 female     -1.06         -1.62         -0.485 0.486     0.448      0.522
2 male       -1.06         -1.63         -0.490 0.486     0.448      0.523

Model comparison

loo_simple <- loo(model_simple_brain_all)
loo_phylo <- loo(model_phylo_brain_all)
loo_compare(loo_simple, loo_phylo)

                       elpd_diff se_diff
model_phylo_brain_all     0.0       0.0 
model_simple_brain_all -733.1      16.3 

Table S2

The next task is to extract the coefficients from all models (8) in order to prepare a table (Table S2 in the paper) that presents the coefficients by sex, global coefficients. Additionally, I will include the Pagel’s \(\lambda\) representing the phylogenetic signal.

# list of models (no brain included, for now)
models_list_simple <- list(
  model_simple_liver,
  model_simple_heart,
  model_simple_pituitary_glands,
  model_simple_spleen,
  model_simple_stomach,
  model_simple_lungs,
  model_simple_adrenal_glands,
  model_simple_kidneys
)

models_list_phylo <- list(
  model_phylo_liver,
  model_phylo_heart,
  model_phylo_pituitary_glands,
  model_phylo_spleen,
  model_phylo_stomach,
  model_phylo_lungs,
  model_phylo_adrenal_glands,
  model_phylo_kidneys
)

# Organ names for simple models
names(models_list_simple) <- c(
  "liver","heart","pituitary_glands","spleen",
  "stomach","lungs","adrenal_glands","kidneys"
)

# Organ names for phylogenetic models
names(models_list_phylo) <- c(
  "liver","heart","pituitary_glands","spleen",
  "stomach","lungs","adrenal_glands","kidneys"
)

# Hypothesis string for lambda (phylogenetic signal), based on the vignette
# available here:
# https://cran.r-project.org/web/packages/brms/vignettes/brms_phylogenetics.html

hyp_lambda <- paste0(
  "sd_phylo__Intercept^2 /",
  "(sd_phylo__Intercept^2 + sd_sex__Intercept^2 + sigma^2) = 0"
)

# Make a function to do the process less repetitive. It will extract from each
# model sex-specific coefficients + global slope + lambda (for phylo models only)

extract_all_coefs <- function(fit, organ, model_type) {
  # Sex-specific coefficients (random effects)
  coef_sex <- coef(fit)$sex %>%
    as_tibble(rownames = "sex") %>%
    select(sex, matches("Estimate|Q2.5|Q97.5")) %>%
    transmute(
      model_type,
      organ,
      sex,
      intercept      = Estimate.Intercept,
      intercept_low  = Q2.5.Intercept,
      intercept_high = Q97.5.Intercept,
      slope          = Estimate.log10_body_size,
      slope_low      = Q2.5.log10_body_size,
      slope_high     = Q97.5.log10_body_size
    )
  
  # Global slope (fixed effect)
  global_slope_vec <- fixef(fit)["log10_body_size", ]
  global_slope <- tibble(
    global_slope      = global_slope_vec["Estimate"],
    global_slope_low  = global_slope_vec["Q2.5"],
    global_slope_high = global_slope_vec["Q97.5"]
  )
  
  # Global intercept (fixed effect)
  global_intercept_vec <- fixef(fit)["Intercept", ]
  global_intercept <- tibble(
    global_intercept      = global_intercept_vec["Estimate"],
    global_intercept_low  = global_intercept_vec["Q2.5"],
    global_intercept_high = global_intercept_vec["Q97.5"]
  )
  
  # Lambda (phylogenetic signal) - only for phylo models
  if(model_type == "phylo") {
    hyp_lambda <- hypothesis(fit, hyp_lambda, class = NULL)
    lambda_vals <- hyp_lambda$hypothesis %>%
      transmute(
        lambda         = Estimate
      )
  } else {
    lambda_vals <- tibble(lambda = NA_real_)
  }
  
  # Combine everything
  coef_sex %>%
    bind_cols(global_slope, global_intercept, lambda_vals)
}

# Create complete data frames
df_simple_all_organs <- imap_dfr(
  models_list_simple,
  ~ extract_all_coefs(.x, organ = .y, model_type = "simple")
)

df_phylo_all_organs <- imap_dfr(
  models_list_phylo,
  ~ extract_all_coefs(.x, organ = .y, model_type = "phylo")
)

# Final data frame with everything
df_all_organs <- bind_rows(df_simple_all_organs, df_phylo_all_organs)

# add R²
r2_simple <- imap(models_list_simple, ~ bayes_R2(.x)["R2", "Estimate"])
r2_phylo  <- imap(models_list_phylo,  ~ bayes_R2(.x)["R2", "Estimate"])

df_all_organs <- df_all_organs %>%
  arrange(model_type,organ) %>%
  mutate(R2 = case_when(
    model_type == "simple" ~ r2_simple[organ],
    model_type == "phylo"  ~ r2_phylo[organ]
  ))

df_all_organs_ci <- df_all_organs %>%
  mutate(
    lambda = as.numeric(lambda),
    R2 = as.numeric(R2),
    intercept_ci    = sprintf("%.2f [%.2f, %.2f]", 
                              round(intercept, 2), round(intercept_low, 2), round(intercept_high, 2)),
    slope_ci        = sprintf("%.2f [%.2f, %.2f]", 
                              round(slope, 2), round(slope_low, 2), round(slope_high, 2)),
    global_slope_ci = sprintf("%.2f [%.2f, %.2f]", 
                              round(global_slope, 2), round(global_slope_low, 2), round(global_slope_high, 2)),
    global_intercept_ci = sprintf("%.2f [%.2f, %.2f]", 
                             round(global_intercept, 2), round(global_intercept_low, 2), round(global_intercept_high, 2)),
    lambda = sprintf("%.2f", round(lambda, 2)),
    R2 = sprintf("%.2f", round(R2, 2))
  ) %>%
  mutate(across(where(is.list), as.character)) %>%
  select(model_type, organ, sex, global_intercept_ci, intercept_ci, slope_ci, global_slope_ci, lambda, R2) %>%
  arrange(organ) %>%
  as.data.frame(stringsAsFactors = FALSE)

# Table with filters
df_all_organs_ci %>%
  datatable(
    filter = 'top',
    options = list(
      columnDefs = list(
        list(width = "5%", targets = 0),
        list(width = "15%", targets = 1), 
        list(width = "7%", targets = 2),
        list(width = "17%", targets = 3),
        list(width = "17%", targets = 4),
        list(width = "17%", targets = 5),
        list(width = "7%", targets = 6),
        list(width = "7%", targets = 7)
      ),
      pageLength = 25,
      scrollX = TRUE,
      autoWidth = FALSE
    ),
    class = 'cell-border stripe hover',
    escape = FALSE
  ) %>%
  formatStyle(columns = 0:7, borderRight = "1px solid #ddd")

Figure 3

What I will do here is copy the list of simple and phylogenetic models, their names, and extract the LOO values. Then, I will compare these using the loo_compare function, extract the estimates from this comparison, and plot them.

# Define model lists
models_list_simple <- list(
  model_simple_liver, model_simple_heart, model_simple_pituitary_glands,
  model_simple_spleen, model_simple_stomach, model_simple_lungs,
  model_simple_adrenal_glands, model_simple_kidneys
)

names(models_list_simple) <- c(
  "liver", "heart", "pituitary_glands", "spleen",
  "stomach", "lungs", "adrenal_glands", "kidneys"
)

names(models_list_phylo) <- names(models_list_simple)

# Compare models and extract elpd_diff + se
diff_by_organ <- map_dfr(names(models_list_simple), function(org) {

  # calculate LOO
  loo_phylo  <- loo(models_list_phylo[[org]])
  loo_simple <- loo(models_list_simple[[org]])

  # Compare
  comp <- loo_compare(loo_phylo, loo_simple)
  comp_df <- as.data.frame(comp)

  # Extract row ref
  row <- comp_df[comp_df$elpd_diff != 0, , drop = FALSE]

  # It should be one
  if (nrow(row) != 1) {
    stop("Unexpected loo_compare output for organ: ", org)
  }

  # Outputs
  tibble(
    organ     = gsub("_", " ", org),
    diff_elpd = -row$elpd_diff,   # phylo − simple
    se_diff   =  row$se_diff,
    Vetahri_score   = abs((-row$elpd_diff) / row$se_diff)
  )
}) %>%
  mutate(
    sig_level = if_else(Vetahri_score > 4, "large", "ns")
  ) %>%
  arrange(desc(diff_elpd))


# Figure 3
Figure_3a <- ggplot(diff_by_organ, aes(x = reorder(organ, diff_elpd), y = diff_elpd, 
                                colour = sig_level)) +
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = diff_elpd - se_diff, ymax = diff_elpd + se_diff),
                width = 0.3, linewidth = 0.8) +
    coord_flip() +
scale_colour_manual(
  values = c("ns" = "grey30", "large" = "#009E73"),
  name = "Model improvement",
  labels = c("ns" = expression(z <= 4), "large" = expression(z > 4))
) +
labs(
  x = NULL, 
  y = expression(Delta ~ "ELPD" ~ "(phylogenetic - simple)")
  ) +
  theme_pubr(base_size = 12,border = TRUE) +
  theme(
    panel.grid.minor = element_blank(),
    axis.text.y  = element_text(size = 12),
    axis.title.y = element_text(size = 12),
    legend.position = c(0.80, 0.15),
    legend.margin = margin(3, 3, 3, 3),
    legend.title = element_text(size = 10)
  )


# list of models (only brain)
models_list_simple_brain <- list(
  model_simple_brain_all,
  model_simple_brain_fishes,
  model_simple_brain_birds,
  model_simple_brain_mammals
)

models_list_phylo_brain <- list(
  model_phylo_brain_all,
  model_phylo_brain_fishes,
  model_phylo_brain_birds,
  model_phylo_brain_mammals
)

# Clades names for simple models
names(models_list_simple_brain) <- c(
  "All","Fishes","Birds","Mammals"
)

# Organ names for phylogenetic models
names(models_list_phylo_brain) <- names(models_list_simple_brain)

# Compare models and extract elpd_diff + se
diff_by_clade <- map_dfr(names(models_list_simple_brain), function(clade) {

  # calculate LOO
  loo_phylo_brain  <- loo(models_list_phylo_brain[[clade]])
  loo_simple_brain <- loo(models_list_simple_brain[[clade]])

  # compare
  comp <- loo_compare(loo_phylo_brain, loo_simple_brain)
  comp_df_brain <- as.data.frame(comp)

  # extract row no reference
  row <- comp_df_brain[comp_df_brain$elpd_diff != 0, , drop = FALSE]

  # doubkle check...should be one row
  if (nrow(row) != 1) {
    stop("Unexpected loo_compare output for organ: ", fam)
  }

  # make output
  tibble(
    clade     = gsub("_", " ", clade),
    diff_elpd = -row$elpd_diff,   # phylo − simple
    se_diff   =  row$se_diff,
    Vetahri_score   = abs((-row$elpd_diff) / row$se_diff)
  )
}) %>%
  mutate(
    sig_level = if_else(Vetahri_score > 4, "large", "ns")
  ) %>%
  arrange(desc(diff_elpd))


Figure_3b <- ggplot(diff_by_clade, aes(x = reorder(clade, diff_elpd), y = diff_elpd, 
                                colour = sig_level)) +
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = diff_elpd - se_diff, ymax = diff_elpd + se_diff),
                width = 0.15, linewidth = 0.8) +
  coord_flip() +
scale_colour_manual(
  values = c("ns" = "grey30", "large" = "#009E73"),
  name = "Model improvement",
  labels = c("ns" = expression(z <= 4), "large" = expression(z > 4))
) +
labs(
  x = NULL, 
  y = expression(Delta ~ "ELPD" ~ "(phylogenetic - simple)")
  ) +
  theme_pubr(base_size = 12, border = TRUE) +
  theme(
    panel.grid.minor = element_blank(),
    axis.text.y  = element_text(size = 12),
    axis.title.y = element_text(size = 12),
    legend.position = c(0.80, 0.15),
    legend.margin = margin(3, 3, 3, 3),
    legend.title = element_text(size = 10)
  )



Figure_3a_no_leg <- Figure_3a + theme(legend.position = "none")
Figure_3b_no_leg <- Figure_3b + theme(legend.position = "none")

# Extract form the first plot
legend_3 <- get_legend(
  Figure_3a +
    theme(
      legend.position = "right",
      legend.margin   = margin(3, 3, 3, 3),
      legend.title    = element_text(size = 10)
    )
)

# combine figures and legend
Figure_3 <- plot_grid(
  plot_grid(Figure_3a_no_leg, Figure_3b_no_leg,
            labels = c("A", "B"),
            ncol   = 1,
            align  = "hv"),
  legend_3,
  ncol        = 2,          
  rel_widths  = c(2, 1)  
)

Figure_3

# Export figure
ggsave("../outputs/Figure_3.pdf", Figure_3, width = 6, height = 6)

Figure 4

This figure comprises two plots. One illustrates the relationships between slopes, while the other shows intercepts for males and females, both derived from phylogenetic models. For clarity, the manuscript discusses results and estimates exclusively from the phylogenetic models; however, summary outputs for both model types are provided in Tables  S2 (Supplementary Material).

# Define common theme (used by all plots)
common_theme <- theme_bw() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    legend.position = c(0.15, 0.85),
    legend.key.size = unit(0.4, "cm"),
    legend.text = element_text(size = 8),
    legend.title = element_blank(),
    legend.spacing.y = unit(0.5, "cm"),
    legend.margin = margin(2, 2, 2, 2),
    legend.background = element_rect(fill = "white", size = 0.3),
    plot.margin = margin(15, 15, 15, 15),
    axis.text = element_text(size = 12),
    axis.title = element_text(size = 12),
    axis.title.x = element_text(margin = margin(t = 12)),
    axis.title.y = element_text(margin = margin(r = 5)),
    panel.border = element_rect(colour = "black", fill = NA, size = 1)
  )

# Define legend theme
no_legend_theme <- theme_bw() +
  theme(
    panel.grid = element_blank(),
    legend.position = "left",
    plot.margin = margin(15, 15, 15, 15),
    axis.text = element_text(size = 12),
    axis.title = element_text(size = 12),
    axis.title.x = element_text(margin = margin(t = 12)),
    axis.title.y = element_text(margin = margin(r = 5)),
    panel.border = element_rect(colour = "black", fill = NA, size = 1)
  )

# Clean organ names and colours
df_all_organs$organ <- gsub("_", " ", df_all_organs$organ)
all_organs <- sort(unique(df_all_organs$organ))
my_cols <- setNames(c("#4169E1", "#DC143C", "#D1600F", "#8B0000", 
                      "#FFFF00", "#8DA715", "#FF69B4", "black"),
                    all_organs)

# Plot 1: slopes - phylogenetic model (with legend)
df1 <- df_all_organs %>% 
  filter(model_type == "phylo", sex %in% c("male", "female")) %>%
  select(organ, sex, slope, slope_low, slope_high) %>%
  pivot_wider(names_from = sex, values_from = c(slope, slope_low, slope_high), names_sep = "_") %>%
  mutate(male_xmin = slope_low_male, male_xmax = slope_high_male,
         female_ymin = slope_low_female, female_ymax = slope_high_female)

p1 <- ggplot(df1, aes(x = slope_male, y = slope_female)) +
  geom_errorbarh(aes(xmin = male_xmin, xmax = male_xmax, colour = organ), height = 0.038, size = 0.8) +
  geom_errorbar(aes(ymin = female_ymin, ymax = female_ymax, colour = organ), width = 0.038, size = 0.8) +
  geom_point(aes(fill = organ), size = 3, shape = 21, stroke = 0.6, colour = "white") +
  geom_abline(intercept = 0, slope = 1, linetype = "dashed", colour = "grey40", size = 0.6, alpha = 0.5) +
  scale_colour_manual(values = my_cols) +
  scale_fill_manual(values = my_cols, guide = "none") +
  guides(colour = guide_legend(
    override.aes = list(shape = 21, colour = "white", stroke = 0.6, size = 3))) +
  labs(x = "Male slope", y = "Female slope") +
  coord_fixed() + common_theme +
  scale_x_continuous(limits = c(0.3, 1.5), breaks = seq(0.3, 1.5, by = 0.3)) +
  scale_y_continuous(limits = c(0.3, 1.5), breaks = seq(0.3, 1.5, by = 0.3))


# Plot 2: intercepts - phylogenetic model
df2 <- df_all_organs %>% 
  filter(model_type == "phylo", sex %in% c("male", "female")) %>%
  select(organ, sex, intercept, intercept_low, intercept_high) %>%
  pivot_wider(names_from = sex, values_from = c(intercept, intercept_low, intercept_high), names_sep = "_") %>%
  mutate(male_xmin = intercept_low_male, male_xmax = intercept_high_male,
         female_ymin = intercept_low_female, female_ymax = intercept_high_female)

p2 <- ggplot(df2, aes(x = intercept_male, y = intercept_female)) +
  geom_errorbarh(aes(xmin = male_xmin, xmax = male_xmax, colour = organ), height = 0.12, size = 0.8) +
  geom_errorbar(aes(ymin = female_ymin, ymax = female_ymax, colour = organ), width = 0.12, size = 0.8) +
  geom_point(aes(fill = organ), size = 3, shape = 21, stroke = 0.6, colour = "white") +
  geom_abline(intercept = 0, slope = 1, linetype = "dashed", colour = "grey40", size = 0.6, alpha = 0.5) +
  scale_colour_manual(values = my_cols, guide = "none") +
  scale_fill_manual(values = my_cols, guide = "none") +
  labs(x = "Male intercept", y = "Female intercept") +
  coord_fixed() + no_legend_theme +
  scale_x_continuous(limits = c(-5, 2), breaks = seq(-5, 2, by = 1)) +
  scale_y_continuous(limits = c(-5, 2), breaks = seq(-5, 2, by = 1))

# temporal plot, just for the legend
legend_plot <- ggplot(df1, aes(x = 1, y = 1, fill = organ)) +
  geom_point(shape = 21, size = 4, stroke = 0.8, colour = "white") +
  scale_fill_manual(values = my_cols) +
  theme_void() +
  guides(fill = guide_legend(title = NULL, override.aes = list(shape = 21, colour = "white", stroke = 0.6, size = 3))) +
  theme(legend.position = "right", legend.spacing.y = unit(0.3, "cm"))

# Extract legend from plot from the previous plot
legend_right <- get_legend(legend_plot+ theme(legend.position = "right", legend.spacing.y = unit(0.3, "cm")))

# Add model labels to panels
p1 <- p1 + theme(legend.position = "none", plot.margin = margin(2, 2, 2, 2))

p2 <- p2 + theme(legend.position = "none", plot.margin = margin(2, 2, 2, 2)) 

# Combine panels and create final figure
combined_plot <- (p1 / p2) + plot_layout(widths = c(1, 1), heights = c(1, 1))

Figure_4 <- plot_grid(combined_plot, legend_right, rel_widths = c(2, 1))

Figure_4

# Export figure
ggsave("../outputs/Figure_4.pdf", Figure_4, width = 5, height = 7)

Table S3

## Extract coefficients from brain models 

# Hypothesis string for lambda (phylogenetic signal), based on the vignette
# available here:
# https://cran.r-project.org/web/packages/brms/vignettes/brms_phylogenetics.html

hyp_lambda <- paste0(
  "sd_phylo__Intercept^2 /",
  "(sd_phylo__Intercept^2 + sd_sex__Intercept^2 + sigma^2) = 0"
)

# To do the process less repetitive I extract from each
# model sex-specific coefficients + global slope + lambda (for phylo models only)
extract_all_coefs <- function(fit, clade, model_type) {
  # Sex-specific coefficients (random effects)
  coef_sex <- coef(fit)$sex %>%
    as_tibble(rownames = "sex") %>%
    select(sex, matches("Estimate|Q2.5|Q97.5")) %>%
    transmute(
      model_type,
      clade,
      sex,
      intercept      = Estimate.Intercept,
      intercept_low  = Q2.5.Intercept,
      intercept_high = Q97.5.Intercept,
      slope          = Estimate.log10_body_size,
      slope_low      = Q2.5.log10_body_size,
      slope_high     = Q97.5.log10_body_size
    )
  
  # Global slope (fixed effect)
  global_slope_vec <- fixef(fit)["log10_body_size", ]
  global_slope <- tibble(
    global_slope      = global_slope_vec["Estimate"],
    global_slope_low  = global_slope_vec["Q2.5"],
    global_slope_high = global_slope_vec["Q97.5"]
  )
  
    # Global intercept (fixed effect)
  global_intercept_vec <- fixef(fit)["Intercept", ]
  global_intercept <- tibble(
    global_intercept      = global_intercept_vec["Estimate"],
    global_intercept_low  = global_intercept_vec["Q2.5"],
    global_intercept_high = global_intercept_vec["Q97.5"]
  )
  
  # Lambda (phylogenetic signal) - only for phylo models
  if(model_type == "phylo") {
    hyp_lambda <- hypothesis(fit, hyp_lambda, class = NULL)
    lambda_vals <- hyp_lambda$hypothesis %>%
      transmute(
        lambda         = Estimate
      )
  } else {
    lambda_vals <- tibble(lambda = NA_real_)
  }
  
  # Combine everything
  coef_sex %>%
    bind_cols(global_slope, global_intercept, lambda_vals)
}

# Create complete data frames
df_simple_all_clade <- imap_dfr(
  models_list_simple_brain,
  ~ extract_all_coefs(.x, clade = .y, model_type = "simple")
)

df_phylo_all_clade<- imap_dfr(
  models_list_phylo_brain,
  ~ extract_all_coefs(.x, clade = .y, model_type = "phylo")
)

# Final data frame with everything
df_all_clade <- bind_rows(df_simple_all_clade, df_phylo_all_clade)

# add R²
r2_simple <- imap(models_list_simple_brain, ~ bayes_R2(.x)["R2", "Estimate"])
r2_phylo  <- imap(models_list_phylo_brain,  ~ bayes_R2(.x)["R2", "Estimate"])

df_all_clade <- df_all_clade %>%
  arrange(model_type,clade) %>%
  mutate(R2 = case_when(
    model_type == "simple" ~ r2_simple[clade],
    model_type == "phylo"  ~ r2_phylo[clade]
  ))

df_all_clade <- df_all_clade %>%
  mutate(
    lambda = as.numeric(lambda),
    R2 = as.numeric(R2),
    intercept_ci    = sprintf("%.2f [%.2f, %.2f]", 
                              round(intercept, 2), round(intercept_low, 2), round(intercept_high, 2)),
    slope_ci        = sprintf("%.2f [%.2f, %.2f]", 
                              round(slope, 2), round(slope_low, 2), round(slope_high, 2)),
    global_slope_ci = sprintf("%.2f [%.2f, %.2f]", 
                              round(global_slope, 2), round(global_slope_low, 2), round(global_slope_high, 2)),
    global_intercept_ci = sprintf("%.2f [%.2f, %.2f]", 
                              round(global_intercept, 2), round(global_intercept_low, 2), round(global_intercept_high, 2)),
    lambda = sprintf("%.2f", round(lambda, 2)),
    R2 = sprintf("%.2f", round(R2, 2))
  ) %>%
  mutate(across(where(is.list), as.character)) %>%
  select(model_type, clade, sex, global_intercept_ci, intercept_ci, slope_ci, global_slope_ci, lambda, R2) %>%
  arrange(clade) %>%
  as.data.frame(stringsAsFactors = FALSE)

# Table with filters remains identical
df_all_clade %>%
  datatable(
    filter = 'top',
    options = list(
      columnDefs = list(
        list(width = "5%", targets = 0),
        list(width = "15%", targets = 1), 
        list(width = "7%", targets = 2),
        list(width = "17%", targets = 3),
        list(width = "17%", targets = 4),
        list(width = "17%", targets = 5),
        list(width = "7%", targets = 6),
        list(width = "7%", targets = 7)
      ),
      pageLength = 25,
      scrollX = TRUE,
      autoWidth = FALSE
    ),
    class = 'cell-border stripe hover',
    escape = FALSE
  ) %>%
  formatStyle(columns = 0:7, borderRight = "1px solid #ddd")

Figure 5

clade_names <- c("All", "Fishes", "Birds", "Mammals")

sex_colours <- c("male" = "#008B8B", "female" = "#9370DB", "Grand mean" = "gray70")

common_theme <- theme_bw(base_size = 12) +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    legend.position = "none",
    panel.spacing.x = unit(0.3, "lines"),
    axis.text.y = element_blank(),
    axis.text.x = element_text(size = 12),
    axis.title.x = element_text(margin = margin(t = 8)),
    axis.title.y = element_text(margin = margin(r = 8)),
    axis.title = element_text(size = 14),
    plot.margin = margin(t = 4, r = 8, b = 4, l = 8, unit = "pt")
  )

plot_slopes <- function(model, clade_name) {
  fixed_effects <- model %>%
    spread_draws(b_log10_body_size)

  random_slopes <- model %>%
    gather_draws(r_sex[sex, term], regex = TRUE) %>%
    filter(term == "log10_body_size")

  slope_draws <- random_slopes %>%
    left_join(
      fixed_effects,
      by = c(".chain", ".iteration", ".draw")
    ) %>%
    mutate(
      slope = .value + b_log10_body_size,
      sex_label = glue::glue("{sex}")
    ) %>%
    select(.draw, sex_label, slope)

  global_slope_draws <- fixed_effects %>%
    mutate(
      slope = b_log10_body_size,
      sex_label = "Grand mean"
    ) %>%
    select(.draw, sex_label, slope)

  slope_draws_combined <- bind_rows(slope_draws, global_slope_draws)

  slope_draws_combined$sex_label <- factor(
    slope_draws_combined$sex_label,
    levels = c(sort(unique(slope_draws$sex_label)), "Grand mean")
  )

  ggplot(
    slope_draws_combined,
    aes(
      x = slope,
      y = sex_label,
      fill = sex_label
    )
  ) +
    annotate(
      "rect",
      xmin = -0.6, xmax = 1.6,
      ymin = 3.69, ymax = 3.99,
      fill = "grey90", color = NA
    ) +
    annotate(
      "text",
      x = 0.5, y = 3.85,
      label = clade_name,
      size = 4.5, fontface = "bold",
      hjust = 0.5, vjust = 0.5
    ) +
    stat_halfeye(
      scale = 0.8,
      justification = 0.08,
      slab_alpha    = 0.8,  
      interval_alpha = 1,   
      point_alpha    = 1,   
      colour         = "black" 
    ) +
    scale_y_discrete(expand = expansion(mult = c(0.02, 0.35))) +
    coord_cartesian(xlim = c(-0.5, 1.5), clip = "off") +
    scale_fill_manual(values = sex_colours) +
    labs(x = "", y = NULL, title = NULL) +
    common_theme
}

plot_intercepts <- function(model, clade_name) {
  fixed_effects <- model %>%
    spread_draws(b_Intercept)

  random_intercepts <- model %>%
    gather_draws(r_sex[sex, term], regex = TRUE) %>%
    filter(term == "Intercept")

  intercept_draws <- random_intercepts %>%
    left_join(
      fixed_effects,
      by = c(".chain", ".iteration", ".draw")
    ) %>%
    mutate(
      intercept = .value + b_Intercept,
      sex_label = glue::glue("{sex}")
    ) %>%
    select(.draw, sex_label, intercept)

  global_intercept_draws <- fixed_effects %>%
    mutate(
      intercept = b_Intercept,
      sex_label = "Grand mean"
    ) %>%
    select(.draw, sex_label, intercept)

  intercept_draws_combined <- bind_rows(intercept_draws, global_intercept_draws)

  intercept_draws_combined$sex_label <- factor(
    intercept_draws_combined$sex_label,
    levels = c(sort(unique(intercept_draws$sex_label)), "Grand mean")
  )

  ggplot(
    intercept_draws_combined,
    aes(
      x = intercept,
      y = sex_label,
      fill = sex_label
    )
  ) +
    annotate(
      "rect",
      xmin = -Inf, xmax = Inf,
      ymin = 3.69, ymax = 3.99,
      fill = "grey90", color = NA
    ) +
    annotate(
      "text",
      x = 0, y = 3.85,
      label = clade_name,
      size = 4.5, fontface = "bold",
      hjust = 0.5, vjust = 0.5
    ) +
    stat_halfeye(
      scale = 0.8,
      justification = 0.08,
      slab_alpha     = 0.8,
      interval_alpha = 1,
      point_alpha    = 1,
      colour         = "black"
    ) +
    scale_x_continuous(breaks = seq(-3, 3, 1.5)) +
    scale_y_discrete(expand = expansion(mult = c(0.02, 0.35))) +
    coord_cartesian(xlim = c(-3.2, 3.2), clip = "off") +
    scale_fill_manual(values = sex_colours) +
    labs(x = "", y = NULL) +
    common_theme
}

# build lists of plots (assuming models_list_phylo_brain is your list of 4 models)
intercept_plots <- lapply(seq_along(models_list_phylo_brain), function(i) {
  plot_intercepts(models_list_phylo_brain[[i]], clade_names[i])
})

slope_plots <- lapply(seq_along(models_list_phylo_brain), function(i) {
  plot_slopes(models_list_phylo_brain[[i]], clade_names[i])
})

# combined figure
Figure_5 <- cowplot::plot_grid(
  plotlist = c(
    intercept_plots[[1]], slope_plots[[1]],
    intercept_plots[[2]], slope_plots[[2]],
    intercept_plots[[3]], slope_plots[[3]],
    intercept_plots[[4]], slope_plots[[4]]
  ),
  ncol = 2,
  nrow = 4,
  align = "hv"
)

Figure_5

ggsave("../outputs/Figure_5.pdf", Figure_5, width = 10, height = 14)

Session information

R version 4.5.2 (2025-10-31)

Platform: aarch64-apple-darwin20

attached base packages: stats, graphics, grDevices, utils, datasets, methods and base

other attached packages: ggtreeExtra(v.1.21.0), ggnewscale(v.0.5.2), DT(v.0.34.0), ggridges(v.0.5.7), RColorBrewer(v.1.1-3), modelr(v.0.1.11), broom(v.1.0.12), viridis(v.0.6.5), viridisLite(v.0.4.3), posterior(v.1.6.1), lubridate(v.1.9.5), forcats(v.1.0.1), readr(v.2.2.0), tidyverse(v.2.0.0), furrr(v.0.3.1), future(v.1.69.0), cmdstanr(v.0.9.0), purrr(v.1.2.1), loo(v.2.9.0), bayestestR(v.0.17.0), tidybayes(v.3.0.7), bayesplot(v.1.15.0), dplyr(v.1.2.0), brms(v.2.23.0), Rcpp(v.1.1.1), tidyr(v.1.3.2), stringr(v.1.6.0), details(v.0.4.0), ggthemes(v.5.2.0), tibble(v.3.3.1), ggtree(v.4.0.5), patchwork(v.1.3.2), cowplot(v.1.2.0), ggpubr(v.0.6.3), ggplot2(v.4.0.2), kableExtra(v.1.4.0), readxl(v.1.4.5), phytools(v.2.5-2), maps(v.3.4.3), ape(v.5.8-1) and knitr(v.1.51)

loaded via a namespace (and not attached): svUnit(v.1.0.8), splines(v.4.5.2), ggplotify(v.0.1.3), cellranger(v.1.1.0), lifecycle(v.1.0.5), rstatix(v.0.7.3), StanHeaders(v.2.32.10), doParallel(v.1.0.17), globals(v.0.19.0), processx(v.3.8.6), lattice(v.0.22-9), MASS(v.7.3-65), insight(v.1.4.6), crosstalk(v.1.2.2), ggdist(v.3.3.3), backports(v.1.5.0), magrittr(v.2.0.4), sass(v.0.4.10), rmarkdown(v.2.31), jquerylib(v.0.1.4), yaml(v.2.3.12), otel(v.0.2.0), pkgbuild(v.1.4.8), multcomp(v.1.4-29), abind(v.1.4-8), expm(v.1.0-0), quadprog(v.1.5-8), yulab.utils(v.0.2.4), TH.data(v.1.1-5), tensorA(v.0.36.2.1), rappdirs(v.0.3.4), sandwich(v.3.1-1), gdtools(v.0.5.0), inline(v.0.3.21), listenv(v.0.10.0), tidytree(v.0.4.7), bridgesampling(v.1.2-1), parallelly(v.1.46.1), svglite(v.2.2.2), codetools(v.0.2-20), xml2(v.1.5.2), tidyselect(v.1.2.1), aplot(v.0.2.9), farver(v.2.1.2), stats4(v.4.5.2), matrixStats(v.1.5.0), jsonlite(v.2.0.0), Formula(v.1.2-5), survival(v.3.8-6), iterators(v.1.0.14), emmeans(v.2.0.2), systemfonts(v.1.3.2), foreach(v.1.5.2), tools(v.4.5.2), ragg(v.1.5.1), treeio(v.1.34.0), glue(v.1.8.0), mnormt(v.2.1.2), gridExtra(v.2.3), xfun(v.0.57), distributional(v.0.7.0), withr(v.3.0.2), numDeriv(v.2016.8-1.1), combinat(v.0.0-8), fastmap(v.1.2.0), digest(v.0.6.39), timechange(v.0.4.0), R6(v.2.6.1), gridGraphics(v.0.5-1), estimability(v.1.5.1), colorspace(v.2.1-2), textshaping(v.1.0.5), utf8(v.1.2.6), generics(v.0.1.4), data.table(v.1.18.2.1), fontLiberation(v.0.1.0), clusterGeneration(v.1.3.8), httr(v.1.4.8), htmlwidgets(v.1.6.4), scatterplot3d(v.0.3-45), pkgconfig(v.2.0.3), gtable(v.0.3.6), S7(v.0.2.1), htmltools(v.0.5.9), fontBitstreamVera(v.0.1.1), carData(v.3.0-6), scales(v.1.4.0), png(v.0.1-8), ggfun(v.0.2.0), rstudioapi(v.0.18.0), reshape2(v.1.4.5), tzdb(v.0.5.0), curl(v.7.0.0), coda(v.0.19-4.1), checkmate(v.2.3.4), nlme(v.3.1-168), cachem(v.1.1.0), zoo(v.1.8-15), DEoptim(v.2.2-8), parallel(v.4.5.2), desc(v.1.4.3), pillar(v.1.11.1), grid(v.4.5.2), vctrs(v.0.7.2), car(v.3.1-5), arrayhelpers(v.1.1-0), xtable(v.1.8-8), evaluate(v.1.0.5), mvtnorm(v.1.3-6), cli(v.3.6.5), compiler(v.4.5.2), rlang(v.1.1.7), rstantools(v.2.6.0), ggsignif(v.0.6.4), labeling(v.0.4.3), ps(v.1.9.1), plyr(v.1.8.9), fs(v.2.0.1), pander(v.0.6.6), rstan(v.2.32.7), ggiraph(v.0.9.6), stringi(v.1.8.7), QuickJSR(v.1.9.0), lazyeval(v.0.2.2), optimParallel(v.1.0-2), V8(v.8.0.1), Brobdingnag(v.1.2-9), fontquiver(v.0.2.1), Matrix(v.1.7-4), hms(v.1.1.4), clipr(v.0.8.0), igraph(v.2.2.2), bslib(v.0.10.0), RcppParallel(v.5.1.11-2), phangorn(v.2.12.1) and fastmatch(v.1.1-8)