Owl duetting analysis

Owl duetting evolution

Author

Marcelo Araya Salas

Published

June 24, 2024

Source code and data found at

1 Purpose

Assess factors shaping the evolution of duettting in owls

2 Load packages

Code

# knitr is require for creating html/pdf/word reports
# formatR is used for soft-wrapping code


if (!"sketchy" %in% installed.packages()[,1])
    remotes::install_github("maRce10/sketchy")

# unload R packages (specific packages required for this code)
x <-
    c(
        "ggplot2",
        "kableExtra",
        "knitr",
        "formatR",
        "rprojroot",
        "xaringanExtra",
        "viridis",
        "brms",
        "ape",
        "brmsish",
        "corrplot",
        "phytools",
        "rstan",
        github = "maRce10/brmsish"
    )

# install/ load packages
sketchy::load_packages(packages = x)

## set function to prune trees
prep_tree <- function(tree, spp) {
    spp <- unique(spp)
    if (is(tree, "multiPhylo"))
        trees <- lapply(tree, function(x)
            drop.tip(phy = x, tip = setdiff(x$tip.label, spp)))
    else
        trees <-
            drop.tip(phy = tree, tip = setdiff(tree$tip.label, spp))
    
    return(tree)
}

rstan_options(auto_write = TRUE)

Code

## Open dataframes

# read data
trees_owl <- read.nexus("./data/raw/BirdTree_Owl_Output.nex")
consensus_tree_owl <- read.nexus("./data/raw/Owl_Consensus_Tree.nex")

set.seed(123)
trees_owl <- trees_owl[sample(1:100, size = 40)]
corr_df <- df <- readRDS("./data/raw/Owl_Duet.RDS")

# Set categorical response variable as ordered factors (ordinal data) and  #USE FOR BRMS
df$tip.labels <-
    as.factor(df$tip.labels) # binary categorical response variable
df$Duetting <-
    as.factor(df$Duetting) # categorical explanatory variable
df$Social_Bond <-
    as.factor(df$Social_Bond) # categorical explanatory variable
df$Territoriality <-
    as.factor(df$Territoriality) # categorical explanatory variable
df$White_Plumage <-
    as.factor(df$White_Plumage) # categorical explanatory variable
df$Nocturnality <-
    as.factor(df$Nocturnality) # categorical explanatory variable
df$Habitat_Density <-
    as.factor(df$Habitat_Density) # categorical explanatory variable
df$Migration <-
    as.factor(df$Migration) # categorical explanatory variable

# scale and center numerical data
df$Central_Latitude <-
    scale(df$Central_Latitude, center = TRUE, scale = TRUE)

# use if you need to reset rownames as sequential numbers
rownames(df) <- 1:nrow(df)

# look at structure
str(df)

'data.frame':   63 obs. of  9 variables:
 $ tip.labels      : Factor w/ 63 levels "Aegolius_acadicus",..: 1 2 3 4 5 6 7 8 9 10 ...
 $ Duetting        : Factor w/ 2 levels "0","1": 2 1 1 2 2 2 1 1 2 2 ...
 $ Social_Bond     : Factor w/ 3 levels "0","1","2": 1 1 1 1 2 1 1 1 2 3 ...
 $ Territoriality  : Factor w/ 2 levels "0","1": 1 1 1 1 2 1 1 1 2 2 ...
 $ White_Plumage   : Factor w/ 2 levels "0","1": 2 1 1 1 1 2 2 2 2 2 ...
 $ Nocturnality    : Factor w/ 2 levels "0","1": 2 2 1 2 2 1 2 1 1 2 ...
 $ Habitat_Density : Factor w/ 3 levels "0","1","2": 2 1 3 2 1 2 2 3 3 2 ...
 $ Migration       : Factor w/ 3 levels "0","1","2": 3 2 3 2 1 1 1 2 1 1 ...
 $ Central_Latitude: num [1:63, 1] 1.06 1.42 1.16 1.2 -1.04 ...
  ..- attr(*, "scaled:center")= num 15.9
  ..- attr(*, "scaled:scale")= num 29.1

Code

# View(df)

3 DATA EXPLORATION AND PRE-ANALYSIS

Code

######################
# Response variables #
######################

# Make the correlation table
### NOTE: DEFAULT CODE DOES NOT REMOVE MISSING DATA; use complete.obs to drop nas.
corrtab <-
    cor(corr_df[,-1], method = "pearson", use = "complete.obs")
round(corrtab, 2)

                 Duetting Social_Bond Territoriality White_Plumage Nocturnality
Duetting             1.00        0.16           0.27         -0.08         0.31
Social_Bond          0.16        1.00           0.41         -0.01         0.20
Territoriality       0.27        0.41           1.00          0.03         0.11
White_Plumage       -0.08       -0.01           0.03          1.00        -0.39
Nocturnality         0.31        0.20           0.11         -0.39         1.00
Habitat_Density     -0.18        0.04          -0.19          0.27        -0.39
Migration           -0.28       -0.45          -0.47         -0.01        -0.21
Central_Latitude    -0.03       -0.13          -0.32          0.09        -0.16
                 Habitat_Density Migration Central_Latitude
Duetting                   -0.18     -0.28            -0.03
Social_Bond                 0.04     -0.45            -0.13
Territoriality             -0.19     -0.47            -0.32
White_Plumage               0.27     -0.01             0.09
Nocturnality               -0.39     -0.21            -0.16
Habitat_Density             1.00      0.25             0.24
Migration                   0.25      1.00             0.36
Central_Latitude            0.24      0.36             1.00

