Colonias de Pheidole bergi en caminos rurales del desierto del Monte central (Mendoza, Argentina)

---

Resumen

Analizamos cómo varían la abundancia y arreglo espacial de colonias activas de Pheidole bergi con la heterogeneidad del hábitat a distintas escalas espaciales y con los cambios embientales de largo plazo en el desierto del Monte central (Reserva de la Biósfera de Ñacuñán y alrededores, Provincia de Mendoza, Argentina). Las densidades de colonias en grillas permanentes dentro de un algarrobal protegido del pastoreo fueron similares a las encontradas sobre caminos rurales (de tierra y muy bajo tránsito) que lo cruzan, pero mostraron una reducción >50% en el periodo 2001–2025. Hubo variaciones significativas de densidad a lo largo de algunos caminos rurales, en especial en los que cruzan algarrobales protegidos. La densidad de colonias fue 40–50% menor en los que cruzan otros hábitats (dunas y algarrobales de campos con pastoreo). El agregamiento de colonias detectado parece ser respuesta a heterogeneidad ambiental (e.g., características edáficas, de la vegetación o de los recursos asociados) más que a interacciones positivas (dependientes de la distancia) entre las colonias.

---

Introducción

• Preguntas:
  • ¿Qué factores afectan la abundancia y disposición espacial de una hormiga generalista en un ambiente árido?
  • ¿Cómo se modifican con un aumento temporario/permanente del estrés (sequía/aridez/ganadería)?
• Contexto
  • Hormiga omnívora/generalista: Pheidole bergi
  • Monte central
    • Heterogeneidad ambiental a dos escalas: macro (suelo-> agua -> vegetación = “ambientes”) y micro (vegetación leñosa/perenne -> estreś abiótico/recursos/procesos = “parches”)
    • Heterogeneidad temporal: sequías
    • Perturbaciones: heterogeneidad ambiental macro (algarrobal pastoreado, caminos disturbados), aridez (cambio climático)
• Disposición espacial
  • como causa (uso de hábitat, impacto, dinámica poblacional)
  • como efecto (indicador de factores/procesos subyacentes)
• Preguntas específicas:
  • ¿Cómo varian la abundancia y disposicion de colonias de P. bergi de acuerdo a la heterogeneidad espacial (natural, antrópica e interacción)?
  • ¿Esos patrones son permanentes o se modifican con el aumento a largo plazo de la aridez [antecedente local: reducción drástica de abundancia de hormigas granívoras especialistas Pogonomyrmex spp.]?

Métodos

• Reserva de la Biosfera de Ñacuñán (Mendoza, Argentina) y campos ganaderos vecinos
• Relevamiento de colonias de Pb en caminos rurales de bajo tránsito y mantenimiento (picadas)
• Febrero 2019 y febrero 2025: dos segmentos de 1500 m sobre sendas picadas de 7 m de ancho en algarrobal (precisión: 1 dm)
• Febrero 2021: nueve segmentos de 500 m sobre picadas de 6 m de ancho, tres por ambiente: algarrobal, algarrobal pastoreado y medanal (precisión: 1 cm) [Nota: una transecta en pastoreado (P2) son dos segmentos de 231 y 269 m (separados por 100 m de otro ambiente)]
• Para comparar densidades se usan datos de relevamiento de colonias en dos grillas (J, V) de 50×50 m [grano/precisión: 1 m²] en febrero 2001 y 2019 [sus patrones espaciales se analizan aparte]
• R + spatstats

knitr::include_graphics("Fig_picadas.png")

Figura 1: Arriba izquierda: Pheidole bergi (obrera). Arriba derecha: entrada típica de colonia de P.bergi. Abajo izquierda: buscando y mapeando colonias a lo largo de un segmento de camino rural (F) que cruza el algarrobal protegido dentro de la Reserva de Ñacuñán (Mendoza, Argentina). Abajo derecha: un camino rural que cruza los médanos y, más adelante, continúa a través de algarrobales.

Análisis y Resultados

library(spatstat) # documented in "Spatial Point Patterns: Methodology and Applications with R" by Baddeley et al (2016)
library(dplyr) # para resumir datos por transecta
library(ggplot2) # para figuras
library(gridExtra) #para combinar figuras ggplot
library(png) # para importar las imágenes con readPNG
library(scales) # para importar las imágenes con readPNG

tema <- theme_bw(base_size = 16) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "none") #tema base para figuras ggplot

#Base de datos de picadas
Caminos <- read.delim("Colonias_caminos.txt", colClasses = c("factor", "factor", "factor", "numeric"))

#Agrupados por transecta (áreas muestreadas: 2019 y 2025: 1500m x 7 m; 2021: 500m x 6m)
Colonias <- Caminos %>%
              group_by(Year, Transect) %>%
              summarise(Year = unique(Year), Habitat = unique(Habitat), Nests = n()) %>%
              mutate(Density = if_else(Year == "2021", Nests / (500*6) * 10000, Nests / (1500*7) * 10000))
Colonias <- Colonias %>% arrange(Year, Habitat)

#Transectas como redes de un solo segmento para análisis con spatstat
#Las dos transectas medidas en 2019 son de 1500 m de largo
Pb_2019D <- lpp(data.frame(cbind(x = with(Caminos, x[Transect == "D" & Year == "2019"]), y = 0)), linnet(ppp(x = c(0, 1500), y = c(0, 0), window = owin(c(0,1500), c(0,1))), m = matrix(c(FALSE, TRUE, TRUE, FALSE), 2, 2)))
Pb_2019F <- lpp(data.frame(cbind(x = with(Caminos, x[Transect == "F" & Year == "2019"]), y = 0)), linnet(ppp(x = c(0, 1500), y = c(0, 0), window = owin(c(0,1500), c(0,1))), m = matrix(c(FALSE, TRUE, TRUE, FALSE), 2, 2)))
#las mismas en 2025
Pb_2025D <- lpp(data.frame(cbind(x = with(Caminos, x[Transect == "D" & Year == "2025"]), y = 0)), linnet(ppp(x = c(0, 1500), y = c(0, 0), window = owin(c(0,1500), c(0,1))), m = matrix(c(FALSE, TRUE, TRUE, FALSE), 2, 2)))
Pb_2025F <- lpp(data.frame(cbind(x = with(Caminos, x[Transect == "F" & Year == "2025"]), y = 0)), linnet(ppp(x = c(0, 1500), y = c(0, 0), window = owin(c(0,1500), c(0,1))), m = matrix(c(FALSE, TRUE, TRUE, FALSE), 2, 2)))

