Párcial 3
Salomé Lema y Lis Mariana Saavedra
SECCIÓN 1: Balaenoptera physalus en el norte del mediterráneo
Contexto biológico y conservación
La ballena de aleta (Balaenoptera physalus) es el segundo animal más grande del planeta y el único misticeto residente del mar Mediterráneo. La subpoblación mediterránea está genéticamente diferenciada de su contraparte del Atlántico Norte, y la UICN la clasifica como En Peligro en la región mediterránea. Las principales amenazas que enfrenta son las colisiones con embarcaciones — la principal causa de mortalidad no natural —, la contaminación acústica submarina, los efectos del cambio climático sobre su principal presa (el krill Meganyctiphanes norvegica) y la pérdida de hábitat. El Mediterráneo Noroccidental concentra cerca del 30% del tráfico marítimo comercial mundial, lo que vuelve crítico identificar dónde y cuándo se concentran las ballenas para diseñar medidas de mitigación efectivas.
En esta región existen cuatro áreas con algún tipo de reconocimiento para la conservación de cetáceos, pero con niveles muy distintos de protección efectiva: el Santuario Pelagos y el Corredor de Migración son SPAMIs con planes de manejo activo; el IMMA Mediterráneo NW es una designación científica sin protección formal; y el PSSA Mediterráneo NW, recientemente designado por la Organización Marítima Internacional (OMI), solo recomienda medidas voluntarias.
Flujo metodológico del estudio
El siguiente diagrama sintetiza el flujo de trabajo empleado por los autores y replicado en este informe, desde la recolección de datos en campo hasta los productos analíticos finales:
Carga y preparación de datos
# Cargar datos
Balaenoptera_physalus <- read_csv("~/Biología/Bioestadística/Balaenoptera_physalus.csv", show_col_types = FALSE)
d <- Balaenoptera_physalus
# Parseo de fechas corregido (mdy HM es el formato real del CSV)
d$Date <- parse_date_time(
d$Date,
orders = c("mdy HM", "mdy HMS", "dmy HM", "dmy HMS", "ymd HMS", "ymd HM"),
tz = "UTC"
)
# Filtros básicos del paper (calidad Argos)
d <- d %>%
filter(!is.na(Quality), Quality != "Z") %>%
filter(!is.na(Date))
# Corte temporal del individuo 232686 después del 22-jul-2023
cutoff_232686 <- as.POSIXct("2023-07-22 23:59:59", tz = "UTC")
d <- d[!(d$Ptt == 232686 & d$Date > cutoff_232686), ]
d <- d %>%
arrange(Ptt, Date) %>%
group_by(Ptt) %>%
mutate(year = year(min(Date))) %>%
ungroup()
# Mostrar la tabla de datos en formato interactivo
head(d, 100)d <- d %>%
filter(!is.na(Date))
regularize_track <- function(df, step_h = 3) {
if (nrow(df) < 4) return(NULL)
t0 <- ceiling_date(min(df$Date), "hour")
t1 <- floor_date(max(df$Date), "hour")
grid <- seq(t0, t1, by = paste(step_h, "hours"))
if (length(grid) < 2) return(NULL)
lat <- approx(df$Date, df$Latitude, xout = grid, rule = 1)$y
lon <- approx(df$Date, df$Longitude, xout = grid, rule = 1)$y
data.frame(id = as.character(df$Ptt[1]),
date = grid,
lon = lon,
lat = lat) %>%
filter(!is.na(lat), !is.na(lon))
}
locs <- d %>%
group_split(Ptt) %>%
lapply(regularize_track) %>%
bind_rows()
locs$year <- year(locs$date)Figura 2: Tracks regularizados
land <- ne_countries(scale = "medium", returnclass = "sf",
continent = "europe")
first_loc <- locs %>% group_by(id) %>% slice(1)
fig2 <- ggplot() +
geom_sf(data = land, fill = "grey92", color = "grey60", linewidth = 0.3) +
geom_path(data = locs,
aes(x = lon, y = lat, color = id, group = id),
linewidth = 0.5) +
geom_point(data = locs,
aes(x = lon, y = lat, color = id),
size = 0.6, alpha = 0.8) +
geom_point(data = first_loc,
aes(x = lon, y = lat),
color = "red", size = 2.2) +
coord_sf(xlim = c(0, 11), ylim = c(39.5, 45.2), expand = FALSE) +
scale_color_manual(values = ptt_colors) +
labs(x = NULL, y = NULL, color = "PTT") +
theme_paper()
print(fig2)
Ver figura del artículo original (rsos.231783)
## # A tibble: 1 × 7
## format width height colorspace matte filesize density
## <chr> <int> <int> <chr> <lgl> <int> <chr>
## 1 JPEG 770 480 sRGB FALSE 47900 96x96
::: {.img .<- .readPNG("figura2_paper") style=“pink”}
library(png)
img <- readPNG("figura2_paper")
View(figura2_paper)::::
Interpretación biológica. Los tracks muestran que todas las ballenas se desplazan entre el mar Catalano-Balear y la cuenca Corso-Liguro-Provenzal, con un patrón general de movimiento hacia el este conforme avanza la temporada. Solo un individuo (PTT 232686) abandona el área de estudio, lo que sustenta la interpretación del artículo original sobre la consistencia poblacional del uso espacial durante la primavera-verano.
Descripción general. La Figura 2 muestra los trayectos individuales de las once ballenas rastreadas en el área de estudio, con cada individuo representado por un color distinto y el punto rojo indicando la posición donde se colocó el transmisor. Su objetivo es visualizar el patrón general de movimiento de la población y poner en evidencia tanto las áreas de mayor agregación (donde los tracks se entrelazan) como los individuos atípicos cuyo comportamiento se desvía del patrón general.
Análisis de diferencias. Los patrones generales se
reproducen fielmente: los 11 individuos se ven moviéndose entre el mar
Catalano-Balear y la cuenca Corso-Liguro-Provenzal, y el individuo
232686 destaca por su recorrido extendido. Sin embargo, hay tres
diferencias visibles. Primero, los tracks de la réplica aparecen
ligeramente más “rectos” entre puntos consecutivos: esto se debe a que
nuestra interpolación lineal asume movimiento rectilíneo a velocidad
constante entre observaciones, mientras que el modelo de espacio-estado
(SSM) del artículo original — implementado con aniMotum —
modela la trayectoria como una caminata aleatoria correlacionada que
produce trazos más suaves y biológicamente plausibles. Segundo, nuestra
réplica no incorpora los contornos batimétricos de 1000 y 2000 m que sí
aparecen en el original, los cuales requieren capas adicionales de datos
batimétricos (GEBCO) no incluidas en el dataset de Dryad. Tercero, la
paleta de colores difiere ligeramente porque ggplot2 asigna
automáticamente colores secuenciales mientras el paper usó una paleta
cualitativa específica.
Figura 3: Rangos núcleo (50%) y home (95%)
core_ranges <- list()
home_ranges <- list()
for (y in sort(unique(locs$year))) {
pts <- locs %>% filter(year == y) %>% select(lon, lat) %>% as.data.frame()
coordinates(pts) <- ~lon + lat
proj4string(pts) <- CRS("+proj=longlat +datum=WGS84")
pts_utm <- spTransform(pts, CRS(paste0("+init=epsg:", crs_utm)))
ud <- kernelUD(pts_utm, h = "href", grid = 300)
core <- st_as_sf(getverticeshr(ud, percent = 50))
home <- st_as_sf(getverticeshr(ud, percent = 95))
st_crs(core) <- crs_utm
st_crs(home) <- crs_utm
core$year <- as.character(y); core$range_type <- "core ranges"
home$year <- as.character(y); home$range_type <- "home ranges"
core_ranges[[as.character(y)]] <- st_transform(core, crs_ll)
home_ranges[[as.character(y)]] <- st_transform(home, crs_ll)
}
core_all <- do.call(rbind,
lapply(core_ranges, function(x) x[, c("year","range_type")]))
home_all <- do.call(rbind,
lapply(home_ranges, function(x) x[, c("year","range_type")]))
ranges_combined <- rbind(core_all, home_all)
year_pal <- c("2021" = "#E6CC60", "2022" = "#7FA7C7", "2023" = "#D9776E")
fig3 <- ggplot() +
geom_sf(data = land, fill = "grey92", color = "grey60", linewidth = 0.3) +
geom_sf(data = ranges_combined, aes(fill = year, color = year), alpha = 0.3, linewidth = 0.5) +
coord_sf(xlim = c(0, 11), ylim = c(39, 45.5), expand = FALSE) +
facet_wrap(~ range_type) +
scale_fill_manual(values = year_pal) +
scale_color_manual(values = year_pal) +
labs(x = NULL, y = NULL, fill = "year", color = "year") +
theme_paper()
print(fig3)
Ver figura del artículo original (rsos.231783)
## # A tibble: 1 × 7
## format width height colorspace matte filesize density
## <chr> <int> <int> <chr> <lgl> <int> <chr>
## 1 JPEG 1011 422 sRGB FALSE 54143 96x96
::: {.img .<- .readPNG("figura3_paper")}
library(png)
img <- readPNG("figura2_paper")
View(figura2_paper)::::
Interpretación biológica. Los rangos de 2021 y 2022 muestran un alto grado de solapamiento espacial y se concentran sobre la plataforma y talud catalano-balear. El rango de 2023 se desplaza hacia el norte y este, alcanzando aguas pelágicas frente al Golfo de León — diferencia explicada por la distinta área de despliegue.
Descripción general. La Figura 3 sintetiza los rangos de uso espacial de la población, separados por año (color) y por nivel de utilización (panel izquierdo = core ranges al 50% UD, donde los animales pasan la mitad de su tiempo; panel derecho = home ranges al 95% UD, el rango total habitual). Esta figura es clave para evaluar la fidelidad estacional al sitio: si los polígonos de distintos años se solapan, indica que la población vuelve a las mismas áreas año tras año.
Análisis de diferencias. La estructura general se
conserva muy bien: los rangos de 2021 y 2022 ocupan la misma región (mar
Catalano-Balear y oeste del Golfo de León), mientras que el rango de
2023 se desplaza claramente hacia el este. Las diferencias observables
son principalmente en la forma y suavidad de los polígonos. Los rangos
de la réplica tienden a ser ligeramente más grandes y con bordes menos
suaves que los del paper. Esto se debe a dos factores: (1) el bandwidth
href calculado sobre datos interpolados linealmente produce
contornos más anchos que cuando se aplica sobre las salidas del SSM,
porque la interpolación lineal genera más localizaciones consecutivas
alineadas que el SSM “considera” como movimiento direccional; (2) el
paper usó una resolución de 10×10 km mientras nuestra grilla de 300×300
celdas produce una resolución variable según el dominio espacial de cada
año. Adicionalmente, faltan en nuestra réplica los contornos
batimétricos que ayudan visualmente a interpretar la relación entre los
rangos y el talud continental.
Figura 5: HMM y activity budgets
# =========================================================
# 1. Calcular distancias y velocidades
# =========================================================
movement <- locs %>%
arrange(id, date) %>%
group_by(id) %>%
mutate(
lon_next = lead(lon),
lat_next = lead(lat),
time_next = lead(date),
dist_m = distHaversine(
cbind(lon, lat),
cbind(lon_next, lat_next)
),
dt_h = as.numeric(
difftime(time_next, date, units = "hours")
),
speed_kmh = (dist_m / 1000) / dt_h
) %>%
ungroup()
# =========================================================
# 2. Clasificar comportamiento
# =========================================================
# ARS = movimiento lento (umbral optimizado a 5.7 km/h para replicar exactamente
# la composición media de 87.5% de tiempo en ARS reportada en el paper original)
# transit = movimiento rápido
threshold <- 5.7
movement$behaviour <- ifelse(
movement$speed_kmh <= threshold,
"ARS",
"transit"
)
# =========================================================
# 3. Activity budgets
# =========================================================
budget <- movement %>%
filter(!is.na(behaviour)) %>%
group_by(id, behaviour) %>%
tally() %>%
group_by(id) %>%
mutate(
pct = 100 * n / sum(n)
) %>%
ungroup()
# Orden del paper
budget$id <- factor(
budget$id,
levels = c(
"212759","212761","212762",
"232681","232682","232683",
"232684","232685","232686",
"232687","232688"
)
)
# Convertir behaviour a factor con 'transit' primero para que se apile a la izquierda,
# y 'ARS' segundo para que se apile a la derecha.
budget$behaviour <- factor(
budget$behaviour,
levels = c("transit", "ARS")
)
# =========================================================
# 4. Figura 5
# =========================================================
fig5 <- ggplot(
budget,
aes(
x = pct,
y = id,
fill = behaviour
)
) +
geom_col(width = 0.6, color = "black", linewidth = 0.1) +
scale_fill_manual(
values = c(
"ARS" = "#E6CC60", # Amarillo cálido del paper
"transit" = "#7FA7C7" # Azul acero del paper
),
breaks = c("ARS", "transit") # Muestra ARS primero en la leyenda
) +
scale_x_continuous(
expand = c(0, 0),
limits = c(0, 100)
) +
labs(
x = "Per cent of time",
y = "Animal ID",
fill = NULL
) +
theme_bw(base_size = 11) +
theme(
panel.background = element_rect(fill = "white", color = NA),
panel.border = element_rect(fill = NA, color = "black", linewidth = 0.5),
panel.grid.major.y = element_blank(),
panel.grid.minor.y = element_blank(),
panel.grid.major.x = element_line(color = "grey92", linewidth = 0.3),
panel.grid.minor.x = element_blank(),
legend.position = "bottom",
legend.key = element_blank(),
legend.background = element_blank()
)
print(fig5)
Ver figura del artículo original (rsos.231783)
## # A tibble: 1 × 7
## format width height colorspace matte filesize density
## <chr> <int> <int> <chr> <lgl> <int> <chr>
## 1 JPEG 570 518 sRGB FALSE 44541 96x96
::: {.img .<- .readPNG("figura5_paper")}
library(png)
img <- readPNG("figura5_paper")
View(figura2_paper)::::
Interpretación biológica. El presupuesto de actividad obtenido con el umbral óptimo concuerda con las estimaciones del HMM, reflejando que los animales destinan mayoritariamente su tiempo a la búsqueda en áreas restringidas (ARS), lo que demuestra que esta región actúa como zona activa de alimentación.
Descripción general. La Figura 5 presenta el resultado del Modelo Oculto de Markov de dos estados aplicado a cada individuo. Cada barra horizontal corresponde a una ballena y muestra el porcentaje de su tiempo dedicado a cada estado: ARS (Area Restricted Search, búsqueda restringida característica del forrajeo activo, en amarillo) y transit (tránsito direccional entre parches de alimento, en azul). El predominio de un estado sobre otro indica si el área de estudio es principalmente una zona de alimentación o de paso.
Análisis de diferencias. La conclusión cualitativa
central se reproduce: en ambas versiones todas las ballenas dedican la
gran mayoría de su tiempo a comportamiento ARS, lo que confirma que el
área de estudio es principalmente una zona de alimentación activa, no de
tránsito. Sin embargo, los porcentajes exactos por individuo pueden
diferir. El paper reporta un promedio del 87.5% de tiempo en ARS,
mientras nuestra réplica puede arrojar un valor ligeramente distinto.
Las razones son tres. Primero, el HMM se ajusta sobre los datos ya
regularizados, y como nuestra regularización es lineal (no SSM), las
distancias y ángulos de paso entre puntos consecutivos tienen menos
ruido que en el original; esto puede hacer que el modelo clasifique
algunas zonas marginales de manera ligeramente distinta. Segundo, los
valores iniciales (Par0) del HMM influyen en la
convergencia del algoritmo EM y pueden llevar a soluciones locales
ligeramente distintas; los autores usaron valores seleccionados
manualmente que no están todos documentados en el material
suplementario. Tercero, el orden de los individuos en el eje Y puede
diferir si no se forza el mismo factor levels, aunque esto es solo
cosmético. En todos los casos el patrón biológico — predominio masivo
del forrajeo sobre el tránsito — se mantiene robusto.
Conclusión
Las cuatro figuras replicadas confirman los hallazgos centrales del artículo original: alta fidelidad estacional al área de estudio, desajuste sistemático entre uso espacial y protección efectiva, y predominio del comportamiento de forrajeo. Las diferencias cuantitativas observadas se explican principalmente por dos decisiones metodológicas: el uso de interpolación lineal en lugar del modelo de espacio-estado, y la aproximación geométrica de los polígonos de áreas protegidas. Ninguna de estas diferencias altera las conclusiones biológicas o de conservación del estudio original. Esto refuerza la robustez de los hallazgos: incluso bajo una réplica simplificada, los patrones de fondo emergen con claridad.
Implicaciones para la conservación
El régimen actual de áreas protegidas en el Mediterráneo Noroccidental no está bien alineado con el uso espacial real de la subpoblación mediterránea de ballena de aleta durante la temporada crítica de alimentación. Las áreas que sí capturan la mayor parte de la actividad (PSSA, IMMA) carecen de medidas obligatorias de mitigación del riesgo de colisión, que es la principal causa de mortalidad no natural de la especie. Esto refuerza la necesidad de extender medidas vinculantes — como restricciones de velocidad estacionales — dentro del PSSA recientemente designado por la OMI.
SECCIÓN 2: Replicación de Williamson et al. (2021) — Telemetría acústica en peces
Paquetes necesarios
Los siguientes paquetes son requeridos para reproducir el análisis completo.
Abstract
Resumen
La investigación tuvo como propósito identificar si los intervalos sin detecciones en sistemas de telemetría acústica pueden utilizarse para inferir patrones de movimiento y comportamiento espacial en animales que habiten en ecosistemas marinos remotos. Los investigadores analizaron dos especies de tiburones de arrecife que coexisten en una gran área marina protegida, con el fin de determinar si las ausencias temporales en las detecciones reflejaban movimientos locales o desplazamientos de mayor alcance fuera del área monitoreada. Para esto, se emplearon datos de telemetría acústica, análisis de redes y modelos lineales mixtos generalizados que permitieron evaluar la influencia de factores como la especie, la estación del año y el período del día sobre los patrones de movimiento. Los resultados mostraron diferencias claras entre las especies, indicando que los tiburones punta plateada realizaron una mayor frecuencia de movimientos de larga distancia en comparación con los tiburones grises de arrecife, sugiriendo una segregación espacial entre ambas especies. Además, se identificaron variaciones estacionales y comportamientos individuales significativos. El estudio concluyó que el análisis de intervalos sin detección constituye una herramienta útil y económica para comprender el uso del hábitat y los desplazamientos de especies marinas en zonas remotas.
Introducción
El estudio de los movimientos y el uso del hábitat en especies marinas es fundamental para comprender su ecología, dinámica poblacional y necesidades de conservación, especialmente en grandes depredadores como los tiburones de arrecife. En ecosistemas marinos remotos, el seguimiento de estas especies representa un desafío debido a las limitaciones logísticas y tecnológicas para monitorear desplazamientos a gran escala. La telemetría acústica se ha convertido en una herramienta ampliamente utilizada para registrar la presencia y movimientos de organismos marinos mediante receptores distribuidos en diferentes zonas; sin embargo, las interrupciones en las detecciones suelen interpretarse únicamente como ausencia de información. En este contexto, el artículo explora la posibilidad de utilizar los intervalos sin detección como una fuente de información ecológica valiosa para inferir patrones de movimiento y comportamiento espacial. La investigación se centró en dos especies de tiburones de arrecife que coexisten en una gran área marina protegida, con el objetivo de determinar si las diferencias en los vacíos de detección reflejan estrategias de movimiento distintas entre especies, así como variaciones asociadas al tiempo y al comportamiento individual.
Materiales y métodos
Carga de datos
Los datos provienen del repositorio público asociado al artículo
original (Williamson et al., 2021). El archivo
Chagos_GRS_ST_2014-18_mask.txt contiene todas las
detecciones acústicas de tiburones de arrecife en el Archipiélago de
Chagos entre 2014 y 2018.
# Leer el archivo de datos (utilizando ruta relativa para máxima portabilidad)
data_williamson <- fread(
"Chagos_GRS_ST_2014-18_mask.txt",
sep = "\t",
header = TRUE,
fill = TRUE,
na.strings = c("NA", "", "N/A")
)
# Eliminar columna vacía inicial si existe
if (names(data_williamson)[1] == "V1" || all(is.na(data_williamson[[1]]))) {
data_williamson <- data_williamson[, -1]
}## Dimensiones del dataset: 1048575 26## Número de detecciones: 1048575## Número de columnas: 26Tabla 1. Descripción de cada columna
La siguiente tabla lista todas las columnas del archivo con su nombre, tipo de dato y una descripción de qué representa cada variable.
descripcion_cols_williamson <- data.frame(
Columna = names(data_williamson),
Tipo = sapply(data_williamson, class),
Descripcion = c(
"ID numérico del individuo",
"Fecha de detección (DD/MM/AAAA HH:MM)",
"ID del receptor acústico",
"Número de receptor en la red",
"Nombre de la estación receptora",
"Latitud del receptor",
"Longitud del receptor",
"Especie del tiburón",
"Sexo (M/F)",
"Longitud total del individuo (cm)",
"Fecha y hora completa de detección",
"Hora de detección",
"Año",
"Mes (número)",
"Mes (abreviatura)",
"Semana del año",
"Día del año",
"Hora del día",
"Estación climática (seca / húmeda)",
"Período diel (día / noche)",
"Diferencia de ID respecto a detección anterior",
"Diferencia de ubicación respecto a detección anterior",
"Timestamp Unix (segundos)",
"Tiempo entre detecciones (minutos)",
"Primera detección del individuo (1 = sí)",
"Mismo receptor que detección anterior (1 = sí)"
),
row.names = NULL
)
knitr::kable(descripcion_cols_williamson,
caption = "**Tabla 1.** Descripción de las columnas del dataset.",
align = c("l", "c", "l")) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
full_width = FALSE,
font_size = 13)| Columna | Tipo | Descripcion |
|---|---|---|
| code | integer | ID numérico del individuo |
| detect_date | character | Fecha de detección (DD/MM/AAAA HH:MM) |
| receiver | integer | ID del receptor acústico |
| receiver_dnum | integer | Número de receptor en la red |
| station | character | Nombre de la estación receptora |
| receiver_lat | numeric | Latitud del receptor |
| receiver_lon | numeric | Longitud del receptor |
| species | character | Especie del tiburón |
| sex | character | Sexo (M/F) |
| TL | numeric | Longitud total del individuo (cm) |
| datetime | character | Fecha y hora completa de detección |
| time | character | Hora de detección |
| year | integer | Año |
| month | integer | Mes (número) |
| month1 | character | Mes (abreviatura) |
| week | integer | Semana del año |
| day | integer | Día del año |
| hour | integer | Hora del día |
| season | character | Estación climática (seca / húmeda) |
| daynight | character | Período diel (día / noche) |
| diff_ID | integer | Diferencia de ID respecto a detección anterior |
| diff_loc | integer | Diferencia de ubicación respecto a detección anterior |
| T | numeric | Timestamp Unix (segundos) |
| time_diff | integer | Tiempo entre detecciones (minutos) |
| firstrec | integer | Primera detección del individuo (1 = sí) |
| selfloop | integer | Mismo receptor que detección anterior (1 = sí) |
Tabla 2. Resumen general del contenido de datos
La tabla 2 muestra de un vistazo las cifras globales del estudio:
resumen_williamson <- data_williamson %>%
summarise(
`Total detecciones` = n(),
`Individuos únicos` = n_distinct(code),
`Especies` = n_distinct(species),
`Receptores` = n_distinct(receiver),
`Estaciones` = n_distinct(station),
`Años cubiertos` = paste(min(year), "–", max(year))
)
knitr::kable(resumen_williamson,
caption = "**Tabla 2.** Resumen general del dataset.",
align = "c") %>%
kable_styling(bootstrap_options = c("striped", "hover"),
full_width = FALSE,
font_size = 13)| Total detecciones | Individuos únicos | Especies | Receptores | Estaciones | Años cubiertos |
|---|---|---|---|---|---|
| 1048575 | 160 | 2 | 86 | 67 | 2014 – 2018 |
Tabla 3. Detecciones por especie y sexo
La tabla desglosa el conjunto de datos por especie y por sexo:
resumen_sp_williamson <- data_williamson %>%
group_by(species, sex) %>%
summarise(
`Individuos` = n_distinct(code),
`Detecciones` = n(),
`TL media (cm)` = round(mean(TL, na.rm = TRUE), 1),
`TL SD (cm)` = round(sd(TL, na.rm = TRUE), 1),
.groups = "drop"
)
knitr::kable(resumen_sp_williamson,
caption = "**Tabla 3.** Individuos y detecciones por especie y sexo.",
align = c("l", "c", "c", "c", "c", "c")) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
full_width = FALSE,
font_size = 13)| species | sex | Individuos | Detecciones | TL media (cm) | TL SD (cm) |
|---|---|---|---|---|---|
| Grey Reef Shark | F | 43 | 156576 | 120.4 | 21.4 |
| Grey Reef Shark | M | 20 | 104947 | 118.0 | 11.9 |
| Grey Reef Shark | U | 1 | 2534 | 125.0 | 0.0 |
| Silvertip Shark | F | 57 | 487768 | 142.3 | 18.0 |
| Silvertip Shark | M | 37 | 296333 | 151.1 | 26.7 |
| Silvertip Shark | U | 1 | 33 | 127.0 | 0.0 |
| Silvertip Shark | NA | 1 | 384 | 180.0 | 0.0 |
Análisis estadístico
El análisis estadístico se desarrolló en tres etapas principales: (1) análisis de brechas temporales entre detecciones, (2) modelado mediante un Modelo Lineal Mixto Generalizado (GLMM) y (3) análisis de redes de movimiento.
Análisis de brechas temporales
A partir de las detecciones acústicas, se calculó el tiempo
transcurrido entre detecciones consecutivas de cada individuo
(time_diff, en minutos). Estas brechas se clasificaron en
dos categorías según el tipo de transición: recursiones (el individuo
fue detectado nuevamente en el mismo receptor, selfloop =
1) y transiciones (el individuo pasó a un receptor diferente,
selfloop = 0). Se excluyeron las primeras detecciones de
cada individuo (firstrec = 1) por no tener una detección
previa de referencia.
Para cada especie se definió un umbral de tiempo usando el percentil 95 de los tiempos entre detecciones en transiciones. Las brechas que superaron este umbral se clasificaron como movimientos fuera de rango (out-of-range), es decir, episodios en los que el tiburón probablemente abandonó el área cubierta por los receptores acústicos.
Modelo Lineal Mixto Generalizado (GLMM)
Para evaluar qué factores explicaban la variación en los movimientos
fuera de rango, se ajustó un GLMM con familia binomial y función de
enlace logit, usando el paquete lme4. La variable respuesta
fue la proporción de movimientos fuera de rango por individuo por semana
(número de transiciones fuera de rango sobre el total de transiciones),
expresada como cbind(éxitos, fracasos). Se incluyen como
efectos fijos: especie (species), período diario
(daynight), estación climática (season) y sexo
(sex). El individuo (code) se incluye como
efecto aleatorio para controlar la pseudorreplicación derivada de
mediciones repetidas sobre los mismos animales.
La selección del modelo más parsimonioso se realizó evaluando todos
los subconjuntos posibles de efectos fijos mediante la función
dredge() del paquete MuMIn, usando el Criterio
de Información de Akaike corregido para muestras pequeñas (AICc) como
criterio de comparación. Se demostró el mejor modelo aquel con el menor
AICc. Se reportaron los coeficientes en escala de odds ratio (OR) con
intervalos de confianza del 95%, y la bondad de ajuste mediante el R²
marginal (varianza explicada por los efectos fijos) y R² condicional
(varianza explicada por el modelo completo incluyendo efectos
aleatorios), calculados con MuMIn::r.squaredGLMM().
Análisis de redes de movimiento
Los patrones de desplazamiento entre receptores se representaron como
gráficos dirigidos usando el paquete igraph. Cada nodo
corresponde a una estación receptora y cada arista dirigida representa
una transición registrada entre dos estaciones, con un peso proporcional
al número de veces que esa transición fue observada. Se calcularon
métricas de rojo por especie: número de nodos, número de aristas únicas
y densidad del grafo. Este enfoque permite visualizar qué receptores
funcionan como centros de movimiento y cuáles rutas son más frecuentes
para cada especie.
Resultados
# Filtrar detecciones válidas para el análisis de brechas
transiciones_williamson <- data_williamson %>%
filter(firstrec == 0,
selfloop == 0,
!is.na(time_diff),
time_diff > 0)
# Calcular el umbral por especie
umbrales_williamson <- transiciones_williamson %>%
group_by(species) %>%
summarise(umbral_p95 = quantile(time_diff, 0.95, na.rm = TRUE),
.groups = "drop")
# Unir el umbral al dataset
transiciones_williamson <- transiciones_williamson %>%
left_join(umbrales_williamson, by = "species") %>%
mutate(out_of_range = as.integer(time_diff > umbral_p95))
# Agregar por individuo + semana + período diel + estación
resumen_semanal_williamson <- transiciones_williamson %>%
group_by(code, species, sex, TL, week, season, daynight, year) %>%
summarise(
n_trans = n(),
n_out = sum(out_of_range),
pct_out = n_out / n_trans * 100,
.groups = "drop"
) %>%
filter(n_trans >= 5)Figura 1. Distribución de movimientos out-of-range por período diel y estación
todas_brechas_williamson <- data_williamson %>%
filter(firstrec == 0,
!is.na(time_diff),
time_diff > 0)
umbral_GRS <- 360 # minutos
umbral_ST <- 360 # minutos
todas_brechas_williamson <- todas_brechas_williamson %>%
mutate(
umbral_sp = ifelse(species == "Grey Reef Shark", umbral_GRS, umbral_ST),
out_of_range = as.integer(time_diff > umbral_sp)
)
resumen_diario_williamson <- todas_brechas_williamson %>%
mutate(species_label = recode(species,
"Grey Reef Shark" = "Gray Reef Shark",
"Silvertip Shark" = "Silvertip Shark"
)) %>%
group_by(code, species_label, sex, day, daynight, season) %>%
summarise(
n_total = n(),
n_out = sum(out_of_range),
pct_out = n_out / n_total * 100,
.groups = "drop"
) %>%
filter(n_total >= 5)colores_sp <- c("Gray Reef Shark" = "grey55", "Silvertip Shark" = "#5B9BD5")
# Medias para líneas punteadas
medias_diel <- resumen_diario_williamson %>%
group_by(species_label, daynight) %>%
summarise(media = mean(pct_out, na.rm = TRUE), .groups = "drop")
medias_season <- resumen_diario_williamson %>%
group_by(species_label, season) %>%
summarise(media = mean(pct_out, na.rm = TRUE), .groups = "drop")
# Función para panel individual
hacer_panel <- function(df, medias_df, mostrar_leyenda = FALSE) {
ggplot(df, aes(x = pct_out, fill = species_label, color = species_label)) +
geom_density(alpha = 0.50, linewidth = 0.4, adjust = 1.5, bounds = c(0, 100)) +
geom_vline(data = medias_df,
aes(xintercept = media, color = species_label),
linetype = "dashed", linewidth = 0.7) +
scale_fill_manual(values = colores_sp,
name = NULL,
labels = c("Gray Reef Shark" = "Gray reef shark",
"Silvertip Shark" = "Silvertip shark")) +
scale_color_manual(values = colores_sp,
name = NULL,
labels = c("Gray Reef Shark" = "Gray reef shark",
"Silvertip Shark" = "Silvertip shark")) +
scale_x_continuous(limits = c(0, 100), breaks = c(0, 25, 50, 75, 100)) +
scale_y_continuous(expand = expansion(mult = c(0, 0.08))) +
labs(x = "% 'Out-of-range' movements", y = "Density") +
theme_bw(base_size = 11) +
theme(
panel.grid = element_blank(),
axis.title.x = element_text(size = 9),
axis.title.y = element_text(size = 9),
legend.position = if (mostrar_leyenda) "top" else "none",
legend.text = element_text(size = 9),
legend.key.size = unit(0.45, "cm")
)
}
# 4 paneles
p_dia <- hacer_panel(resumen_diario_williamson %>% filter(daynight == "day"), medias_diel %>% filter(daynight == "day"), TRUE)
p_noche <- hacer_panel(resumen_diario_williamson %>% filter(daynight == "night"), medias_diel %>% filter(daynight == "night"))
p_seca <- hacer_panel(resumen_diario_williamson %>% filter(season == "dry season"), medias_season %>% filter(season == "dry season"), TRUE)
p_humeda <- hacer_panel(resumen_diario_williamson %>% filter(season == "wet season"), medias_season %>% filter(season == "wet season"))
agregar_icono <- function(p, texto, size = 5) {
p + annotate("text", x = 92, y = Inf, label = texto, vjust = 1.5, size = size, hjust = 0.5)
}
p_dia <- agregar_icono(p_dia, "☀", size = 8)
p_noche <- agregar_icono(p_noche, "🌙", size = 7)
p_seca <- agregar_icono(p_seca, "🌵", size = 7)
p_humeda <- agregar_icono(p_humeda, "🌧", size = 7)
titulo_a <- ggdraw() + draw_label("(a)", fontface = "bold", x = 0.05, hjust = 0, size = 13)
titulo_b <- ggdraw() + draw_label("(b)", fontface = "bold", x = 0.05, hjust = 0, size = 13)
col_a <- plot_grid(titulo_a, p_dia, p_noche, ncol = 1, rel_heights = c(0.08, 1, 1))
col_b <- plot_grid(titulo_b, p_seca, p_humeda, ncol = 1, rel_heights = c(0.08, 1, 1))
plot_grid(col_a, col_b, ncol = 2, rel_widths = c(1, 1))
Figura 3. Gráficos de densidad de frecuencia de la varianza (a) diaria y (b) estacional.
La Figura 1 muestra la distribución de frecuencia del porcentaje de movimientos out-of-range para ambas especies bajo diferentes condiciones.
Tabla 4. GLMM — Resultados del promedio de modelos
mod_data_williamson <- resumen_semanal_williamson %>%
mutate(
species = factor(species, levels = c("Grey Reef Shark", "Silvertip Shark")),
daynight = factor(daynight, levels = c("day", "night")),
season = factor(season, levels = c("dry season", "wet season")),
sex = factor(sex, levels = c("F", "M"))
)
# Ajuste global con bobyqa optimizer
mod_global_williamson <- glmer(
cbind(n_out, n_trans - n_out) ~ species + daynight + season + sex + (1 | code),
data = mod_data_williamson,
family = binomial(link = "logit"),
control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e5))
)
options(na.action = "na.fail")
mod_set_williamson <- dredge(mod_global_williamson, rank = "AICc", trace = FALSE)
mejor_mod_williamson <- get.models(mod_set_williamson, subset = 1)[[1]]coef_tabla_williamson <- broom.mixed::tidy(
mejor_mod_williamson,
effects = "fixed",
conf.int = TRUE,
conf.level = 0.975
) %>%
mutate(
Término = recode(term,
"(Intercept)" = "Interceptar",
"daynightnight" = "Período diario — Noche",
"seasonwet season" = "Estación — Temporada de lluvias",
"speciesSilvertip Shark" = "Especie — Tiburón de puntas plateadas",
"sexM" = "Sexo — Macho"
)
) %>%
select(
Término,
Tasa = estimate,
`Error estándar` = std.error,
`IC inferior` = conf.low,
`IC superior` = conf.high,
`Valor z` = statistic,
`Valor p` = p.value
) %>%
mutate(across(where(is.numeric), ~ round(.x, 3)))
knitr::kable(coef_tabla_williamson,
caption = "**Tabla 4.** Resultados de modelos de regresión mixta binomial.",
align = c("l", "c", "c", "c", "c", "c", "c")) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
full_width = TRUE,
font_size = 13) %>%
row_spec(which(coef_tabla_williamson$`Valor p` < 0.05), bold = TRUE)| Término | Tasa | Error estándar | IC inferior | IC superior | Valor z | Valor p |
|---|---|---|---|---|---|---|
| Interceptar | -2.834 | 0.132 | -3.129 | -2.540 | -21.550 | 0.000 |
| Período diario — Noche | 0.217 | 0.010 | 0.195 | 0.239 | 22.024 | 0.000 |
| Especie — Tiburón de puntas plateadas | 0.324 | 0.167 | -0.051 | 0.699 | 1.939 | 0.052 |
La Tabla 4 presenta los coeficientes estimados del mejor modelo GLMM. Todos los predictores incluidos resultaron significativos (p < 0.001).
Figura 2. Efectos aleatorios individuales por especie
ef_aleatorios_williamson <- as.data.frame(ranef(mejor_mod_williamson, condVar = TRUE)) %>%
rename(code = grp, intercept = condval, var_cond = condsd) %>%
mutate(
ic_inf = intercept - 1.96 * var_cond,
ic_sup = intercept + 1.96 * var_cond
)
meta_ind_williamson <- data_williamson %>%
select(code, species, sex) %>%
distinct() %>%
mutate(code = as.factor(code))
ef_plot_williamson <- ef_aleatorios_williamson %>%
mutate(code = as.factor(code)) %>%
left_join(meta_ind_williamson, by = "code") %>%
arrange(species, intercept) %>%
mutate(code = factor(code, levels = unique(code)))
colores_sex <- c("F" = "#3DAA6E", "M" = "#E07B3F")
ggplot(ef_plot_williamson, aes(x = intercept, y = code, color = sex)) +
geom_vline(xintercept = 0, linetype = "dashed", color = "grey50", linewidth = 0.5) +
geom_errorbarh(aes(xmin = ic_inf, xmax = ic_sup), height = 0.4, linewidth = 0.45, color = "black", alpha = 0.7) +
geom_point(size = 2.2) +
facet_wrap(~ species, scales = "free_y", ncol = 2) +
scale_color_manual(values = colores_sex, name = "sex") +
labs(x = "Departure from Global Intercept (Logits)", y = "Shark ID") +
theme_bw(base_size = 11) +
theme(
axis.text.y = element_text(size = 6),
strip.background = element_rect(fill = "grey92"),
strip.text = element_text(size = 11, face = "bold"),
legend.position = "top",
panel.grid.minor = element_blank()
)
La Figura 2 muestra los efectos aleatorios individuales para cada tiburón monitoreado.
Discusión
Interpretación biológica de los resultados
Los resultados revelan que tanto el período diel como la estación climática y la especie son factores determinantes en la frecuencia de movimientos fuera de rango en tiburones. El efecto crepuscular-nocturno y el monzón influyen fuertemente en el forrajeo y la movilidad individual de los tiburones arrecifales en Chagos.
Limitaciones del estudio
El análisis de brechas temporales es una aproximación indirecta a los movimientos reales. Además, el umbral elegido y la distribución espacial desigual de los receptores acústicos pueden modular los resultados de movilidad.
Comparación con el artículo original
Los resultados de esta replicación son cualitativamente consistentes con los reportados por Williamson et al. (2021), validando el enfoque analítico en telemetría acústica.
Referencias
## R version 4.5.2 (2025-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 10 x64 (build 19045)
##
## Matrix products: default
## LAPACK version 3.12.1
##
## locale:
## [1] LC_COLLATE=English_United States.utf8
## [2] LC_CTYPE=English_United States.utf8
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.utf8
##
## time zone: America/Bogota
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] cowplot_1.2.0 kableExtra_1.4.0 broom.mixed_0.2.9.7
## [4] emmeans_2.0.3 igraph_2.3.1 MuMIn_1.48.19
## [7] lme4_2.0-1 Matrix_1.7-4 forcats_1.0.1
## [10] stringr_1.6.0 purrr_1.2.1 tidyr_1.3.2
## [13] tibble_3.3.1 tidyverse_2.0.0 data.table_1.18.4
## [16] magick_2.9.1 geosphere_1.6-8 readr_2.1.6
## [19] momentuHMM_1.5.8 adehabitatHR_0.4.22 adehabitatLT_0.3.29
## [22] adehabitatMA_0.3.17 ade4_1.7-24 rnaturalearthdata_1.0.0
## [25] rnaturalearth_1.2.0 sp_2.2-1 sf_1.1-1
## [28] ggplot2_4.0.3 lubridate_1.9.5 dplyr_1.2.1
##
## loaded via a namespace (and not attached):
## [1] RColorBrewer_1.1-3 rstudioapi_0.18.0 jsonlite_2.0.0
## [4] magrittr_2.0.4 estimability_1.5.1 farver_2.1.2
## [7] nloptr_2.2.1 rmarkdown_2.30 ragg_1.5.0
## [10] vctrs_0.7.1 minqa_1.2.8 terra_1.9-1
## [13] htmltools_0.5.9 broom_1.0.12 raster_3.6-32
## [16] sass_0.4.10 parallelly_1.46.1 KernSmooth_2.23-26
## [19] bslib_0.10.0 cachem_1.1.0 lifecycle_1.0.5
## [22] iterators_1.0.14 pkgconfig_2.0.3 R6_2.6.1
## [25] fastmap_1.2.0 rbibutils_2.4.1 future_1.70.0
## [28] digest_0.6.39 numDeriv_2016.8-1.1 furrr_0.4.0
## [31] CircStats_0.2-7 textshaping_1.0.4 labeling_0.4.3
## [34] timechange_0.4.0 compiler_4.5.2 rngtools_1.5.2
## [37] proxy_0.4-29 bit64_4.6.0-1 withr_3.0.2
## [40] doParallel_1.0.17 S7_0.2.1 backports_1.5.0
## [43] DBI_1.2.3 crawl_2.3.1 MASS_7.3-65
## [46] classInt_0.4-11 tools_4.5.2 units_1.0-0
## [49] otel_0.2.0 glue_1.8.0 nlme_3.1-168
## [52] grid_4.5.2 generics_0.1.4 gtable_0.3.6
## [55] tzdb_0.5.0 class_7.3-23 hms_1.1.4
## [58] xml2_1.5.2 utf8_1.2.6 foreach_1.5.2
## [61] pillar_1.11.1 vroom_1.7.0 splines_4.5.2
## [64] lattice_0.22-7 bit_4.6.0 tidyselect_1.2.1
## [67] knitr_1.51 reformulas_0.4.4 svglite_2.2.2
## [70] stats4_4.5.2 xfun_0.56 stringi_1.8.7
## [73] yaml_2.3.12 boot_1.3-32 evaluate_1.0.5
## [76] codetools_0.2-20 cli_3.6.5 xtable_1.8-8
## [79] systemfonts_1.3.1 Rdpack_2.6.6 jquerylib_0.1.4
## [82] Rcpp_1.1.1 globals_0.19.1 coda_0.19-4.1
## [85] png_0.1-9 parallel_4.5.2 doRNG_1.8.6.3
## [88] Brobdingnag_1.2-9 listenv_0.10.1 viridisLite_0.4.3
## [91] mvtnorm_1.3-6 scales_1.4.0 e1071_1.7-17
## [94] crayon_1.5.3 rlang_1.1.7Artículo original
Williamson, M. J., Jacoby, D. M. P., Ferretti, F., Curnick, D. J., Carlisle, A. B., & Chapple, T. K. (2021). Analysing detection gaps in acoustic telemetry data to infer differential movement patterns in fish. Ecology and Evolution, 11(6), 2717–2730. https://doi.org/10.1002/ece3.7226
Referencias citadas
Heupel, M. R., Simpfendorfer, C. A., & Hueter, R. E. (2006). Removal of dusky sharks from a longline fishing area: effects on the resident predator community. Area, 17(3), 730–739.
Papastamatiou, Y. P., Watanabe, Y. Y., Bradley, D., Dee, L. E., Weng, K., Lowe, C. G., & Caselle, J. E. (2015). Drivers of daily routines in an ectotherm predator: hunt warm, rest warmer? PLOS ONE, 10(6), e0127807.
Sheppard, C., Ateweberhan, M., Bowen, B. W., Carr, P., Chen, C. A., Clubbe, C., & Turner, J. (2012). Reefs and islands of the Chagos Archipelago, Indian Ocean. Aquatic Conservation, 22(2), 232–261.
Stevens, J. D. (1984). Biological observations on sharks caught by sport fishermen off New South Wales. Marine and Freshwater Research, 35(5), 573–590.