Code

# Investigate visually correlations using library(ellipse)
corrplot(
    corrtab,
    mar = c(0.1, 0.1, 0.1, 0.1),
    method = "number",
    tl.cex = 0.25,
    tl.offset = 2
)

Code

corrplot(
    corrtab,
    mar = c(0.1, 0.1, 0.1, 0.1),
    tl.cex = 0.25,
    tl.offset = 2
)

Code

# Set categorical data to analyze as numbers
df$duet <- as.numeric(df$Duetting)


### phylogenetic signal for response variable
trait_duet <- df$duet
length(trait_duet)

[1] 63

Code

names(trait_duet) <- df$tip.labels


tree_duet <- phangorn::mcc(trees_owl)
phylosig(
    tree_duet,
    trait_duet,
    method = "K",
    test = TRUE,
    nsim = 1000
)


Phylogenetic signal K : 0.386789 
P-value (based on 1000 randomizations) : 0.023

Code

phylosig(
    tree_duet,
    trait_duet,
    method = "lambda",
    test = TRUE,
    nsim = 1000
)


Phylogenetic signal lambda : 0.251294 
logL(lambda) : -37.722 
LR(lambda=0) : 1.00417 
P-value (based on LR test) : 0.316303

4 Statistical analysis

5 Taking into account phylogenetic uncertainty

5.1 Full model

Code

### See https://urldefense.proofpoint.com/v2/url?u=https-3A__psyarxiv.com_x8swp_&d=DwIGAg&c=t1KneP5CA95qj21Jnj2QJw&r=uXrXEa22P1XTtZOkZJzsx8pAQCDanQ2sMUk_AfMLQ0I&m=Th_rFO1G4mV3f039Te1VuvpmcjwF2KXjHAfcU0UhAb8vpbbU-nEL3HwCTlICHqS3&s=qolX0qBhJIjMDpanzBhiCR3QCeDrzGH4A9eyiqlhW2g&e=  for explanation of types of ordinal models implemented in brms.
# And https://urldefense.proofpoint.com/v2/url?u=https-3A__solomonkurz.netlify.app_post_2021-2D12-2D29-2Dnotes-2Don-2Dthe-2Dbayesian-2Dcumulative-2Dprobit_&d=DwIGAg&c=t1KneP5CA95qj21Jnj2QJw&r=uXrXEa22P1XTtZOkZJzsx8pAQCDanQ2sMUk_AfMLQ0I&m=Th_rFO1G4mV3f039Te1VuvpmcjwF2KXjHAfcU0UhAb8vpbbU-nEL3HwCTlICHqS3&s=PQ-Sgi8kE3P2Yamj_rtG-r2rZUbEB5iytYjn9NLil08&e= 

# set trees
#A <- vcv.phylo(trees)

# set parameters
set.seed(123)
parallel <-
    7 # set this to the number of cores you want to use in the cluster
thin <- 100
priors <- c(prior(normal(0, 2), class = "b"))
chains <- 4


## run model

### Owl duets: present/absent vs explanatory variables
mod <-
    phylogenetic_uncertainty(
        formula = Duetting ~ Social_Bond + Territoriality + Nocturnality + Habitat_Density + Migration + Central_Latitude + White_Plumage,
        data = df,
        phylos = prep_tree(trees_owl, spp = df$tip.labels),
        family = bernoulli(link = "logit"),
        chains = chains,
        cores = chains,
        iter = 20000,
        prior = priors,
        control = list(adapt_delta = 0.99, max_treedepth = 15),
        init = 0,
        thin = thin,
        sp.id.column = "tip.labels"
    )

saveRDS(mod, "./data/processed/phylogenetic_uncertainty_model.RDS")

Code

extended_summary(read.file = "./data/processed/phylogenetic_uncertainty_model.RDS", plot.area.prop = 0.5)

Warning: Removed 170 rows containing missing values or values outside the scale range
(`stat_slabinterval()`).

5.2 phylogenetic_uncertainty_model

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 2) Intercept-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5)	Duetting ~ Social_Bond + Territoriality + Nocturnality + Habitat_Density + Migration + Central_Latitude + White_Plumage + (1 \| gr(tip.labels, cov = vcv.phylo))	20000	400	1	10000	225 (5.625e-05%)	0	26240.5	79221.7	72886991

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	-1.287	-8.943	6.704	1.002	118134.14	837992.68
b_Social_Bond1	0.910	-2.619	4.414	1.006	34212.31	135590.81
b_Social_Bond2	-0.327	-3.784	3.150	1.006	39120.93	129019.67
b_Territoriality1	1.579	-1.976	4.959	1.008	26865.05	85685.49
b_Nocturnality1	0.820	-2.894	4.453	1.001	278630.75	2815847.06
b_Habitat_Density1	0.550	-2.838	3.869	1.006	39403.44	154156.28
b_Habitat_Density2	-0.797	-4.465	2.904	1.001	328170.98	950503.07
b_Migration1	-0.893	-4.393	2.732	1.002	156732.86	2362920.23
b_Migration2	-1.089	-4.750	2.681	1.003	93845.56	489123.09
b_Central_Latitude	-0.521	-3.623	2.230	1.008	26240.50	79221.70
b_White_Plumage1	-0.127	-3.637	3.421	1.004	55221.58	445307.11

5.3 Specific models

Code