Pb2019y25 <- list("Algarrobal (D) - 2019" = Pb_2019D,
               "Algarrobal (D) - 2025" = Pb_2025D,
               "Algarrobal (F) - 2019" = Pb_2019F,
               "Algarrobal (F) - 2025" = Pb_2025F)

#Las medidas en 2021 son de 500 m de largo
#Algarrobal
Pb_2021_AlgD3 <- lpp(data.frame(cbind(x = with(Caminos, x[Transect == "D3"]), y = 0)), linnet(ppp(x = c(0, 500), y = c(0, 0), window = owin(c(0,500), c(0,1))), m = matrix(c(FALSE, TRUE, TRUE, FALSE), 2, 2)))
Pb_2021_AlgF1 <- lpp(data.frame(cbind(x = with(Caminos, x[Transect == "F1"]), y = 0)), linnet(ppp(x = c(0, 500), y = c(0, 0), window = owin(c(0,500), c(0,1))), m = matrix(c(FALSE, TRUE, TRUE, FALSE), 2, 2)))
Pb_2021_AlgT <- lpp(data.frame(cbind(x = with(Caminos, x[Transect == "T"]), y = 0)), linnet(ppp(x = c(0, 500), y = c(0, 0), window = owin(c(0,500), c(0,1))), m = matrix(c(FALSE, TRUE, TRUE, FALSE), 2, 2)))

#Pastoreado
Pb_2021_PastL <- lpp(data.frame(cbind(x = with(Caminos, x[Transect == "L"]), y = 0)), linnet(ppp(x = c(0, 500), y = c(0, 0), window = owin(c(0,500), c(0,1))), m = matrix(c(FALSE, TRUE, TRUE, FALSE), 2, 2)))
Pb_2021_PastP1 <- lpp(data.frame(cbind(x = with(Caminos, x[Transect == "P1"]), y = 0)), linnet(ppp(x = c(0, 500), y = c(0, 0), window = owin(c(0,500), c(0,1))), m = matrix(c(FALSE, TRUE, TRUE, FALSE), 2, 2)))
Pb_2021_PastP2 <- lpp(data.frame(cbind(x = with(Caminos, x[Transect == "P2"]), y = 0)), linnet(ppp(x = c(0, 500), y = c(0, 0), window = owin(c(0,500), c(0,1))), m = matrix(c(FALSE, TRUE, TRUE, FALSE), 2, 2)))
#[ojo que la transecta P2 es discontinua (en el medio tiene parche de otro ambiente): no analizar espacial]

#Medanal
Pb_2021_MedD1 <- lpp(data.frame(cbind(x = with(Caminos, x[Transect == "D1"]), y = 0)), linnet(ppp(x = c(0, 500), y = c(0, 0), window = owin(c(0,500), c(0,1))), m = matrix(c(FALSE, TRUE, TRUE, FALSE), 2, 2)))
Pb_2021_MedD2 <- lpp(data.frame(cbind(x = with(Caminos, x[Transect == "D2"]), y = 0)), linnet(ppp(x = c(0, 500), y = c(0, 0), window = owin(c(0,500), c(0,1))), m = matrix(c(FALSE, TRUE, TRUE, FALSE), 2, 2)))
Pb_2021_MedM <- lpp(data.frame(cbind(x = with(Caminos, x[Transect == "M"]), y = 0)), linnet(ppp(x = c(0, 500), y = c(0, 0), window = owin(c(0,500), c(0,1))), m = matrix(c(FALSE, TRUE, TRUE, FALSE), 2, 2)))

Pb2021 <- list("Algarrobal (D3)" = Pb_2021_AlgD3,
               "Algarrobal (F1)" = Pb_2021_AlgF1,
               "Algarrobal (T)" = Pb_2021_AlgT,
               "Grazed algarrobal (L)" = Pb_2021_PastL,
               "Grazed algarrobal (P1)" = Pb_2021_PastP1,
               "Grazed algarrobal (P2)" = Pb_2021_PastP2,
               "Sand dunes (D1)" = Pb_2021_MedD1,
               "Sand dunes (D2)" = Pb_2021_MedD2,
               "Sand dunes (M)" = Pb_2021_MedM)

#Base de datos de parcelas (=grillas) para estimar densidades
# Imágenes raster (png) binarias en grillas 50x50 (pixel = 1 m^2^ => 2500 pixeles)
# nidos como objetos spatstat::im (imágenes raster de pixeles regulares; con nido: pixel = 0 | sin nido: pixel = 1)
# spatstat::transmat cambia convención de x e y de la matriz de puntos raster al estándar interno de spatstat
J2001_ants_50 <- im(transmat(readPNG("J2001_nidos_50x50.png"), from = "European", to = "spatstat"), unitname = "m")
J2019_ants_50 <- im(transmat(readPNG("J2019_nidos_50x50.png"), from = "European", to = "spatstat"), unitname = "m")
V2001_ants_50 <- im(transmat(readPNG("V2001_nidos_50x50.png"), from = "European", to = "spatstat"), unitname = "m")
V2019_ants_50 <- im(transmat(readPNG("V2019_nidos_50x50.png"), from = "European", to = "spatstat"), unitname = "m")

#Armar la base de datos de densidades de grillas y de picadas
#Ojo que al tomar el número de pixeles con colonias, en V para ambos años se tiene una colonia menos, ya que en verdad los mapas muestran presencia/ausencia y no abundancia => agregado "a mano" (+1)
Grillas <- tibble(Id = c("J","J", "V", "V"),
                  Year = c("2001", "2019", "2001", "2019"),
                  Habitat = rep("Algarrobal", 4),
                  Nests = c(sum(!J2001_ants_50$v), sum(!J2019_ants_50$v), sum(!V2001_ants_50$v) + 1, sum(!V2019_ants_50$v) + 1),
                  Density = Nests / 2500 * 10000)

