Mobbing ontogeny

Published

January 18, 2024

Source code and data found at https://github.com/maRce10/mobbing_ontogeny

Cargar paquetes

Code

sketchy::load_packages(c(
  "googlesheets4",
  "httpuv",
  "ggplot2",
  "kableExtra",
  "viridis",
  "brms",
  github = "maRce10/brmsish",
  "emmeans"
))

.print_df <- function(x, highlight = NULL, height = "300px",...) {
  kbl <- kableExtra::kable(
    as.data.frame(x)
    ,
    align = "c",
    row.names = F,
    format = "html",
    escape = F
  )
  
  if (!is.null(highlight))
  kbl <- column_spec(kbl, column = which(names(x) %in% highlight), background = "#ccebff", bold = TRUE)
    
  kbl <-
    kableExtra::kable_styling(kbl, bootstrap_options = "striped", font_size = 14)

  
  kbl <-
    kableExtra::scroll_box(kbl, width = "800px", height = height)
  
  knitr::asis_output(kbl)
}

1 Predicciones

1.1 Aprendizaje

Code

acoso <- c(4, 2, 2, 0)
estimulo <- c("Simpátrico", "Alopátrico", "Alopátrico +", "Control")

sim_df <- data.frame(estimulo, acoso)

set.seed(123)
out <- lapply(seq_len(nrow(sim_df)), function(x){
  
  acoso <- rnorm(n = 30, mean = sim_df$acoso[x])
  individuos <- rpois(n = 30, lambda = sim_df$acoso[x] / 2)

  
  data.frame(estimulo = sim_df$estimulo[x], acoso, individuos)
})

sim_data <- do.call(rbind, out)

sim_data$estimulo <- factor(sim_data$estimulo, levels = c("Simpátrico", "Alopátrico", "Alopátrico +", "Control"))

ggplot(sim_data, aes(x = estimulo, y = acoso)) +
  geom_boxplot(outlier.alpha = 0, fill = viridis(10, alpha = 0.6)[7]) +
  labs(x = "Estímulo", y = "Intensidad del acoso") + 
  theme_classic(base_size = 24) +
    theme(axis.text.y=element_blank(),
        axis.ticks.y=element_blank())

Code

agg_data <- data.frame(estimulo = c("Simpátrico", "Alopátrico", "Alopátrico +", "Control"), individuos = c(15, 9, 9, 3))

agg_data$estimulo <- factor(agg_data$estimulo, levels = c("Simpátrico", "Alopátrico", "Alopátrico +", "Control"))

ggplot(agg_data, aes(x = estimulo, y = individuos)) +
  geom_bar(stat="identity", fill = viridis(10, alpha = 0.4)[7], color = "black") +
  labs(x = "Estímulo", y = "Número de individuos") + 
  theme_classic(base_size = 24) +
    theme(axis.text.y=element_blank(),
        axis.ticks.y=element_blank())

1.2 Innato

Code

acoso <- c(4, 4, 4, 1)
estimulo <- c("Simpátrico", "Alopátrico", "Alopátrico +", "Control")

sim_df <- data.frame(estimulo, acoso)

set.seed(1)
out <- lapply(seq_len(nrow(sim_df)), function(x){
  
  acoso <- rnorm(n = 30, mean = sim_df$acoso[x])
  
  individuos <- rpois(n = 30, lambda = sim_df$acoso[x] / 2)

  data.frame(estimulo = sim_df$estimulo[x], acoso, individuos)
})

sim_data <- do.call(rbind, out)

sim_data$estimulo <- factor(sim_data$estimulo, levels = c("Simpátrico", "Alopátrico", "Alopátrico +", "Control"))

ggplot(sim_data, aes(x = estimulo, y = acoso)) +
  geom_boxplot(outlier.alpha = 0, fill = viridis(10, alpha = 0.6)[7]) +
  labs(x = "Estímulo", y = "Intensidad del acoso") + 
  theme_classic(base_size = 24) +
    theme(axis.text.y=element_blank(),
        axis.ticks.y=element_blank())

Code