models <-
    c(
        "Duetting ~Social_Bond + Territoriality + Nocturnality + Habitat_Density + Migration + Central_Latitude",
        "White_Plumage ~ Social_Bond + Territoriality + Nocturnality + Habitat_Density + Migration + Central_Latitude",
        "Duetting ~ White_Plumage",
        "White_Plumage ~ Duetting",
        "Duetting ~Social_Bond",
        "Duetting ~ Territoriality",
        "Duetting ~ Nocturnality",
        "Duetting ~ Habitat_Density",
        "Duetting ~ Migration",
        "Duetting ~ Central_Latitude",
        "White_Plumage ~Social_Bond",
        "White_Plumage ~ Territoriality",
        "White_Plumage ~ Nocturnality",
        "White_Plumage ~ Habitat_Density",
        "White_Plumage ~ Migration",
        "White_Plumage ~ Central_Latitude"
    )

Code

# set parameters
set.seed(123)
parallel <- 1
thin <- 100
priors <- c(prior(normal(0, 2), class = "b"))
chains <- 2


## run model

### Owl duets: present/absent vs explanatory variables

for (i in models) {
    
    print(i)
    mod_name <- gsub("~", "_by_", i)
    mod_name <- gsub(" \\+ ", "_", mod_name)
    mod_name <- paste0(mod_name, ".RDS")
    mod_name <- file.path("./data/processed", mod_name)

    if (!file.exists(mod_name)){
    
    mod <-
        phylogenetic_uncertainty(
            formula = as.formula(i),
            data = df,
            phylos = prep_tree(trees_owl, spp = df$tip.labels),
            family = bernoulli(link = "logit"),
            chains = chains,
            cores = 20,
            iter = 10000,
            prior = priors,
            control = list(adapt_delta = 0.99, max_treedepth = 15),
            init = 0,
            thin = thin,
            sp.id.column = "tip.labels"
        )
    

    saveRDS(mod, mod_name)
    rm("mod")
    }
}

Code

# mod <- readRDS("./data/processed/phylogenetic_uncertainty_model.RDS")

for (i in models) {

     mod_name <- gsub("~", "_by_", i)
    mod_name <- gsub(" \\+ ", "_", mod_name)
    mod_name <- paste0(mod_name, ".RDS")
    mod_name <- file.path("./data/processed", mod_name)
    
 extended_summary(read.file = mod_name, plot.area.prop = 0.5, fit.name = i, highlight = TRUE)
}

Warning: Removed 51 rows containing missing values or values outside the scale range
(`stat_slabinterval()`).

5.4 Duetting ~Social_Bond + Territoriality + Nocturnality + Habitat_Density + Migration + Central_Latitude

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 2) Intercept-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5)	Duetting ~ Social_Bond + Territoriality + Nocturnality + Habitat_Density + Migration + Central_Latitude + (1 \| gr(tip.labels, cov = vcv.phylo))	10000	400	1	5000	18 (9e-06%)	0	27375.21	88754.48	455137784

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	-1.332	-8.747	6.411	1.002	171036.77	905924.50
b_Social_Bond1	0.900	-2.581	4.369	1.006	36116.16	128389.51
b_Social_Bond2	-0.323	-3.750	3.122	1.006	39757.33	150010.59
b_Territoriality1	1.581	-1.918	4.918	1.007	29952.71	88754.48
b_Nocturnality1	0.845	-2.840	4.445	1.001	325691.63	1392603.62
b_Habitat_Density1	0.544	-2.786	3.818	1.005	42682.24	190722.96
b_Habitat_Density2	-0.807	-4.446	2.862	1.001	710815.58	1510430.18
b_Migration1	-0.908	-4.368	2.671	1.002	188455.41	1381627.57
b_Migration2	-1.122	-4.750	2.625	1.003	102268.76	1221411.01
b_Central_Latitude	-0.477	-3.545	2.190	1.008	27375.21	92001.37

Warning: Removed 69 rows containing missing values or values outside the scale range
(`stat_slabinterval()`).

5.5 White_Plumage ~ Social_Bond + Territoriality + Nocturnality + Habitat_Density + Migration + Central_Latitude

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 2) Intercept-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5)	White_Plumage ~ Social_Bond + Territoriality + Nocturnality + Habitat_Density + Migration + Central_Latitude + (1 \| gr(tip.labels, cov = vcv.phylo))	10000	200	1	5000	28 (2.8e-05%)	0	6325.64	45201.65	1579423574

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	-1.233	-10.426	6.895	1.002	64771.29	449506.85
b_Social_Bond1	0.085	-3.593	3.768	1.007	16391.97	154004.83
b_Social_Bond2	0.016	-3.675	3.705	1.008	14000.80	85615.63
b_Territoriality1	0.404	-3.263	4.028	1.004	29707.15	578829.90
b_Nocturnality1	-0.120	-3.924	3.696	1.006	22353.09	561225.09
b_Habitat_Density1	0.449	-3.214	4.011	1.005	22666.55	630329.40
b_Habitat_Density2	0.408	-3.453	4.226	1.003	56983.98	752541.98
b_Migration1	-0.499	-4.291	3.326	1.001	1168381.54	754200.09
b_Migration2	0.296	-3.547	4.110	1.002	119868.53	721720.66
b_Central_Latitude	-0.533	-3.758	2.744	1.017	6325.64	45201.65

Warning: Removed 194 rows containing missing values or values outside the scale range
(`stat_slabinterval()`).

5.6 Duetting ~ White_Plumage

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 2) Intercept-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5)	Duetting ~ White_Plumage + (1 \| gr(tip.labels, cov = vcv.phylo))	10000	200	1	5000	22 (2.2e-05%)	0	24265.79	46765.81	195874151

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	0.395	-5.476	6.37	1.009	161960.86	160159.87
b_White_Plumage1	-0.028	-3.294	3.28	1.008	24265.79	46765.81