Densidades <- bind_rows("Areal (grids)" = Grillas, "Lineal (dirt roads)" = rename(Colonias, Id = Transect), .id = "Sample")
#para agrupar en figura (~ Año * Sample * Habitat pero incompleto)
Densidades$Grupo <- as.factor(c(1, 2, 1, 2, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7, 7))

Densidad de colonias activas

La densidad de colonias de P. bergi en Ñacuñán fue estimada en diferentes ambientes y con distintos muestreadores a lo largo de los años (Fig. 2).

#Densidad de colonias: todas las muestras
ggplot(Densidades, aes(x = interaction(Year, Id), y = Density, fill = factor(Habitat, levels = c("Algarrobal", "Grazed", "Dunes")), colour = Sample)) +
  geom_col(width = .75) +
  labs(y = "Density (nests/ha)", x = "Sample site") +
  scale_x_discrete(limits = c("2001.J", "2001.V", "2019.J", "2019.V", "2019.D", "2019.F", "2021.T", "2021.F1", "2021.D3", "2021.L", "2021.P1", "2021.P2", "2021.D1", "2021.D2", "2021.M", "2025.D", "2025.F"), labels = function(x) {substr(x, 6, nchar(x))}) +
  scale_colour_manual(name = "Sample type", values = c("black", "transparent"), labels = c("Areal (grids)", "Lineal (dirt roads)")) +
  scale_fill_manual(name = "Habitat", values = c("forestgreen", "olivedrab3", "darkkhaki"), labels = c("Algarrobal", "Grazed algarrobal", "Sand dunes")) +
  annotate("text", x = c("2001.J", "2019.V", "2021.P1", "2025.D"), y = c(380, 380, 380, 380), label = c("2001", "2019", "2021", "2025"), hjust = c(0, 0, 0.5, 0)) +
  geom_vline(xintercept = c(2.5, 6.5, 15.5), lty = 2, col = "darkgrey") +
  tema + theme(legend.position = "right")

ggplot(Densidades, aes(x = Grupo, y = Density, colour = factor(Habitat, levels = c("Algarrobal", "Grazed", "Dunes")), shape = Sample)) +
  geom_point(size = 3, stroke = 1.2) +
  geom_line(aes(group = Id), lty = 3, lwd = 1) +
  stat_summary(fun = mean, geom = "point", size = 6, shape = "━") +
  labs(y = "Density (nests/ha)", x = "Sample sites") +
  scale_x_discrete(labels = c("J|V", "J|V", "D|F", "D3|F1|T", "D1|D2|M", "L|P1|P2", "D|F")) +
  scale_shape_manual(name = "Sample type", values = c(15, 21), labels = c( "Areal (grids)", "Lineal (dirt roads)")) +
  scale_colour_manual(name = "Habitat", values = c("forestgreen", "olivedrab3", "darkkhaki"), labels = c("Algarrobal", "Grazed algarrobal", "Sand dunes")) +
  annotate("text", x = c(1, 2.5, 5, 7), y = c(380, 380, 380, 380), label = c("2001", "2019", "2021", "2025"), size = 6) +
  geom_vline(xintercept = c(1.5, 3.5, 6.5), lty = 2, lwd = 0.75, col = "darkgrey") +
  tema + theme(legend.position = "right", axis.text.x = element_text(size = 10))

Figura 2: Arriba: Densidad de colonias de P. bergi en Ñacuñán en diferentes ambientes, estimada con muestreadores areales (parcelas permanentes de 50 × 50 m) y lineales (segmentos de 500–1500 m de caminos rurales) en el periodo 2001–2025. Abajo: Densidad promedio (barras horizontales) en cada condición y fecha; los valores de cada réplica espacial se indican con círculos (en dos de los grupos de 2021 hay pares de réplicas con el mismo valor estimado). Las líneas punteadas unen a las mismas réplicas espaciales que fueron muestreadas en distintos momentos.

Cuando fue medida en un mismo año (2019), la densidad de colonias de P. bergi no mostró diferencias relevantes entre la estimada con dos grillas cuadradas ubicadas dentro del algarrobal (J y V) o con dos segmentos de 1500 m a lo largo de caminos rurales (picadas) que cruzan el mismo algarrobal (D y F; Fig. 3, izquierda). Combinando, entonces, las estimaciones realizadas con distintos muestreadores en un mismo ambiente (Algarrobal) el cambio interanual en la densidad de colonias activas fue notorio, mostrando un importante decrecimiento entre 2001 y 2025 (Fig. 3, centro). Comparando entre ambientes medidos en un mismo año (2021), las densidades de nidos activos fueron más altas (y algo más variables) en los caminos que cruzan algarrobales que en los de algarrobales pastoreados y médanos (Fig. 3, derecha).

grid.arrange(
  #Densidad promedio de transectas en caminos que cruzan el algarrobal y grillas dentro del algarrobal
  ggplot(subset(Densidades, Year == "2019"), aes(x = Sample, y = Density, group = Sample)) +
    stat_summary(fun.data = mean_se, geom = "pointrange") +
    ylim(c(5, 200)) +
    labs(y = "Density (nests/ha)", x = "Sample type", subtitle = "Algarrobal (feb 2019)") +
    tema + theme(axis.text.x = element_text(size = 10)),
  #Densidad promedio por año
ggplot(subset(Densidades, Habitat == "Algarrobal"), aes(x = as.numeric(Year), y = Density, group = Year)) +
  stat_summary(fun.data = mean_se, geom = "pointrange") +
  ylim(c(10, 360)) + labs(y = "", x = "Year", subtitle = "Algarrobal (roads and grids)") +
  tema + theme(axis.text.x = element_text(size = 10)),
  #Promedio densidad de colonias por ambiente (solo 2021)
  ggplot(subset(Densidades, Year == "2021"), aes(x = Habitat, y = Density, group = Habitat, colour = Habitat)) + 
    stat_summary(fun.data = mean_se, geom = "pointrange") +
    labs(y = "", subtitle = "Rural roads (feb 2021)") +
    ylim(c(3, 150)) +
    scale_x_discrete(limits = c("Algarrobal", "Grazed", "Dunes")) +
    scale_colour_manual(values = c("forestgreen", "darkkhaki", "olivedrab3")) +
    tema + theme(axis.text.x = element_text(size = 10)),
  nrow = 1, widths = c(3, 4, 3))