agg_data <- data.frame(estimulo = c("Simpátrico", "Alopátrico", "Alopátrico +", "Control"), individuos = c(4, 4, 4, 1))

agg_data$estimulo <- factor(agg_data$estimulo, levels = c("Simpátrico", "Alopátrico", "Alopátrico +", "Control"))

ggplot(agg_data, aes(x = estimulo, y = individuos)) +
  geom_bar(stat="identity", fill = viridis(10, alpha = 0.4)[7], color = "black") +
  labs(x = "Estímulo", y = "Número de individuos") + 
  theme_classic(base_size = 24) +
    theme(axis.text.y=element_blank(),
        axis.ticks.y=element_blank())

2 Leer datos desde google docs

Code

mob_data <- read_sheet('https://docs.google.com/spreadsheets/d/19MjQaAXfoMa59AchxTL4RC0TjGLVwHqvkotBXwn8G30/edit#gid=1949444895', na = "NA")

mob_data$HORA <- format(as.POSIXct(mob_data$HORA), format = "%H:%M")

.print_df(mob_data, height = "600px")

FECHA	SITIO	EQUIPO	HORA	PUNTO	ESTIMULO	INDIVIDUOS	MORFOESPECIE	OBSERVACIONES
2024-01-14	LAS_CRUCES	1	05:50	1	GLBR	0	0	En un claro. Presencia de neblina.
2024-01-14	LAS_CRUCES	1	06:11	2	GLCR	4	2	Poca neblina. Borde del camino.
2024-01-14	LAS_CRUCES	1	06:32	3	BADE	0	0	Cielo despejado. Borde del camino.
2024-01-14	LAS_CRUCES	1	06:50	4	GLEU	0	0	Dosel más o menos cerrado. Borde del camino. Sonido de cigarras
2024-01-14	LAS_CRUCES	1	07:19	5	GLEU	1	1	Área de dosel más o menos cerrado. Borde del camino.
2024-01-14	LAS_CRUCES	1	07:36	6	GLBR	2	2	La higuera. En un claro. Sonido de cigarras
2024-01-14	LAS_CRUCES	1	07:55	7	BADE	0	0	Dosel muy cerrado. Borde del camino. Sonido de cigarras.
2024-01-14	LAS_CRUCES	1	08:19	8	GLCR	2	2	Cruzando el río. Terreno con pendiente, en el límite de un claro. Sonido de cigarras
2024-01-14	LAS_CRUCES	1	08:40	9	GLCR	14	6	Zona plana. Dosel cerrado. Soleado.
2024-01-14	LAS_CRUCES	1	09:05	10	GLEU	3	3	Dosel cerrado. Soleado. Poca pendiente.
2024-01-14	LAS_CRUCES	1	09:25	11	GLBR	2	1	Dosel cerrado. Cielo depejado. Sonido de cigarras. Cerca de un barranco.
2024-01-14	LAS_CRUCES	1	09:44	12	BADE	0	0	Cielo despejado. Poca inclinación. Sonido de cigarras.
2024-01-14	LAS_CRUCES	2	06:06	1	GLCR	0	0	SE GRABÓ UN MINUTO DESPUÉS
2024-01-14	LAS_CRUCES	2	06:27	2	GLEU	0	0	SE COLOCÓ EL PARLANTE A 125M DEL PUNTO ANTERIOR, POR ACCESIBILIDAD DE VÍA.
2024-01-14	LAS_CRUCES	2	06:42	3	GLBR	0	0	NA
2024-01-14	LAS_CRUCES	2	06:58	4	BADE	0	0	BANDADA MIXTA A LO LEJOS
2024-01-14	LAS_CRUCES	2	07:12	5	GLEU	1	1	SE ACERCÓ A MENOS DE UN METRO, MUCHO RUIDO DE CHICHARRAS
2024-01-14	LAS_CRUCES	2	07:27	6	GLBR	0	0	SOLEADO, MUCHAS CHICHARRAS
2024-01-14	LAS_CRUCES	2	07:39	7	BADE	0	0	NA
2024-01-14	LAS_CRUCES	2	07:51	8	GLCR	0	0	NA
2024-01-14	LAS_CRUCES	2	08:04	9	GLBR	0	0	NA
2024-01-14	LAS_CRUCES	2	08:18	10	GLCR	0	0	NA
2024-01-14	LAS_CRUCES	2	08:33	11	BADE	0	0	NA
2024-01-14	LAS_CRUCES	2	08:46	12	GLEU	0	0	NA
2024-01-14	LAS_CRUCES	2	09:00	13	BADE	0	0	NA
2024-01-14	LAS_CRUCES	2	09:16	14	GLBR	0	0	NA
2024-01-14	LAS_CRUCES	2	15:34	15	GLEU	0	0	Realizado por la tarde. Proximidad a la calle
2024-01-14	LAS_CRUCES	2	15:50	16	GLCR	1	1	Junto al bosque
2024-01-14	LAS_CRUCES	2	16:15	17	GLCR	6	4	Dosel cerrado
2024-01-14	LAS_CRUCES	2	16:38	18	BADE	0	0	Orilla del rio
2024-01-14	LAS_CRUCES	2	17:17	19	GLBR	2	1	Dosel mas abierto. Ruido de autos
2024-01-15	LAS_CRUCES	1	05:43	1	GLCR	3	2	Aumentó la vocalización
2024-01-15	LAS_CRUCES	1	05:58	2	GLEU	2	2	Sitio con pendiente
2024-01-15	LAS_CRUCES	1	06:10	3	BADE	0	0	NA
2024-01-15	LAS_CRUCES	1	06:35	4	GLBR	2	2	Aumentó la vocalización
2024-01-15	LAS_CRUCES	2	05:50	2	BADE	0	0	NA
2024-01-15	LAS_CRUCES	2	06:02	3	GLBR	1	1	Vocalizaciones cercanas al parlante
2024-01-15	LAS_CRUCES	2	06:16	4	GLEU	3	2	Arremon se acercó mucho al parlante
2024-01-15	LAS_CRUCES	2	06:35	5	GLEU	0	0	NA
2024-01-15	LAS_CRUCES	2	06:47	6	GLBR	0	0	NA