Warning: Removed 239 rows containing missing values or values outside the scale range
(`stat_slabinterval()`).

5.7 White_Plumage ~ Duetting

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 2) Intercept-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5)	White_Plumage ~ Duetting + (1 \| gr(tip.labels, cov = vcv.phylo))	10000	200	1	5000	776 (0.000776%)	0	11349.12	81841.61	788126693

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	-0.974	-8.808	5.716	1.003	59190.93	439411.63
b_Duetting1	0.004	-3.516	3.524	1.01	11349.12	81841.61

Warning: Removed 44 rows containing missing values or values outside the scale range
(`stat_slabinterval()`).

5.8 Duetting ~Social_Bond

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 2) Intercept-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5)	Duetting ~ Social_Bond + (1 \| gr(tip.labels, cov = vcv.phylo))	10000	80	1	5000	7 (1.75e-05%)	0	6680.199	21960.48	939067187

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	-0.066	-5.962	6.210	1.005	248112.232	165298.26
b_Social_Bond1	1.162	-2.256	4.350	1.006	9749.434	21960.48
b_Social_Bond2	0.076	-3.395	3.174	1.008	6680.199	45653.67

Warning: Removed 42 rows containing missing values or values outside the scale range
(`stat_slabinterval()`).

5.9 Duetting ~ Territoriality

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 2) Intercept-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5)	Duetting ~ Territoriality + (1 \| gr(tip.labels, cov = vcv.phylo))	10000	80	1	5000	16 (4e-05%)	0	173076.9	16733.73	928459507

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	-0.461	-5.710	4.743	1.005	173076.9	143994.51
b_Territoriality1	1.790	-1.028	4.671	1.004	207067.1	16733.73

Warning: Removed 163 rows containing missing values or values outside the scale range
(`stat_slabinterval()`).

5.10 Duetting ~ Nocturnality

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 2) Intercept-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5)	Duetting ~ Nocturnality + (1 \| gr(tip.labels, cov = vcv.phylo))	10000	80	1	5000	21 (5.25e-05%)	0	30944.54	18294.2	1020194856

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	-0.585	-6.139	5.272	1.008	209548.99	87619.82
b_Nocturnality1	1.300	-2.114	4.353	1.007	30944.54	18294.20

Warning: Removed 71 rows containing missing values or values outside the scale range
(`stat_slabinterval()`).

5.11 Duetting ~ Habitat_Density

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 2) Intercept-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5)	Duetting ~ Habitat_Density + (1 \| gr(tip.labels, cov = vcv.phylo))	10000	80	1	5000	7 (1.75e-05%)	0	12622.07	47167.33	1786541651

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	0.345	-5.560	6.308	1.004	207755.93	173252.72
b_Habitat_Density1	0.385	-2.619	3.496	1.005	12622.07	48535.52
b_Habitat_Density2	-1.238	-4.578	2.351	1.003	24761.27	47167.33

Warning: Removed 113 rows containing missing values or values outside the scale range
(`stat_slabinterval()`).

5.12 Duetting ~ Migration

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 2) Intercept-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5)	Duetting ~ Migration + (1 \| gr(tip.labels, cov = vcv.phylo))	10000	80	1	5000	10 (2.5e-05%)	0	14977.95	20431.69	1674747264

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	0.842	-4.190	5.814	1.005	52467.48	136958.01
b_Migration1	-1.283	-4.167	1.916	1.004	39626.47	21129.03
b_Migration2	-1.721	-4.853	1.729	1.004	14977.95	20431.69

Warning: Removed 113 rows containing missing values or values outside the scale range
(`stat_slabinterval()`).

5.13 Duetting ~ Central_Latitude

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 2) Intercept-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5)	Duetting ~ Central_Latitude + (1 \| gr(tip.labels, cov = vcv.phylo))	10000	80	1	5000	30 (7.5e-05%)	0	5906.333	17215.08	1993484585

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	0.348	-5.620	6.573	1.003	226598.904	166997.56
b_Central_Latitude	-0.833	-3.726	1.862	1.009	5906.333	17215.08

Warning: Removed 72 rows containing missing values or values outside the scale range
(`stat_slabinterval()`).

5.14 White_Plumage ~Social_Bond

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 2) Intercept-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5)	White_Plumage ~ Social_Bond + (1 \| gr(tip.labels, cov = vcv.phylo))	10000	80	1	5000	3 (7.5e-06%)	0	5511.218	66964.02	739964077

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	-1.013	-8.974	5.644	1.001	283338.504	161889.31
b_Social_Bond1	0.212	-3.310	3.742	1.008	5868.794	67353.00
b_Social_Bond2	0.003	-3.585	3.550	1.009	5511.218	66964.02

Warning: Removed 40 rows containing missing values or values outside the scale range
(`stat_slabinterval()`).

5.15 White_Plumage ~ Territoriality

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 2) Intercept-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5)	White_Plumage ~ Territoriality + (1 \| gr(tip.labels, cov = vcv.phylo))	10000	80	1	5000	5 (1.25e-05%)	0	14311.52	169418.4	724727198

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	-1.205	-8.933	5.344	1.001	253445.20	169418.4
b_Territoriality1	0.470	-2.999	3.858	1.004	14311.52	281795.3

Warning: Removed 61 rows containing missing values or values outside the scale range
(`stat_slabinterval()`).