Figura 3: Densidad promedio de colonias activas de P.bergi estimada en diferentes condiciones y momentos en el Monte central. Izquierda: en parcelas y caminos del algarrobal en febrero 2019. Centro: En el algarrobal a lo largo del tiempo (independientemente del tipo de muestreador). Derecha: en caminos rurales cruzando distintos ambientes durante febrero 2021.

Patrones espaciales de colonias en caminos

En la Figura 4 se presenta la disposición de las colonias activas a lo largo de las dos transectas de 1500 m censadas en 2019 y 2025, y de las nueve de 500 m censadas en 2021.

#disposición de colonias 2019 y 2025
par(mfrow = c(4, 1), mar = c(0, 0, 1, 0))

for(i in 1:length(Pb2019y25)){
 plot(Pb2019y25[[i]], main = names(Pb2019y25)[i], size = 20, col = "forestgreen", shape="crossticks")
 }

#disposición de colonias 2021
par(mfrow = c(9, 1), mar = c(0, 0, 1, 0))

temp.clrs <- data.frame(Name = names(Pb2021), Clr = c(rep("forestgreen", 3), rep("olivedrab3", 3), rep("darkkhaki", 3)))

for(i in names(Pb2021)) {
  plot(Pb2021[[i]], main = i, size = 7, shape="crossticks", col = temp.clrs$Clr[temp.clrs$Name == i])
  if(i == "Grazed (P2)") points(231, 0, pch = 15, col = "white")
}

Figura 4: Disposición de colonias activas de P.bergi a lo largo de segmentos de 1500 m (2019 y 2025, arriba) o 500 m (2021, abajo) de caminos rurales cruzando distintos ambientes del Monte central.

¿Homogeneidad de la densidad a lo largo de las transectas (CSR)?

La densidad de colonias sobre las transectas parece más espacialmente heterogénea en los algarrobales (3/5 transectas presentan fuertes evidencias de no-homogeneidad de la intensidad del proceso espacial unidimensional [CDF tests]) que en los médanos (1/3) y en los algarrobales pastoreados (0/3) (Fig. 5). En todos esos casos las colonias aparecen más agregadas que lo esperado por azar.

El segmento largo (1500 m) del algarrobal con más colonias (D) resultó marcadamente heterogéneo en 2019 y en 2025, con agregamientos de colonias asociados a escalas de al menos hasta 150 m (Fig. 5). La otra (F) no mostró heterogeneidad relevante en 2019, aunque fue algo más heterogénea al reducirse el número de colonias activas en 2025 (Fig. 5).

#intensidad a lo largo de las transectas D y F 2019 y 2025, tests de uniformidad, centered K (two-tailed, mad.test rank 5 of 200 => p = 0.05) y pcf (idem)
par(mar = c(2, 2, 2, .5))
layout(matrix(1:12, 4, 3, byrow = TRUE), widths = c(3, 1, 1))

#loop of transects with vector of values for vertical text position of test results within plots
loop.texts <- data.frame(Lambda = c(0.06, 0.03, 0.12, 0.08), K = c(-5, -8, -10, -15), Pcf = c(1.25, 1.6, 0.95, 1.2))

for(i in 1:length(Pb2019y25)) {
  plot(rhohat(Pb2019y25[[i]], "x"), legend = FALSE, main = names(Pb2019y25)[i], cex.axis = 1, col = "forestgreen", shadecol = alpha("forestgreen", 0.3))
  text(1525, y = loop.texts$Lambda[i], paste("p = ", round(cdf.test(Pb2019y25[[i]], "x", test = "ks", jitter = F)$p.value, 3)), adj = c(1, 1))
  
  #temp guarda los valores de la función para poder usarlos dos veces (plot + mad.test) y los patrones de puntos para linearpcf. fix.n = TRUE ver Baddeley et al 2014 Ecol.Monog.
  temp <- envelope(Pb2019y25[[i]], linearK, nsim = 199, nrank = 5, r = seq(0, 350, 10), verbose = FALSE, fix.n = TRUE, savefuns = TRUE, savepatterns = TRUE)
  plot(temp, cbind(obs - theo, lo = lo - theo, hi = hi - theo) ~ r, main = expression("K"^L*"(r) - r"), legend = FALSE, cex.axis = 1, col = "forestgreen", shadecol = alpha("forestgreen", 0.3))
  text(350, y = loop.texts$K[i], paste("p = ", round(mad.test(temp, use.theory = TRUE)$p.value, 3)), adj = c(1, 1))
  abline(h = 0, lty = 2, col = "forestgreen", lwd = 0.5)
    
  temp2 <- envelope(Pb2019y25[[i]], linearpcf, nsim = 199, nrank = 5, r = seq(0, 350, 10), verbose = FALSE, fix.n = TRUE, simulate = temp, savefuns = TRUE)
  plot(temp2, cbind(obs - mmean + 1, lo = lo - mmean + 1, hi = hi - mmean + 1) ~ r, main = expression("g"^L*"(r)"), legend = FALSE, cex.axis = 1, col = "forestgreen", shadecol = alpha("forestgreen", 0.3))
  text(350, y = loop.texts$Pcf[i], paste("p = ", round(mad.test(temp2, use.theory = TRUE)$p.value, 3)), adj = c(1, 1))
  abline(h = 1, lty = 2, col = "forestgreen", lwd = 0.5)
}

#idem para las transectas 2021
par(mar = c(2, 2, 2, .5))
layout(matrix(1:9, 3, 3, byrow = TRUE), widths = c(3, 1, 1))

