Mobbing ontogeny

Published

May 4, 2026

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

Cargar paquetes

Code

sketchy::load_packages(c(
  "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 = mako(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 = mako(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 = mako(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 = mako(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

Code

mob_data <- read.csv("./data/raw/data_v2.csv",na.strings = "NA")

mob_data$HORA <- format(strptime(mob_data$HORA, "%H:%M"), "%H:%M")

3 Descripción de los datos

49 puntos de muestreo
63 individuos observados
Un máximo de 43 morfoespecies observadas
41% de muestreos tuvieron respuesta (54% 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 = mako(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 row containing non-finite outside the scale range
(`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 = mako(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 = mako(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 = mako(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 = mako(10)[7], trace.palette = mako,
    highlight = TRUE, remove.intercepts = TRUE, fit.name = "regresion")

5.2 regresion

	formula	family	priors	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	INDIVIDUOS ~ estimulo	poisson (log)	b-normal(0, 6) Intercept-student_t(3, -2.3, 2.5)	5000	4	1	2500	0 (0%)	0	2535.577	2569.933	1968743918

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_estimuloAlopátrico	-1.437	-2.136	-0.799	1.001	5333.999	5052.401
b_estimuloAlopátricoP	-1.465	-2.226	-0.798	1	5183.378	4819.503
b_estimuloControl	-7.942	-15.430	-3.443	1.002	2535.577	2569.933

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 = mako(10)[7])

print(coef_table)

	contrast	estimate	l-95% CI	u-95% CI
1	Simpátrico - Alopátrico	1.429	0.784	2.114
2	Simpátrico - (Alopátrico +)	1.444	0.752	2.170
3	Simpátrico - Control	7.442	3.074	14.377
4	Alopátrico - (Alopátrico +)	0.031	-0.876	0.888
5	Alopátrico - Control	5.999	1.417	12.883
6	(Alopátrico +) - Control	5.967	1.223	12.631

5.3 Numero de morfoespecies

Code

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

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

5.4 regresion

	formula	family	priors	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	MORFOESPECIE ~ estimulo	poisson (log)	b-normal(0, 6) Intercept-student_t(3, -2.3, 2.5)	5000	4	1	2500	0 (0%)	0	2667.026	2601.382	2112136511

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_estimuloAlopátrico	-1.292	-2.141	-0.531	1.001	5217.811	5090.511
b_estimuloAlopátricoP	-1.092	-1.908	-0.330	1.001	5339.278	4577.381
b_estimuloControl	-7.372	-14.587	-2.922	1.001	2667.026	2601.382

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 = mako(10)[7])

print(coef_table)

	contrast	estimate	l-95% CI	u-95% CI
1	Simpátrico - Alopátrico	1.277	0.507	2.110
2	Simpátrico - (Alopátrico +)	1.086	0.302	1.866
3	Simpátrico - Control	6.824	2.443	13.640
4	Alopátrico - (Alopátrico +)	-0.196	-1.228	0.749
5	Alopátrico - Control	5.554	1.024	12.358
6	(Alopátrico +) - Control	5.775	1.132	12.458

Información de la sesión

R version 4.5.2 (2025-10-31)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 22.04.5 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  LAPACK version 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       

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.11.1    brmsish_1.0.0     brms_2.23.0       Rcpp_1.1.1-1.1   
[5] viridis_0.6.5     viridisLite_0.4.3 kableExtra_1.4.0  ggplot2_4.0.3    

loaded via a namespace (and not attached):
 [1] svUnit_1.0.8          tidyselect_1.2.1      dplyr_1.1.4          
 [4] farver_2.1.2          loo_2.9.0             tidybayes_3.0.7      
 [7] S7_0.2.2              fastmap_1.2.0         TH.data_1.1-3        
[10] tensorA_0.36.2.1      digest_0.6.39         estimability_1.5.1   
[13] lifecycle_1.0.5       StanHeaders_2.32.10   survival_3.8-6       
[16] magrittr_2.0.5        posterior_1.6.1       compiler_4.5.2       
[19] rlang_1.2.0           tools_4.5.2           yaml_2.3.12          
[22] knitr_1.51            labeling_0.4.3        bridgesampling_1.2-1 
[25] htmlwidgets_1.6.4     curl_7.1.0            pkgbuild_1.4.8       
[28] plyr_1.8.9            xml2_1.5.2            RColorBrewer_1.1-3   
[31] abind_1.4-8           multcomp_1.4-28       purrr_1.2.0          
[34] withr_3.0.2           stats4_4.5.2          grid_4.5.2           
[37] inline_0.3.21         xtable_1.8-8          scales_1.4.0         
[40] MASS_7.3-65           cli_3.6.6             mvtnorm_1.3-3        
[43] rmarkdown_2.31        crayon_1.5.3          generics_0.1.4       
[46] remotes_2.5.0         otel_0.2.0            RcppParallel_5.1.11-1
[49] rstudioapi_0.17.1     reshape2_1.4.5        pbapply_1.7-4        
[52] ape_5.8-1             rstan_2.32.7          stringr_1.6.0        
[55] splines_4.5.2         bayesplot_1.15.0      parallel_4.5.2       
[58] matrixStats_1.5.0     vctrs_0.7.3           V8_6.0.4             
[61] Matrix_1.7-4          sandwich_3.1-1        jsonlite_2.0.0       
[64] arrayhelpers_1.1-0    systemfonts_1.3.1     ggdist_3.3.3         
[67] packrat_0.9.3         tidyr_1.3.2           glue_1.8.1           
[70] codetools_0.2-20      cowplot_1.2.0         distributional_0.5.0 
[73] xaringanExtra_0.8.0   stringi_1.8.7         gtable_0.3.6         
[76] QuickJSR_1.8.1        tibble_3.3.0          pillar_1.11.1        
[79] htmltools_0.5.9       Brobdingnag_1.2-9     R6_2.6.1             
[82] textshaping_1.0.4     evaluate_1.0.5        lattice_0.22-9       
[85] backports_1.5.1       sketchy_1.0.7         rstantools_2.5.0     
[88] svglite_2.2.2         coda_0.19-4.1         gridExtra_2.3        
[91] nlme_3.1-168          checkmate_2.3.4       xfun_0.57            
[94] zoo_1.8-14            pkgconfig_2.0.3