5.16 White_Plumage ~ Nocturnality

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 2) Intercept-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5)	White_Plumage ~ Nocturnality + (1 \| gr(tip.labels, cov = vcv.phylo))	10000	80	1	5000	35 (8.75e-05%)	0	8055.175	180470.5	659207736

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	-0.789	-8.688	5.927	1.001	277757.194	180470.5
b_Nocturnality1	-0.206	-3.826	3.493	1.006	8055.175	289179.6

Warning: Removed 71 rows containing missing values or values outside the scale range
(`stat_slabinterval()`).

5.17 White_Plumage ~ Habitat_Density

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 2) Intercept-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5)	White_Plumage ~ Habitat_Density + (1 \| gr(tip.labels, cov = vcv.phylo))	10000	80	1	5000	2 (5e-06%)	0	7393.309	174724.2	1973767809

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	-1.266	-8.929	5.263	1.002	250510.953	174724.2
b_Habitat_Density1	0.493	-2.992	3.832	1.007	7393.309	207127.5
b_Habitat_Density2	0.492	-3.317	4.186	1.004	16045.716	262318.7

Warning: Removed 69 rows containing missing values or values outside the scale range
(`stat_slabinterval()`).

5.18 White_Plumage ~ Migration

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 2) Intercept-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5)	White_Plumage ~ Migration + (1 \| gr(tip.labels, cov = vcv.phylo))	10000	80	1	5000	2 (5e-06%)	0	40365.97	169679.6	763570319

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	-0.886	-8.326	5.237	1.001	302064.71	169679.6
b_Migration1	-0.704	-4.334	3.033	1.001	168037.93	296337.3
b_Migration2	0.301	-3.403	3.951	1.003	40365.97	299018.4

Warning: Removed 80 rows containing missing values or values outside the scale range
(`stat_slabinterval()`).

5.19 White_Plumage ~ Central_Latitude

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 2) Intercept-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5)	White_Plumage ~ Central_Latitude + (1 \| gr(tip.labels, cov = vcv.phylo))	10000	80	1	5000	6 (1.5e-05%)	0	3000.62	21275.02	1573392887

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	-0.912	-8.628	5.445	1.001	323996.82	172921.67
b_Central_Latitude	-0.511	-3.548	2.518	1.015	3000.62	21275.02

6 Using a single consensus tree

6.1 Full model

Code

### See https://urldefense.proofpoint.com/v2/url?u=https-3A__psyarxiv.com_x8swp_&d=DwIGAg&c=t1KneP5CA95qj21Jnj2QJw&r=uXrXEa22P1XTtZOkZJzsx8pAQCDanQ2sMUk_AfMLQ0I&m=Th_rFO1G4mV3f039Te1VuvpmcjwF2KXjHAfcU0UhAb8vpbbU-nEL3HwCTlICHqS3&s=qolX0qBhJIjMDpanzBhiCR3QCeDrzGH4A9eyiqlhW2g&e=  for explanation of types of ordinal models implemented in brms.
# And https://urldefense.proofpoint.com/v2/url?u=https-3A__solomonkurz.netlify.app_post_2021-2D12-2D29-2Dnotes-2Don-2Dthe-2Dbayesian-2Dcumulative-2Dprobit_&d=DwIGAg&c=t1KneP5CA95qj21Jnj2QJw&r=uXrXEa22P1XTtZOkZJzsx8pAQCDanQ2sMUk_AfMLQ0I&m=Th_rFO1G4mV3f039Te1VuvpmcjwF2KXjHAfcU0UhAb8vpbbU-nEL3HwCTlICHqS3&s=PQ-Sgi8kE3P2Yamj_rtG-r2rZUbEB5iytYjn9NLil08&e= 

# set trees
#A <- vcv.phylo(trees)

# set parameters
set.seed(123)
parallel <-
    7 # set this to the number of cores you want to use in the cluster
thin <- 1
priors <- c(prior(normal(0, 2), class = "b"))
chains <- 4
A <- ape::vcv.phylo(consensus_tree_owl)

df$phylo <- df$tip.label
## run model

model_simple <- brm(
  Duetting ~ Social_Bond + Territoriality + Nocturnality + Habitat_Density + Migration + Central_Latitude + White_Plumage + (1|gr(phylo, cov = A)),
  data = df,
    family = bernoulli(link = "logit"),
    chains = chains,
    cores = chains,
    iter = 20000,
    prior = priors,
    control = list(adapt_delta = 0.99, max_treedepth = 15),
    init = 0,
    thin = thin,
  data2 = list(A = A),
)

saveRDS(model_simple, "./data/processed/consensus_tree_model.RDS")

Code

extended_summary(read.file = "./data/processed/consensus_tree_model.RDS", plot.area.prop = 0.5)

Warning: Removed 29 rows containing missing values or values outside the scale range
(`stat_slabinterval()`).

6.2 consensus_tree_model

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 2) Intercept-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5)	Duetting ~ Social_Bond + Territoriality + Nocturnality + Habitat_Density + Migration + Central_Latitude + White_Plumage + (1 \| gr(phylo, cov = A))	20000	4	1	10000	1 (2.5e-05%)	0	25702.21	21990.09	2068286766

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	-1.320	-9.063	6.861	1	32485.96	21990.09
b_Social_Bond1	0.961	-2.614	4.452	1	40780.90	29493.88
b_Social_Bond2	-0.349	-3.789	3.187	1	37364.62	25765.24
b_Territoriality1	1.375	-2.219	4.770	1	34832.39	28967.53
b_Nocturnality1	0.870	-2.832	4.524	1	43053.88	29800.62
b_Habitat_Density1	0.445	-2.963	3.765	1	36199.48	25325.08
b_Habitat_Density2	-0.851	-4.510	2.863	1	45939.89	30815.65
b_Migration1	-0.807	-4.353	2.822	1	37535.51	25977.66
b_Migration2	-1.127	-4.808	2.672	1	39574.19	27112.58
b_Central_Latitude	-0.553	-3.694	2.276	1	25702.21	28882.91
b_White_Plumage1	-0.087	-3.656	3.543	1	41341.44	30498.61