#loop of transects with vector of values for vertical text position of test results within plots
loop.texts <- data.frame(Lambda = c(0.12, 0.23, 0.14, 0.07, 0.07, 0.06, 0.06, 0.07, 0.17), K = c(-10, 50, 10, 25, 25, 40, 35, 25, 40), Pcf = c(1.6, 2, 1.3, 0.6, 1.6, 0.4, 1.6, 1.5, 0.7), Clr = c(rep("forestgreen", 3), rep("olivedrab3", 3), rep("darkkhaki", 3)))

for(i in 1:length(Pb2021)) {
  plot(rhohat(Pb2021[[i]], "x"), legend = FALSE, main = paste(names(Pb2021)[i], "- 2021"), cex.axis = 1, col = loop.texts$Clr[i], shadecol = alpha(loop.texts$Clr[i], 0.3))
  text(510, y = loop.texts$Lambda[i], paste("p = ", round(cdf.test(Pb2021[[i]], "x", test = "ks", jitter = F)$p.value, 3)), adj = c(1, 1))
  
  temp <- envelope(Pb2021[[i]], linearK, nsim = 199, nrank = 5, r = seq(0, 120, 5), verbose = FALSE, fix.n = TRUE, savefuns = TRUE, savepatterns = TRUE)
  plot(temp, cbind(obs - theo, lo = lo - theo, hi = hi - theo) ~ r, main = expression("K"^L*"(r) - r"), legend = FALSE, cex.axis = 1, col = loop.texts$Clr[i], shadecol = alpha(loop.texts$Clr[i], 0.3))
  text(120, y = loop.texts$K[i], paste("p = ", round(mad.test(temp, use.theory = TRUE)$p.value, 3)), adj = c(1, 1))
  abline(h = 0, lty = 2, col = loop.texts$Clr[i], lwd = 0.5)
    
  temp2 <- envelope(Pb2021[[i]], linearpcf, nsim = 199, nrank = 5, r = seq(0, 120, 5), verbose = FALSE, fix.n = TRUE, simulate = temp, savefuns = TRUE)
  plot(temp2, cbind(obs - mmean + 1, lo = lo - mmean + 1, hi = hi - mmean + 1) ~ r, main = expression("g"^L*"(r)"), legend = FALSE, cex.axis = 1, col = loop.texts$Clr[i], shadecol = alpha(loop.texts$Clr[i], 0.3))
  text(120, y = loop.texts$Pcf[i], paste("p = ", round(mad.test(temp2, use.theory = TRUE)$p.value, 3)), adj = c(1, 1))
  abline(h = 1, lty = 2, col = loop.texts$Clr[i], lwd = 0.5)
}

Figura 5: Densidad de colonias de P.bergi a lo largo de segmentos de 1500 m (2019 y 2025) o 500 m (2021) de caminos rurales cruzando distintos ambientes del Monte central. Izquierda: intensidad estimada del proceso espacial de patrones de puntos Poisson unidimensional si se asume homogéneo (línea horizontal rayada) o variable (estimada mediante suavizado no paramétrico con kernels gaussianos; sombreado: intervalo de confianza del 95%); el p-valor corresponde a una prueba de homogeneidad basada en la distribución acumulada. Derecha: correlación espacial de las colonias según la función K de Ripley centrada y la función de correlación de pares (g) (líneas: valores observados; sombreado: intervalo crítico del 95% asumiendo procesos Poisson homogéneos) en función de la distancia r entre colonias (en metros); el p-valor corresponde a una prueba global de Diferencia Máxima Absoluta poniendo a prueba que el proceso es Poisson (i.e. sin correlación espacial).

¿Heterogeneidad a lo largo de las transectas sin interacciones entre colonias cercanas (= modelos Poisson inhomogéneos)?

La variación detectada de la densidad de colonias a lo largo de los segmentos de caminos rurales puede ser modelada como un proceso aleatorio sin interacciones (positivas o negativas) entre colonias cercanas (= modelos Poisson inhomogéneos). La intensidad se modificaría, por ejemplo, por variaciones (desconocidas) en características ambientales como la composición, textura y permeabilidad del suelo o la abundancia de alimento.

Todos los casos de agregamientos de colonias no aleatorios que se detectaron (3/5 transectas en Algarrobal, 1/3 transectas en médanos) son compatibles con descripciones no-paramétricas (kernels gaussianos) de la función de intensidad del proceso espacial a lo largo de las transectas (Fig. 6). Es decir, no es necesario invocar interacciones positivas o negativas (distancia-dependientes) entre colonias cercanas para explicar los agregamientos detectados.

### NOTE: Before spatstat.linnet was fixed [still development version by 20230327] for proper estimation of lambda in each iteration, we built different attempts at simulating from the heterogeneous Poisson model (i.e. variable intensities) to build the null models, with density() or rhohat(), either fixing or estimating λ (not true MonteCarlo with λ fixed instead of estimated for each random pattern; n for each iteration was not fixed [fix.n = TRUE didn't work for simulate != NULL). 
### From v.3.1, 'lambda=NULL' in linearKinhom and linearpcfinhom has changed: functions linearK and linearpcf, respectively, are no longer invoked and the intensity lambda is, instead, estimated by kernel smoothing (calling density.lpp with the arguments present in …).

#intensidad a lo largo de las transectas 2019 y 2025 bajo el modelos de suavizado usado por density.lpp, gráficos de modelo con pruebas K-S de bondad de ajuste, gráficos de sobres basados en linearKinhom y linearpcfinhom (a dos colas, centrados en 0 y 1 respectivamente, rank 5 of 200 => p = 0.05) con pruebas MAD pointwise (conservativas) y MAD global corregidas (MonteCarlo de dos etapas)
par(mar = c(2, 2, 1, .5))
layout(matrix(1:12, 4, 3, byrow = TRUE), widths = c(3, 1, 1))

#Before: alternative smoothed λ(x) formulas for Poisson lppm models
#AltModels <- c("log(est_density)", "polynom(x, 5)", "splines::bs(x, df = 5)")
#After spatstat.linnet was fixed (version ≥3.1)
AltModels <- c("log(est_density)")

#vector of values for vertical text position of test results within plots for each transect
  loop.texts <- data.frame(K = c(7, 25, 7, 15), Pcf = c(1.06, 1.12, 0.95, 0.92))