3 Descripción de los datos

40 puntos de muestreo
49 individuos observados
Un máximo de 33 morfoespecies observadas
40% de muestreos tuvieron respuesta (62% de los muestreos del equipo 1 y 25% de los muestreos del equipo 2)

4 Resultados

4.1 Numero de individuos

Code

mob_data$estimulo <- factor(mob_data$ESTIMULO)

mob_data$estimulo <- factor(mob_data$estimulo, labels = c("Control", "Alopátrico", "Simpátrico",  "Alopátrico +"))

mob_data$estimulo <- factor(mob_data$estimulo, levels = c("Simpátrico", "Alopátrico", "Alopátrico +", "Control"))

ggplot(mob_data, aes(x = estimulo, y = INDIVIDUOS)) +
  geom_boxplot(outlier.alpha = 0, fill = viridis(10, alpha = 0.6)[7]) +
  labs(x = "Estímulo", y = "Número de individuos") + 
  theme_classic(base_size = 24) + ylim(c(0, 8))

Warning: Removed 1 rows containing non-finite values (`stat_boxplot()`).

Code

agg_data <- aggregate(INDIVIDUOS ~ estimulo, data = mob_data, FUN = sum)

ggplot(agg_data, aes(x = estimulo, y = INDIVIDUOS)) +
  geom_bar(stat="identity", fill = viridis(10, alpha = 0.6)[7], color = "black") +
  labs(x = "Estímulo", y = "Número de individuos") + 
  theme_classic(base_size = 24)

4.2 Numero de morfoespecies

Code

ggplot(mob_data, aes(x = estimulo, y = MORFOESPECIE)) +
  geom_boxplot(outlier.alpha = 0, fill = viridis(10, alpha = 0.6)[7]) +
  labs(x = "Estímulo", y = "Número de morfoespecies") + 
  theme_classic(base_size = 24)

Code

agg_data <- aggregate(MORFOESPECIE ~ estimulo, data = mob_data, FUN = sum)

agg_data$MORFOESPECIE <- c(9, 5, 6, 0)