6.3 Specific models

Code

# set parameters
set.seed(123)
parallel <- 1
thin <- 1
priors <- c(prior(normal(0, 2), class = "b"))
chains <- 2


## run model

### Owl duets: present/absent vs explanatory variables

for (i in models) {
    
    print(i)
    mod_name <- gsub("~", "_by_", i)
    mod_name <- gsub(" \\+ ", "_", mod_name)
    mod_name <- paste0(mod_name, "consensus_tree.RDS")
    mod_name <- file.path("./data/processed", mod_name)

    if (!file.exists(mod_name)){
    A <- ape::vcv.phylo(consensus_tree_owl)

        df$phylo <- df$tip.label
## run model

mod <- brm(
    formula = as.formula(paste(i, " + (1|gr(phylo, cov = A))")),
  data = df,
    family = bernoulli(link = "logit"),
    chains = chains,
    cores = chains,
    iter = 20000,
    prior = priors,
    control = list(adapt_delta = 0.99, max_treedepth = 15),
    init = 0,
    thin = thin,
  data2 = list(A = A),
)
    

    saveRDS(mod, mod_name)
    rm("mod")
    }
}

Code

# mod <- readRDS("./data/processed/phylogenetic_uncertainty_model.RDS")

for (i in models) {

     mod_name <- gsub("~", "_by_", i)
    mod_name <- gsub(" \\+ ", "_", mod_name)
    mod_name <- paste0(mod_name, "consensus_tree.RDS")
    mod_name <- file.path("./data/processed", mod_name)
    
 extended_summary(read.file = mod_name, plot.area.prop = 0.5, fit.name = i, highlight = TRUE)
}

Warning: Removed 83 rows containing missing values or values outside the scale range
(`stat_slabinterval()`).

6.4 Duetting ~Social_Bond + Territoriality + Nocturnality + Habitat_Density + Migration + Central_Latitude

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 2) Intercept-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5)	Duetting ~ Social_Bond + Territoriality + Nocturnality + Habitat_Density + Migration + Central_Latitude + (1 \| gr(phylo, cov = A))	20000	2	1	10000	0 (0%)	0	19384.81	12647.8	284639929

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	-1.415	-9.247	6.658	1	23802.78	12647.80
b_Social_Bond1	0.970	-2.505	4.422	1	30787.82	15905.94
b_Social_Bond2	-0.355	-3.792	3.068	1	30986.53	15675.11
b_Territoriality1	1.407	-2.146	4.776	1	26566.96	16638.31
b_Nocturnality1	0.900	-2.795	4.516	1	30696.68	14717.45
b_Habitat_Density1	0.451	-2.946	3.744	1	27398.94	15803.71
b_Habitat_Density2	-0.870	-4.531	2.828	1	34158.01	14765.10
b_Migration1	-0.832	-4.279	2.740	1	25514.87	16887.09
b_Migration2	-1.143	-4.865	2.605	1	26502.08	15324.90
b_Central_Latitude	-0.507	-3.584	2.254	1	19384.81	15985.84

Warning: Removed 21 rows containing missing values or values outside the scale range
(`stat_slabinterval()`).

6.5 White_Plumage ~ Social_Bond + Territoriality + Nocturnality + Habitat_Density + Migration + Central_Latitude

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 2) Intercept-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5)	White_Plumage ~ Social_Bond + Territoriality + Nocturnality + Habitat_Density + Migration + Central_Latitude + (1 \| gr(phylo, cov = A))	20000	2	1	10000	0 (0%)	0	14993.37	8862.774	962619197

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	-0.867	-10.315	7.426	1	16612.47	8862.774
b_Social_Bond1	0.057	-3.649	3.789	1	20135.63	15966.641
b_Social_Bond2	0.025	-3.694	3.695	1	20123.52	15688.070
b_Territoriality1	0.409	-3.272	4.000	1	19566.13	15500.040
b_Nocturnality1	-0.270	-4.058	3.584	1	21700.46	14992.959
b_Habitat_Density1	0.352	-3.394	3.990	1	19085.87	15063.381
b_Habitat_Density2	0.283	-3.540	4.064	1	21773.23	15498.388
b_Migration1	-0.496	-4.331	3.386	1	21530.07	14927.408
b_Migration2	0.240	-3.592	4.156	1	21571.41	15064.492
b_Central_Latitude	-0.503	-3.741	2.819	1	14993.37	14548.308

Warning: Removed 40 rows containing missing values or values outside the scale range
(`stat_slabinterval()`).

6.6 Duetting ~ White_Plumage

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 2) Intercept-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5)	Duetting ~ White_Plumage + (1 \| gr(phylo, cov = A))	20000	2	1	10000	0 (0%)	0	12992.49	8749.731	1368161037

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	0.224	-6.370	7.109	1	12992.49	8749.731
b_White_Plumage1	0.000	-3.381	3.464	1	17979.86	14585.969

Warning: Removed 74 rows containing missing values or values outside the scale range
(`stat_slabinterval()`).