#loop for alternative Poisson models with m smoothed function of x
for(m in AltModels) {
  
  #loop for transects
  for(i in 1:length(Pb2019y25)) {
    
    #estimate density with gaussian kernel and fixed sigma in case the model asks for it
    est_density <- densityfun(Pb2019y25[[i]], sigma = 100, dx = 5)
    
    #define model for transect i with m model formula
    templppm <- lppm(as.formula(paste("Pb2019y25[[", i, "]] ~ ", m)))
    
    #plot the m model, with K-S testing fit of observed against model
    plot(templppm, legend = FALSE, main = names(Pb2019y25)[i], cex.axis = .8, style = "w", adjust = 1000, col = "forestgreen")
    text(1525, y = 0, paste("p = ", round(cdf.test(templppm, "x", test = "ks", jitter = F)$p.value, 3)), adj = c(1, 1))
    
    #pointwise envelope of inhomogeneous linear K function (λ estimated for each simulated pattern)
    temp <- envelope(templppm, linearKinhom, r = seq(0, 350, 10), normpower = 2, nsim = 199, nrank = 5, verbose = FALSE, fix.n = TRUE, leaveoneout = FALSE, savefuns = TRUE, savepatterns = TRUE)
    
    #plot pointwise envelope with (single-staged) MAD test according to linearK [conservative according to Baddeley et al 2017]
    plot(temp, cbind(obs - mmean, lo = lo - mmean, hi = hi - mmean) ~ r, main = "Kinhom-r ~ r", legend = FALSE, cex.axis = 0.8, col = "forestgreen", shadecol = alpha("forestgreen", 0.3))
    text(350, y = loop.texts$K[i], paste("p = ", round(mad.test(temp, use.theory = FALSE)$p.value, 3)), adj = c(1, 0), col = "forestgreen")
    
    #(two-staged) MAD global test based on linearKinhom [Baddeley et al 2017]
    text(350, y = loop.texts$K[i], paste("p = ", round(bits.test(templppm, fun = linearKinhom, exponent = Inf, rinterval = c(0, 350), nsim = 39, leaveoneout = FALSE, use.theory = FALSE, verbose = FALSE)$p.value, 3)), adj = c(1, 1.5))
    abline(h = 0, lty = 2, col = "red", lwd = 0.5)
    
    #idem with linearpcfinhom
    temp2 <- envelope(templppm, linearpcfinhom, r = seq(0, 350, 10), normpower = 2, simulate = temp, nsim = 199, nrank = 5, verbose = FALSE, fix.n = TRUE, leaveoneout = FALSE, savefuns = TRUE)
    plot(temp2, cbind(obs - mmean + 1, lo = lo - mmean + 1, hi = hi - mmean + 1) ~ r, main = "pcfinhom ~ r", legend = FALSE, cex.axis = 0.8, col = "forestgreen", shadecol = alpha("forestgreen", 0.3))
    text(350, y = loop.texts$Pcf[i], paste("p = ", round(mad.test(temp2, use.theory = FALSE)$p.value, 3)), adj = c(1, 0), col = "forestgreen")
    text(350, y = loop.texts$Pcf[i], paste("p = ", round(bits.test(templppm, fun = linearKinhom, exponent = Inf, rinterval = c(0, 350), nsim = 39, leaveoneout = FALSE, use.theory = FALSE, verbose = FALSE)$p.value, 3)), adj = c(1, 1.5))
    abline(h = 1, lty = 2, col = "red", lwd = 0.5)
  }
}

#idem para las transectas 2021
par(mar = c(2, 2, 1, .5))
layout(matrix(1:9, 3, 3, byrow = TRUE), widths = c(3, 1, 1))

#loop of transects with vector of values for vertical text position of test results within plots
loop.texts <- data.frame(Lambda = c(0.12, 0.23, 0.14, 0.07, 0.07, 0.06, 0.06, 0.07, 0.17), K = c(40, -20, 10, 25, 25, 40, 25, 25, 30), Pcf = c(1.7, 1.5, 1.2, 1.4, 1.4, 0.4, 1.5, 1.4, 0.6), Clr = c(rep("forestgreen", 3), rep("olivedrab3", 3), rep("darkkhaki", 3)))

for(m in AltModels) {
  for(i in 1:length(Pb2021)) {
    est_density <- densityfun(Pb2021[[i]], sigma = 50, dx = 2.5)
    templppm <- lppm(as.formula(paste("Pb2021[[", i, "]] ~", m)))
    
    plot(templppm, legend = FALSE, main = names(Pb2021)[i], cex.axis = .8, style = "w", adjust = 300, col = loop.texts$Clr[i], shadecol = alpha(loop.texts$Clr[i], 0.3))
    text(510, y = loop.texts$Lambda[i], paste("p = ", round(cdf.test(templppm, "x", test = "ks", jitter = F)$p.value, 3)), adj = c(1, 1))
    
    temp <- envelope(templppm, linearKinhom, nsim = 199, nrank = 5, r = seq(0, 120, 5), normpower = 2, verbose = FALSE, fix.n = TRUE, leaveoneout = FALSE, savefuns = TRUE, savepatterns = TRUE)
    plot(temp, cbind(obs - mmean, lo = lo - mmean, hi = hi - mmean) ~ r, main = "Kinhom-r ~ r", legend = FALSE, cex.axis = 0.8, col = loop.texts$Clr[i], shadecol = alpha(loop.texts$Clr[i], 0.3))
    text(120, y = loop.texts$K[i], paste("p = ", round(mad.test(temp, use.theory = FALSE)$p.value, 3)), adj = c(1, 0), col = loop.texts$Clr[i])
    text(120, y = loop.texts$K[i], paste("p = ", round(bits.test(templppm, fun = linearKinhom, exponent = Inf, rinterval = c(0, 120), nsim = 39, leaveoneout = FALSE, use.theory = FALSE, verbose = FALSE)$p.value, 3)), adj = c(1, 1.5))
    abline(h = 0, lty = 2, col = "red", lwd = 0.5)
    
    temp2 <- envelope(templppm, linearpcfinhom, nsim = 199, nrank = 5, r = seq(0, 120, 5), normpower = 2, simulate = temp, verbose = FALSE, fix.n = TRUE, leaveoneout = FALSE, savefuns = TRUE)
    plot(temp2, cbind(obs - mmean + 1, lo = lo - mmean + 1, hi = hi - mmean + 1) ~ r, main = "pcfinhom ~ r", legend = FALSE, cex.axis = 0.8, col = loop.texts$Clr[i], shadecol = alpha(loop.texts$Clr[i], 0.3))
    text(120, y = loop.texts$Pcf[i], paste("p = ", round(mad.test(temp2, use.theory = FALSE)$p.value, 3)), adj = c(1, 0), col = loop.texts$Clr[i])
    text(120, y = loop.texts$Pcf[i], paste("p = ", round(bits.test(templppm, fun = linearKinhom, exponent = Inf, rinterval = c(0, 120), nsim = 39, leaveoneout = FALSE, use.theory = FALSE, verbose = FALSE)$p.value, 3)), adj = c(1, 1.5))
    abline(h = 1, lty = 2, col = "red", lwd = 0.5)
  }
}

