Owl duetting analysis

Owl duetting evolution

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",
        github = "maRce10/brmsish"
    )

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

Code

## Open dataframes

# read data
trees_owl <- read.nexus("./data/raw/BirdTree_Owl_Output.nex")
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)

2.0.1 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.325715 
P-value (based on 1000 randomizations) : 0.023 

Code

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

Phylogenetic signal lambda : 0.428795 
logL(lambda) : -37.3961 
LR(lambda=0) : 1.65596 
P-value (based on LR test) : 0.19815 

2.0.2 PHYLOGENETICALLY CONTROLLED COMPARATIVE ANALYSES USING JETZ TREE

2.0.3 USE TO COMPARE DUET PRESENCE/ABSENCE TO LIFE HISTORY TRAITS

2.1 BRMS V1 - original

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 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)
}

# 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

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

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

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

2.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	47	0	27358.26	83531.01	38507154

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	-1.292	-8.937	6.678	1.002	128141.86	997518.31
b_Social_Bond1	0.911	-2.618	4.409	1.006	34934.72	135293.54
b_Social_Bond2	-0.328	-3.787	3.146	1.006	39456.17	143891.96
b_Territoriality1	1.582	-1.975	4.957	1.008	27566.95	83531.01
b_Nocturnality1	0.820	-2.894	4.453	1.001	288782.17	2794811.29
b_Habitat_Density1	0.550	-2.845	3.871	1.006	39240.25	164944.73
b_Habitat_Density2	-0.798	-4.458	2.897	1.001	376886.05	2950301.26
b_Migration1	-0.893	-4.391	2.724	1.002	158111.77	769514.37
b_Migration2	-1.092	-4.756	2.678	1.003	94908.56	691197.39
b_Central_Latitude	-0.522	-3.623	2.228	1.008	27358.26	86425.84
b_White_Plumage1	-0.126	-3.633	3.422	1.004	55252.75	333703.97

Takeaways

 

---

 

Session information

R version 4.3.2 (2023-10-31)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.2 LTS

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

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

time zone: America/Costa_Rica
tzcode source: system (glibc)

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

other attached packages:
 [1] phytools_2.1-1      maps_3.4.2          corrplot_0.92      
 [4] brmsish_1.0.0       ape_5.7-1           brms_2.20.4        
 [7] Rcpp_1.0.12         viridis_0.6.5       viridisLite_0.4.2  
[10] xaringanExtra_0.7.0 rprojroot_2.0.4     formatR_1.14       
[13] knitr_1.45          kableExtra_1.4.0    ggplot2_3.5.0      

loaded via a namespace (and not attached):
  [1] tensorA_0.36.2.1        rstudioapi_0.15.0       jsonlite_1.8.8         
  [4] magrittr_2.0.3          TH.data_1.1-2           sketchy_1.0.2          
  [7] estimability_1.4.1      farver_2.1.1            rmarkdown_2.25         
 [10] vctrs_0.6.5             base64enc_0.1-3         htmltools_0.5.7        
 [13] distributional_0.3.2    curl_5.2.0              StanHeaders_2.32.5     
 [16] htmlwidgets_1.6.4       plyr_1.8.9              sandwich_3.1-0         
 [19] emmeans_1.10.0          zoo_1.8-12              igraph_1.6.0           
 [22] mime_0.12               lifecycle_1.0.4         iterators_1.0.14       
 [25] pkgconfig_2.0.3         colourpicker_1.3.0      Matrix_1.6-5           
 [28] R6_2.5.1                fastmap_1.1.1           shiny_1.8.0            
 [31] digest_0.6.34           numDeriv_2016.8-1.1     colorspace_2.1-0       
 [34] crosstalk_1.2.1         labeling_0.4.3          clusterGeneration_1.3.8
 [37] fansi_1.0.6             abind_1.4-5             compiler_4.3.2         
 [40] remotes_2.4.2.1         withr_3.0.0             doParallel_1.0.17      
 [43] backports_1.4.1         inline_0.3.19           shinystan_2.6.0        
 [46] optimParallel_1.0-2     QuickJSR_1.1.3          pkgbuild_1.4.3         
 [49] MASS_7.3-60.0.1         scatterplot3d_0.3-44    gtools_3.9.5           
 [52] loo_2.6.0               tools_4.3.2             httpuv_1.6.14          
 [55] threejs_0.3.3           quadprog_1.5-8          glue_1.7.0             
 [58] nlme_3.1-164            promises_1.2.1          grid_4.3.2             
 [61] checkmate_2.3.1         reshape2_1.4.4          generics_0.1.3         
 [64] gtable_0.3.4            xml2_1.3.6              utf8_1.2.4             
 [67] foreach_1.5.2           pillar_1.9.0            ggdist_3.3.1           
 [70] markdown_1.12           stringr_1.5.1           posterior_1.5.0        
 [73] later_1.3.2             splines_4.3.2           dplyr_1.1.4            
 [76] lattice_0.22-5          survival_3.5-7          tidyselect_1.2.0       
 [79] miniUI_0.1.1.1          pbapply_1.7-2           gridExtra_2.3          
 [82] V8_4.4.2                svglite_2.1.3           stats4_4.3.2           
 [85] xfun_0.42               expm_0.999-9            bridgesampling_1.1-2   
 [88] matrixStats_1.2.0       DT_0.31                 rstan_2.32.5           
 [91] stringi_1.8.3           yaml_2.3.8              evaluate_0.23          
 [94] codetools_0.2-19        tibble_3.2.1            cli_3.6.2              
 [97] RcppParallel_5.1.7      shinythemes_1.2.0       xtable_1.8-4           
[100] systemfonts_1.0.5       munsell_0.5.0           coda_0.19-4.1          
[103] parallel_4.3.2          rstantools_2.4.0        ellipsis_0.3.2         
[106] dygraphs_1.1.1.6        bayesplot_1.11.0        Brobdingnag_1.2-9      
[109] phangorn_2.11.1         mvtnorm_1.2-4           scales_1.3.0           
[112] xts_0.13.2              packrat_0.9.2           crayon_1.5.2           
[115] combinat_0.0-8          rlang_1.1.3             fastmatch_1.1-4        
[118] cowplot_1.1.3           multcomp_1.4-25         mnormt_2.1.1           
[121] shinyjs_2.1.0