6.7 White_Plumage ~ Duetting

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 2) Intercept-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5)	White_Plumage ~ Duetting + (1 \| gr(phylo, cov = A))	20000	2	1	10000	0 (0%)	0	11646.6	7734.528	2137258029

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	-0.798	-8.922	6.203	1	11646.60	7734.528
b_Duetting1	0.009	-3.633	3.596	1	16373.55	14646.132

Warning: Removed 21 rows containing missing values or values outside the scale range
(`stat_slabinterval()`).

6.8 Duetting ~Social_Bond

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 2) Intercept-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5)	Duetting ~ Social_Bond + (1 \| gr(phylo, cov = A))	20000	2	1	10000	1 (5e-05%)	0	11639.24	9862.38	879040948

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	-0.202	-6.558	6.497	1	13320.91	9862.38
b_Social_Bond1	1.244	-2.283	4.551	1	13707.97	13168.16
b_Social_Bond2	-0.072	-3.581	3.154	1	11639.24	13644.80

Warning: Removed 18 rows containing missing values or values outside the scale range
(`stat_slabinterval()`).

6.9 Duetting ~ Territoriality

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 2) Intercept-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5)	Duetting ~ Territoriality + (1 \| gr(phylo, cov = A))	20000	2	1	10000	1 (5e-05%)	0	11267.19	7785.08	396010916

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	-0.631	-6.420	5.246	1	11267.19	7785.08
b_Territoriality1	1.708	-1.473	4.686	1	11803.72	10537.49

Warning: Removed 22 rows containing missing values or values outside the scale range
(`stat_slabinterval()`).

6.10 Duetting ~ Nocturnality

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 2) Intercept-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5)	Duetting ~ Nocturnality + (1 \| gr(phylo, cov = A))	20000	2	1	10000	0 (0%)	0	12224.85	8966.268	1575032726

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	-0.795	-6.886	5.77	1	12531.77	8966.268
b_Nocturnality1	1.337	-2.231	4.61	1	12224.85	10797.038

Warning: Removed 45 rows containing missing values or values outside the scale range
(`stat_slabinterval()`).

6.11 Duetting ~ Habitat_Density

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 2) Intercept-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5)	Duetting ~ Habitat_Density + (1 \| gr(phylo, cov = A))	20000	2	1	10000	0 (0%)	0	16302.21	10734.33	1129355307

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	0.220	-5.996	6.611	1	16302.21	10734.33
b_Habitat_Density1	0.306	-2.850	3.416	1	18803.72	12981.83
b_Habitat_Density2	-1.293	-4.691	2.379	1	17247.18	12233.93

Warning: Removed 27 rows containing missing values or values outside the scale range
(`stat_slabinterval()`).

Warning: Removed 1 row containing missing values or values outside the scale range
(`stat_slabinterval()`).

6.12 Duetting ~ Migration

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 2) Intercept-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5)	Duetting ~ Migration + (1 \| gr(phylo, cov = A))	20000	2	1	10000	0 (0%)	0	12972.23	8086.402	994363304

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	0.598	-5.021	5.876	1	12972.23	8086.402
b_Migration1	-1.182	-4.221	2.149	1	13636.41	11860.885
b_Migration2	-1.787	-5.034	1.748	1	13467.61	10445.956

Warning: Removed 43 rows containing missing values or values outside the scale range
(`stat_slabinterval()`).

6.13 Duetting ~ Central_Latitude

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 2) Intercept-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5)	Duetting ~ Central_Latitude + (1 \| gr(phylo, cov = A))	20000	2	1	10000	1 (5e-05%)	0	9416.819	7932.898	1473776745

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	0.174	-6.356	6.785	1	12667.038	7932.898
b_Central_Latitude	-0.929	-3.917	1.922	1	9416.819	11396.119

Warning: Removed 23 rows containing missing values or values outside the scale range
(`stat_slabinterval()`).

6.14 White_Plumage ~Social_Bond

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 2) Intercept-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5)	White_Plumage ~ Social_Bond + (1 \| gr(phylo, cov = A))	20000	2	1	10000	0 (0%)	0	20467.94	10116.42	1264837356

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	-0.799	-9.150	6.112	1	20467.94	10116.42
b_Social_Bond1	0.145	-3.465	3.708	1	27730.06	16044.76
b_Social_Bond2	0.050	-3.561	3.641	1	29166.55	15331.93

Warning: Removed 17 rows containing missing values or values outside the scale range
(`stat_slabinterval()`).

6.15 White_Plumage ~ Territoriality

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 2) Intercept-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5)	White_Plumage ~ Territoriality + (1 \| gr(phylo, cov = A))	20000	2	1	10000	0 (0%)	0	24145.79	10901.86	717072421

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	-1.076	-8.895	5.599	1	24145.79	10901.86
b_Territoriality1	0.541	-3.010	3.980	1	26242.72	16704.93

Warning: Removed 78 rows containing missing values or values outside the scale range
(`stat_slabinterval()`).

6.16 White_Plumage ~ Nocturnality

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 2) Intercept-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5)	White_Plumage ~ Nocturnality + (1 \| gr(phylo, cov = A))	20000	2	1	10000	1 (5e-05%)	0	17469.54	10582.14	1056711716

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	-0.499	-8.367	6.481	1	17469.54	10582.14
b_Nocturnality1	-0.300	-3.996	3.414	1	24293.85	15397.24

Warning: Removed 29 rows containing missing values or values outside the scale range
(`stat_slabinterval()`).