Figura 6: Variación estimada de la densidad de colonias de P.bergi a lo largo de los segmentos de 1500 m (2019 y 2025) y 500 m (2021) de caminos rurales cruzando distintos ambientes del Monte central, suponiendo que no hay interacciones relevantes (positivas o negativas) entre colonias cercanas (modelos Poisson heterogéneos). Izquierda: intensidad estimada mediante suavizado no paramétrico con kernels gaussianos; el p-valor corresponde a una prueba K-S de bondad de ajuste. Derecha: autocorrelación espacial residual de las colonias (i.e., descontando la variación de la intensidad estimada) según la función K de Ripley centrada y la función de correlación de pares (g) (líneas: valores observados; sombreado: intervalos críticos punto por punto del 95%) en función de la distancia r entre colonias (en metros); los p-valores corresponden a pruebas globales de Máxima Desviación Absoluta (en color, probablemente conservativas) y pruebas globales corregidas (en negro; MonteCarlo de dos etapas).

#gráficos de sobres basados en linearKinhom y linearpcfinhom (a dos colas, centrados en 0 y 1 respectivamente, rank 5 of 200 => p = 0.05) con pruebas MAD pointwise (conservativas) y MAD global corregidas (MonteCarlo de dos etapas) SOLO PARA LAS QUE SE RECHAZÓ CSR
par(mar = c(4, 2, 2, .5), oma = c(2, 0, 2, 0))
layout(matrix(c(0, 1:4, 0, 5:10), 2, 6, byrow = T))

#vector of values for vertical text position of test results within plots for each transect
  loop.texts <- data.frame(K = c(7, 25, 7, 15), Pcf = c(1.06, 1.12, 0.95, 0.92))
#loop for D (2019 & 2025) & F(2025) transects
for(i in c(1, 2)) {
  #estimate density with gaussian kernel and fixed sigma in case the model asks for it
  est_density <- densityfun(Pb2019y25[[i]], sigma = 100, dx = 5)
  #define model for transect i as the estimated density
  templppm <- lppm(Pb2019y25[[i]] ~ log(est_density))
   
  #pointwise envelope of inhomogeneous linear K function (λ estimated for each simulated pattern)
  temp <- envelope(templppm, linearKinhom, nsim = 199, nrank = 5, r = seq(0, 350, 10), normpower = 2, verbose = FALSE, fix.n = TRUE, leaveoneout = FALSE, savefuns = TRUE, savepatterns = TRUE)
  plot(temp, cbind(obs - mmean, lo = lo - mmean, hi = hi - mmean) ~ r, main = expression("K"[inh]^L*"(r) - r"), xlab = "", legend = FALSE, cex.axis = 1, col = "forestgreen", shadecol = alpha("forestgreen", 0.3))
  text(350, y = loop.texts$K[i], paste("p = ", round(mad.test(temp, use.theory = FALSE)$p.value, 3)), adj = c(1, 0), col = "forestgreen")
  text(350, y = loop.texts$K[i], paste("p = ", round(bits.test(templppm, fun = linearKinhom, exponent = Inf, rinterval = c(0, 350), nsim = 39, leaveoneout = FALSE, use.theory = FALSE, verbose = FALSE)$p.value, 3)), adj = c(1, 1.5))
  abline(h = 0, lty = 2, col = "red", lwd = 0.5)
  mtext(names(Pb2019y25)[i], side = 3, at = 250, adj = 0, line = 2, cex = 1)
    
  temp2 <- envelope(templppm, linearpcfinhom, nsim = 199, nrank = 5, r = seq(0, 350, 10), normpower = 2, simulate = temp, verbose = FALSE, fix.n = TRUE, leaveoneout = FALSE, savefuns = TRUE)
  plot(temp2, cbind(obs - mmean + 1, lo = lo - mmean + 1, hi = hi - mmean + 1) ~ r, main = expression("g"[inh]^L*"(r)"), xlab = "", legend = FALSE, cex.axis = 1, col = "forestgreen", shadecol = alpha("forestgreen", 0.3))
  text(350, y = loop.texts$Pcf[i], paste("p = ", round(mad.test(temp2, use.theory = FALSE)$p.value, 3)), adj = c(1, 0), col = "forestgreen")
  text(350, y = loop.texts$Pcf[i], paste("p = ", round(bits.test(templppm, fun = linearKinhom, exponent = Inf, rinterval = c(0, 350), nsim = 39, leaveoneout = FALSE, use.theory = FALSE, verbose = FALSE)$p.value, 3)), adj = c(1, 1.5))
  abline(h = 1, lty = 2, col = "red", lwd = 0.5)
}

#loop of transects with vector of values for vertical text position of test results within plots
loop.texts <- data.frame(Lambda = c(0.12, 0.23, 0.14, 0.07, 0.07, 0.06, 0.06, 0.07, 0.17), K = c(40, -20, 10, 25, 25, 40, 25, 25, 30), Pcf = c(1.7, 1.5, 1.2, 1.4, 1.4, 0.4, 1.5, 1.4, 0.6), Clr = c(rep("forestgreen", 3), rep("olivedrab3", 3), rep("darkkhaki", 3)))