ggplot(agg_data, aes(x = estimulo, y = MORFOESPECIE)) +
  geom_bar(stat="identity", fill = viridis(10, alpha = 0.6)[7], color = "black") +
  labs(x = "Estímulo", y = "Número de morfoespecies") + 
  theme_classic(base_size = 24)

5 Análisis estadístico

5.1 Numero de individuos

Code

iter <- 5000
chains <- 4
priors <- c(prior(normal(0, 6), class = "b"))

set.seed(123)

mod_indiv <- brm(formula = INDIVIDUOS ~ estimulo, family = poisson, 
    data = mob_data, prior = priors, iter = iter, chains = chains,
    cores = chains, control = list(adapt_delta = 0.99, max_treedepth = 15),
    file = "./data/raw/individuos.RDS", file_refit = "always")


mod_mor <- brm(formula = MORFOESPECIE ~ estimulo, family = poisson, 
    data = mob_data, prior = priors, iter = iter, chains = chains,
    cores = chains, control = list(adapt_delta = 0.99, max_treedepth = 15),
    file = "./data/raw/morpfoespecie.RDS", file_refit = "always")

Code

mod <- readRDS("./data/raw/individuos.RDS")

extended_summary(fit = mod,
    n.posterior = 1000, fill = viridis(10)[7], trace.palette = viridis,
    highlight = TRUE, remove.intercepts = TRUE, fit.name = "regresion")

5.2 regresion

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 6) Intercept-student_t(3, -2.3, 2.5)	INDIVIDUOS ~ estimulo	5000	4	1	2500	0	0	2879.695	3393.804	211051743

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_estimuloAlopátrico	-1.438	-2.220	-0.719	1	5074.592	4384.992
b_estimuloAlopátricoP	-1.234	-2.032	-0.524	1	5573.468	4803.804
b_estimuloControl	-7.696	-15.037	-3.220	1.001	2879.695	3393.804

5.2.1 Contrastes

Code

cntrs <- as.data.frame(emmeans(mod, pairwise ~ estimulo)$contrasts)
names(cntrs)[3:4] <- c("l-95% CI", "u-95% CI")

coef_table <- brmsish:::html_format_coef_table(cntrs, highlight = TRUE, fill = viridis(10)[7])

print(coef_table)

	contrast	estimate	l-95% CI	u-95% CI
1	Simpátrico - Alopátrico	1.425	0.709	2.202
2	Simpátrico - (Alopátrico +)	1.223	0.499	1.996
3	Simpátrico - Control	7.143	2.655	13.867
4	Alopátrico - (Alopátrico +)	-0.198	-1.136	0.739
5	Alopátrico - Control	5.704	1.257	12.604
6	(Alopátrico +) - Control	5.899	1.356	12.721

5.3 Numero de morfoespecies

Code

mod <- readRDS("./data/raw/morpfoespecie.RDS")

extended_summary(fit = mod,
    n.posterior = 1000, fill = viridis(10)[7], trace.palette = viridis,
    highlight = TRUE, remove.intercepts = TRUE, fit.name = "regresion")

5.4 regresion

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 6) Intercept-student_t(3, -2.3, 2.5)	MORFOESPECIE ~ estimulo	5000	4	1	2500	0	0	3195.343	3277.219	1132684083

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_estimuloAlopátrico	-1.110	-2.044	-0.227	1.001	5091.331	5214.812
b_estimuloAlopátricoP	-0.759	-1.609	0.059	1.001	5157.734	4789.173
b_estimuloControl	-7.011	-13.895	-2.599	1.001	3195.343	3277.219

5.4.1 Contrastes

Code

cntrs <- as.data.frame(emmeans(mod, pairwise ~ estimulo)$contrasts)
names(cntrs)[3:4] <- c("l-95% CI", "u-95% CI")

coef_table <- brmsish:::html_format_coef_table(cntrs, highlight = TRUE, fill = viridis(10)[7])

print(coef_table)

	contrast	estimate	l-95% CI	u-95% CI