6.17 White_Plumage ~ Habitat_Density

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 2) Intercept-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5)	White_Plumage ~ Habitat_Density + (1 \| gr(phylo, cov = A))	20000	2	1	10000	0 (0%)	0	16403.45	8691.624	2017218430

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	-1.000	-8.919	5.585	1	16403.45	8691.624
b_Habitat_Density1	0.372	-3.199	3.799	1	18176.38	15145.540
b_Habitat_Density2	0.361	-3.463	4.072	1	22583.70	13959.035

Warning: Removed 44 rows containing missing values or values outside the scale range
(`stat_slabinterval()`).

Warning: Removed 2 rows containing missing values or values outside the scale range
(`stat_slabinterval()`).

6.18 White_Plumage ~ Migration

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 2) Intercept-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5)	White_Plumage ~ Migration + (1 \| gr(phylo, cov = A))	20000	2	1	10000	0 (0%)	0	21407.36	9365.559	2069551840

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	-0.715	-8.468	5.787	1	21407.36	9365.559
b_Migration1	-0.641	-4.309	3.135	1	29752.83	15026.168
b_Migration2	0.214	-3.493	3.911	1	32718.15	15996.854

Warning: Removed 16 rows containing missing values or values outside the scale range
(`stat_slabinterval()`).

6.19 White_Plumage ~ Central_Latitude

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 2) Intercept-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5)	White_Plumage ~ Central_Latitude + (1 \| gr(phylo, cov = A))	20000	2	1	10000	0 (0%)	0	14256.11	7863.32	422709779

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	-0.749	-8.981	6.004	1.001	16649.11	7863.32
b_Central_Latitude	-0.483	-3.556	2.776	1	14256.11	13223.96

Takeaways

No significant effects but a important trend in territoriality

Session information

R version 4.1.2 (2021-11-01)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 22.04.4 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0

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

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

other attached packages:
 [1] rstan_2.32.6        StanHeaders_2.32.6  phytools_2.1-1     
 [4] maps_3.4.2          corrplot_0.92       brmsish_1.0.0      
 [7] ape_5.7-1           brms_2.21.0         Rcpp_1.0.12        
[10] viridis_0.6.5       viridisLite_0.4.2   xaringanExtra_0.7.0
[13] rprojroot_2.0.4     formatR_1.14        knitr_1.46         
[16] kableExtra_1.4.0    ggplot2_3.5.1      

loaded via a namespace (and not attached):
  [1] TH.data_1.1-0           colorspace_2.1-0        estimability_1.5.1     
  [4] QuickJSR_1.1.3          rstudioapi_0.15.0       farver_2.1.2           
  [7] remotes_2.5.0           svUnit_1.0.6            optimParallel_1.0-2    
 [10] fansi_1.0.6             mvtnorm_1.2-4           xml2_1.3.6             
 [13] bridgesampling_1.1-2    codetools_0.2-18        splines_4.1.2          
 [16] mnormt_2.1.1            doParallel_1.0.17       bayesplot_1.11.1       
 [19] jsonlite_1.8.8          packrat_0.9.2           ggdist_3.3.2           
 [22] compiler_4.1.2          emmeans_1.10.2          backports_1.4.1        
 [25] Matrix_1.6-5            fastmap_1.1.1           cli_3.6.2              
 [28] htmltools_0.5.8.1       tools_4.1.2             igraph_2.0.3           
 [31] coda_0.19-4.1           gtable_0.3.5            glue_1.7.0             
 [34] reshape2_1.4.4          clusterGeneration_1.3.8 dplyr_1.1.4            
 [37] posterior_1.5.0         V8_4.4.2                fastmatch_1.1-4        
 [40] vctrs_0.6.5             svglite_2.1.3           nlme_3.1-155           
 [43] iterators_1.0.14        tensorA_0.36.2.1        xfun_0.43              
 [46] stringr_1.5.1           lifecycle_1.0.4         phangorn_2.11.1        
 [49] tidybayes_3.0.6         MASS_7.3-55             zoo_1.8-12             
 [52] scales_1.3.0            Brobdingnag_1.2-9       parallel_4.1.2         
 [55] sandwich_3.0-1          inline_0.3.19           expm_0.999-9           
 [58] yaml_2.3.8              curl_5.2.1              pbapply_1.7-2          
 [61] gridExtra_2.3           loo_2.7.0               stringi_1.8.4          
 [64] foreach_1.5.2           checkmate_2.3.1         pkgbuild_1.4.4         
 [67] rlang_1.1.4             pkgconfig_2.0.3         systemfonts_1.0.6      
 [70] matrixStats_1.2.0       distributional_0.4.0    evaluate_0.23          
 [73] lattice_0.20-45         purrr_1.0.2             labeling_0.4.3         
 [76] rstantools_2.4.0        htmlwidgets_1.6.4       cowplot_1.1.3          
 [79] tidyselect_1.2.1        plyr_1.8.9              magrittr_2.0.3         
 [82] R6_2.5.1                generics_0.1.3          multcomp_1.4-18        
 [85] combinat_0.0-8          pillar_1.9.0            withr_3.0.0            
 [88] scatterplot3d_0.3-44    survival_3.2-13         abind_1.4-5            
 [91] tibble_3.2.1            crayon_1.5.2            arrayhelpers_1.1-0     
 [94] utf8_1.2.4              rmarkdown_2.26          grid_4.1.2             
 [97] sketchy_1.0.3           digest_0.6.35           xtable_1.8-4           
[100] tidyr_1.3.1             numDeriv_2016.8-1.1     RcppParallel_5.1.7     
[103] stats4_4.1.2            munsell_0.5.1           quadprog_1.5-8