#loop for D3, F1 & M (2021) transects
  for(i in c(1, 2, 9)) {
    #estimate density with gaussian kernel and fixed sigma in case the model asks for it
    est_density <- densityfun(Pb2021[[i]], sigma = 50, dx = 2.5)
    #define model for transect i as the estimated density
    templppm <- lppm(Pb2021[[i]] ~ log(est_density))
     
    #pointwise envelope of inhomogeneous linear K function (λ estimated for each simulated pattern)
    temp <- envelope(templppm, linearKinhom, nsim = 199, nrank = 5, r = seq(0, 120, 5), normpower = 2, verbose = FALSE, fix.n = TRUE, leaveoneout = FALSE, savefuns = TRUE, savepatterns = TRUE)
    plot(temp, cbind(obs - mmean, lo = lo - mmean, hi = hi - mmean) ~ r, main = expression("K"[inh]^L*"(r) - r"), xlab = "", legend = FALSE, cex.axis = 1, col = loop.texts$Clr[i], shadecol = alpha(loop.texts$Clr[i], 0.3))
    text(120, y = loop.texts$K[i], paste("p = ", round(mad.test(temp, use.theory = FALSE)$p.value, 3)), adj = c(1, 0), col = loop.texts$Clr[i])
    text(120, y = loop.texts$K[i], paste("p = ", round(bits.test(templppm, fun = linearKinhom, exponent = Inf, rinterval = c(0, 120), nsim = 39, leaveoneout = FALSE, use.theory = FALSE, verbose = FALSE)$p.value, 3)), adj = c(1, 1.5))
    abline(h = 0, lty = 2, col = "red", lwd = 0.5)
    mtext(paste(names(Pb2021)[i], "- 2021"), side = 3, at = 80, adj = 0, line = 2, cex = 1)
    
    temp2 <- envelope(templppm, linearpcfinhom, nsim = 199, nrank = 5, r = seq(0, 120, 5), normpower = 2, simulate = temp, verbose = FALSE, fix.n = TRUE, leaveoneout = FALSE, savefuns = TRUE)
    plot(temp2, cbind(obs - mmean + 1, lo = lo - mmean + 1, hi = hi - mmean + 1) ~ r, main = expression("g"[inh]^L*"(r)"), xlab = "", legend = FALSE, cex.axis = 1, col = loop.texts$Clr[i], shadecol = alpha(loop.texts$Clr[i], 0.3))
    text(120, y = loop.texts$Pcf[i], paste("p = ", round(mad.test(temp2, use.theory = FALSE)$p.value, 3)), adj = c(1, 0), col = loop.texts$Clr[i])
    text(120, y = loop.texts$Pcf[i], paste("p = ", round(bits.test(templppm, fun = linearKinhom, exponent = Inf, rinterval = c(0, 120), nsim = 39, leaveoneout = FALSE, use.theory = FALSE, verbose = FALSE)$p.value, 3)), adj = c(1, 1.5))
    abline(h = 1, lty = 2, col = "red", lwd = 0.5)
  }

Figura 7: Idem Fig. 6 derecha, pero solo para las transectas que se consideraron heterogéneas (ver Fig. 5).

Conclusiones

• densidades similares en algarrobales protegidos y caminos rurales (disturbados) que los cruzan
  
• densidades menores en caminos que cruzan algarrobales pastoreados y dunas que en los del algarrobal protegido
  
• reducción notoria de la densidad de colonias entre 2001 y 2025
  
• hay importantes variaciones de densidad a lo largo de los segmentos de caminos rurales, en particular de los que cruzan el algarrobal. Esto podría deberse a variaciones ambientales como características del suelo, de la vegetación o de los recursos alimentarios asociados
  
• el agregamiento de colonias detectado en varias de los segmentos censados no implica necesariamente la existencia de interacciones positivas (dependientes de la distancia) entre colonias cercanas

---

sessionInfo()

## R version 4.5.0 (2025-04-11)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.2 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.12.0 
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0  LAPACK version 3.12.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_AU.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_AU.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_AU.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_AU.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: America/Argentina/Buenos_Aires
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] scales_1.4.0           png_0.1-8              gridExtra_2.3         
##  [4] ggplot2_3.5.2          dplyr_1.1.4            spatstat_3.3-2        
##  [7] spatstat.linnet_3.2-5  spatstat.model_3.3-5   rpart_4.1.24          
## [10] spatstat.explore_3.4-2 nlme_3.1-168           spatstat.random_3.3-3 
## [13] spatstat.geom_3.3-6    spatstat.univar_3.1-2  spatstat.data_3.1-6   
## 
## loaded via a namespace (and not attached):
##  [1] sass_0.4.10           generics_0.1.3        tensor_1.5           
##  [4] lattice_0.22-5        digest_0.6.37         magrittr_2.0.3       
##  [7] spatstat.utils_3.1-3  RColorBrewer_1.1-3    evaluate_1.0.3       
## [10] grid_4.5.0            bookdown_0.43         fastmap_1.2.0        
## [13] jsonlite_2.0.0        Matrix_1.7-3          spatstat.sparse_3.1-0
## [16] mgcv_1.9-1            jquerylib_0.1.4       abind_1.4-8          
## [19] cli_3.6.5             rlang_1.1.6           polyclip_1.10-7      
## [22] splines_4.5.0         withr_3.0.2           cachem_1.1.0         
## [25] yaml_2.3.10           tools_4.5.0           deldir_2.0-4         
## [28] vctrs_0.6.5           R6_2.6.1              lifecycle_1.0.4      
## [31] pkgconfig_2.0.3       bslib_0.9.0           pillar_1.10.2        
## [34] gtable_0.3.6          glue_1.8.0            xfun_0.52            
## [37] tibble_3.2.1          tidyselect_1.2.1      knitr_1.50           
## [40] farver_2.1.2          goftest_1.2-3         htmltools_0.5.8.1    
## [43] labeling_0.4.3        rmarkdown_2.29        compiler_4.5.0