1	Simpátrico - Alopátrico	1.097	0.178	1.977
2	Simpátrico - (Alopátrico +)	0.752	-0.078	1.581
3	Simpátrico - Control	6.511	2.205	13.140
4	Alopátrico - (Alopátrico +)	-0.350	-1.378	0.659
5	Alopátrico - Control	5.425	0.844	11.997
6	(Alopátrico +) - Control	5.752	1.302	12.355

Información de la sesión

R version 4.3.1 (2023-06-16)
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] emmeans_1.8.9       brmsish_1.0.0       brms_2.20.4        
 [4] Rcpp_1.0.12         viridis_0.6.4       viridisLite_0.4.2  
 [7] kableExtra_1.3.4    ggplot2_3.4.4       httpuv_1.6.12      
[10] googlesheets4_1.1.1

loaded via a namespace (and not attached):
  [1] pbapply_1.7-2        gridExtra_2.3        remotes_2.4.2.1     
  [4] inline_0.3.19        rlang_1.1.3          magrittr_2.0.3      
  [7] matrixStats_1.1.0    compiler_4.3.1       loo_2.6.0           
 [10] systemfonts_1.0.5    vctrs_0.6.5          reshape2_1.4.4      
 [13] rvest_1.0.3          stringr_1.5.1        pkgconfig_2.0.3     
 [16] crayon_1.5.2         fastmap_1.1.1        backports_1.4.1     
 [19] ellipsis_0.3.2       labeling_0.4.3       utf8_1.2.4          
 [22] threejs_0.3.3        promises_1.2.1       rmarkdown_2.25      
 [25] markdown_1.11        purrr_1.0.2          xfun_0.41           
 [28] jsonlite_1.8.8       later_1.3.1          parallel_4.3.1      
 [31] R6_2.5.1             dygraphs_1.1.1.6     StanHeaders_2.26.28 
 [34] stringi_1.8.3        estimability_1.4.1   cellranger_1.1.0    
 [37] rstan_2.32.3         knitr_1.45           zoo_1.8-12          
 [40] base64enc_0.1-3      bayesplot_1.10.0     Matrix_1.6-2        
 [43] igraph_1.6.0         tidyselect_1.2.0     rstudioapi_0.15.0   
 [46] abind_1.4-5          yaml_2.3.8           codetools_0.2-19    
 [49] miniUI_0.1.1.1       curl_5.1.0           pkgbuild_1.4.3      
 [52] lattice_0.22-5       tibble_3.2.1         plyr_1.8.9          
 [55] shiny_1.7.5.1        withr_3.0.0          bridgesampling_1.1-2
 [58] askpass_1.2.0        posterior_1.5.0      coda_0.19-4         
 [61] evaluate_0.23        sketchy_1.0.2        RcppParallel_5.1.7  
 [64] ggdist_3.3.0         xts_0.13.1           xml2_1.3.5          
 [67] pillar_1.9.0         tensorA_0.36.2       packrat_0.9.2       
 [70] checkmate_2.3.1      DT_0.30              stats4_4.3.1        
 [73] shinyjs_2.1.0        distributional_0.3.2 generics_0.1.3      
 [76] rstantools_2.3.1.1   munsell_0.5.0        scales_1.3.0        
 [79] gtools_3.9.4         xtable_1.8-4         glue_1.7.0          
 [82] tools_4.3.1          xaringanExtra_0.7.0  shinystan_2.6.0     
 [85] colourpicker_1.3.0   webshot_0.5.5        fs_1.6.3            
 [88] mvtnorm_1.2-3        cowplot_1.1.1        grid_4.3.1          
 [91] ape_5.7-1            QuickJSR_1.0.7       crosstalk_1.2.0     
 [94] colorspace_2.1-0     nlme_3.1-163         googledrive_2.1.1   
 [97] cli_3.6.2            rappdirs_0.3.3       fansi_1.0.6         
[100] gargle_1.5.2         svglite_2.1.2        Brobdingnag_1.2-9   
[103] dplyr_1.1.3          gtable_0.3.4         digest_0.6.34       
[106] htmlwidgets_1.6.2    farver_2.1.1         htmltools_0.5.7     
[109] lifecycle_1.0.4      httr_1.4.7           mime_0.12           
[112] openssl_2.1.1        shinythemes_1.2.0