Pogos en parcela

---

Pogonomyrmex spp. in ~15 ha of algarrobal [Ñacuñán 1999–2003]

Methods

We searched extensively for Pogonomyrmex spp. colonies in ~15 ha of algarrobal during three consecutive growing seasons (spring to early autumn) of 1999–2002, when colonies show external activity (Pol & Lopez de Casenave 2004). During four times per activity season (October, December, February and April) we walked the area at times of day when Pogonomyrmex foragers are known to be active (cita Pol). Most colonies are only recognizable either by detecting workers at the openings doing nest maintenance or by detecting and following returning foragers. Colonies were marked permanently with a numbered stick and mapped; colonies on dirt roads or roadsides were ignored. Ten colonies detected on the first sampling season were not mapped adequately; they were counted but removed for spatial analyses. During each sampling period all previously detected colonies were checked, recording their status (with or without external activity) and, if they have moved the colony entrance, the distance of its displacement. Any colony >5 m from any previously detected colony was considered a new colony. Inactive colonies were checked in the following sampling periods (including the summer 2002-2003) and considered “dead” if they showed no further external activity.

Libraries and databases

require(scales)
require(tidyr)
require(dplyr)
require(ggplot2)
require(gridExtra)
require(spatstat)
tema <- theme_bw(base_size = 12) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "none")

#disables scientific notation
options(scipen = 6)

# importación de planilla de destino de nidos; rangos de distancia; cada fila es una observación [d: descubierto, a: activo, n: no activo, c: cerrado ("para siempre")] de cada nido detectado en algún momento
Destino_nidos <- read.delim("~/Documents/Ecodes/Hormigas/Espacial/Destino_nidos.csv", comment.char = "#")
Destino_nidos$Dist <- as.factor(Destino_nidos$Dist)
Destino_nidos$Sp <- as.factor(Destino_nidos$Sp)
levels(Destino_nidos$Dist) <- c("0–25 m", "25–40 m", "40–90 m", "90–140 m", ">140 m")
# identificación Temporada activo (TRUE if "d" or "a") en $Temp1, $Temp2 y $Temp3
Destino_nidos$Temp1 <- apply(Destino_nidos[4:7], 1, function(x) {any("d" %in% x, "a" %in% x)})
Destino_nidos$Temp2 <- apply(Destino_nidos[8:11], 1, function(x) {any("d" %in% x, "a" %in% x)})
Destino_nidos$Temp3 <- apply(Destino_nidos[12:15], 1, function(x) {any("d" %in% x, "a" %in% x)})

# a formato largo y eliminación de celdas vacías (sacando filas): cada fila es un nido observado en un muestreo
Dest_nidos_long <- Destino_nidos %>%
  select(!(starts_with("Temp"))) %>%
  pivot_longer(cols = starts_with(c("Oct", "Dic", "Feb", "Abr")), names_to = "Fecha", names_transform = as.factor, values_to = "Status")
Dest_nidos_long <- subset(Dest_nidos_long, Status != "")
Dest_nidos_long$Status <- as.factor(Dest_nidos_long$Status)
levels(Dest_nidos_long$Dist) <- c("0–25 m", "25–40 m", "40–90 m", "90–140 m", ">140 m")
# identificación temporada (1 a 3) en $Season
Dest_nidos_long$Season <- as.factor(Dest_nidos_long$Fecha)
levels(Dest_nidos_long$Season) <- list("1" = c("Oct99", "Dic99", "Feb00", "Abr00"), "2" = c("Oct00", "Dic00", "Feb01", "Abr01"), "3" = c("Oct01", "Dic01", "Feb02", "Abr02"))
# identificación de mes (Oct, Dic, Feb, Abr) en $Month
Dest_nidos_long$Month <- as.factor(Dest_nidos_long$Fecha)
levels(Dest_nidos_long$Month) <- list("Oct" = c("Oct99", "Oct00", "Oct01"), "Dic" = c("Dic99", "Dic00", "Dic01"), "Feb" = c("Feb00", "Feb01", "Feb02"), "Abr" = c("Abr00", "Abr01", "Abr02"))

# importación de planilla de distancia de desplazamiento de entrada de colonias: cada fila es una observación de una colonia que desplazó su entrada entre periodos
Desplazamiento_nidos <- read.delim("~/Documents/Ecodes/Hormigas/Espacial/Desplazamiento_nidos.csv", comment.char = "#")
Desplazamiento_nidos$Month <- ordered(factor(Desplazamiento_nidos$Fecha, labels = c("Apr", "Apr", "Apr", "Dec", "Dec", "Dec", "Feb", "Feb", "Feb", "Oct")), levels = c("Oct", "Dec", "Feb", "Apr"))

# importación de planillas de coordenadas de nidos en mapa-esquema y de coordenadas de vértices de área de estudio y picadas en el área; no todos los nidos observados fueron mapeados
Mapa_nidos <- read.delim("~/Documents/Ecodes/Hormigas/Espacial/Mapanidos.txt", comment.char = "#")
Extra_mapa_nidos <- read.delim("~/Documents/Ecodes/Hormigas/Espacial/Mapanidosextra.txt", comment.char = "#")
# desplazamiento de coordenadas relativas (x - 2300; y - 1200) para trabajar con números más chicos
Mapa_nidos$x <- Mapa_nidos$x - 2300
Mapa_nidos$y <- Mapa_nidos$y - 1200
Extra_mapa_nidos$x <- Extra_mapa_nidos$x - 2300
Extra_mapa_nidos$y <- Extra_mapa_nidos$y - 1200
# creates an object of class "ppp" representing a spatial point pattern [spatstat] (marked point pattern with species as marks [factor])
xyNidos <- ppp(x = Mapa_nidos$x, y = Mapa_nidos$y, poly = Extra_mapa_nidos[1:6, 1:2], unitname = "m", marks = as.factor(Mapa_nidos$Sp))

#creates an object of class "psp" representing a line segment pattern [spatstat] (segmentos de picada, marca es el nombre, misma ventana de observación que xyNidos)
xyPicadas <- psp(x0 = Extra_mapa_nidos[c(7,8,10), 1], y0 = Extra_mapa_nidos[c(7,8,10), 2], x1 = Extra_mapa_nidos[c(8,9,11), 1], y1 = Extra_mapa_nidos[c(8,9,11), 2], window = owin(poly = Extra_mapa_nidos[1:6, 1:2]), marks = droplevels(as.factor(Extra_mapa_nidos$Objeto[8:10])), check = TRUE)

#función de distancias a Picadas
DistPicadas <- distfun(xyPicadas)

#cálculo de distancias de cada nido mapeado a las picadas
Mapa_nidos$distpic <- with(Mapa_nidos, round(DistPicadas(x, y), digits = 2))

#cálculo de "número de muestreos activo" de cada colonia para incluir en mapa [el R50 está en "mapa" pero no tiene "destino" (fue descubierto en 2002-2003 => Active = 0)]
Mapa_nidos <- Dest_nidos_long %>%
  group_by(Nido) %>%
  filter(Status == "d" | Status == "a") %>%
  summarise(Active = n(), Season1 = any(Season == 1), Season2 = any(Season == 2), Season3 = any(Season == 3)) %>%
  right_join(Mapa_nidos)

Mapa_nidos$Active[is.na(Mapa_nidos$Active)] <- 0

#identificación de nidos que se confirmaron cerrados en algún momento del muestreo 
Mapa_nidos <- Dest_nidos_long %>%
  group_by(Nido) %>%
  summarise(Closed = any(Status == "c")) %>%
  right_join(Mapa_nidos)

# "ppp" objects for spatial point patterns [spatstat] by Season (i.e., nests active at least once in that Season), with marks for species
NidosS1 <- with(subset(Mapa_nidos, Season1), ppp(x = x, y = y, poly = Extra_mapa_nidos[1:6, 1:2], unitname = "m", marks = as.factor(Sp)))

#to add "PI" as mark (with no colonies)
NidosS1$marks <- factor(NidosS1$marks, levels = c("PI", levels(NidosS1$marks)))

NidosS2 <- with(subset(Mapa_nidos, Season2), ppp(x = x, y = y, poly = Extra_mapa_nidos[1:6, 1:2], unitname = "m", marks = as.factor(Sp)))
NidosS3 <- with(subset(Mapa_nidos, Season3), ppp(x = x, y = y, poly = Extra_mapa_nidos[1:6, 1:2], unitname = "m", marks = as.factor(Sp)))
NidosSeasons <- list(NidosS1, NidosS2, NidosS3)

#hyperframe to model seasons as replicates
xyNidosbySeason <- hyperframe(Nidos = solist(NidosS1, NidosS2, NidosS3),
                              Picadas = solist(DistPicadas, DistPicadas, DistPicadas), 
                              Season = as.factor(c("1", "2", "3")))

# function to average the Kcross values ij and ji (they may be asymmetrical due to edge corrections; see p. 596 in Baddeley et al(2016))
Kcross.avrg <- function(X, i, j, r, sigma) {
  Kij <- Kcross(X = X, i = i, j = j, r = r, sigma = sigma)
  Kji <- Kcross(X = X, i = j, j = i, r = r, sigma = sigma)
  eval.fv((Kij + Kji) / 2)
}

# function to average the pcfcross values ij and ji (they may be asymmetrical due to edge corrections; see p. 596 in Baddeley et al(2016))
pcfcross.avrg <- function(X, i, j, r, sigma) {
  pcfij <- pcfcross(X = X, i = i, j = j, r = r, sigma = sigma, lambdaI = NULL, lambdaJ = NULL)
  pcfji <- pcfcross(X = X, i = j, j = i, r = r, sigma = sigma, lambdaI = NULL, lambdaJ = NULL)
  eval.fv((pcfij + pcfji) / 2)
}

Results

Abundance and spatio-temporal dynamics of Pogonomyrmex colonies

We detected 105 Pogonomyrmex spp colonies in ~15 ha of algarrobal during three consecutive years, though no more than half of those were active at any particular time (max = 51 in Feb01, Fig. 1A). Of those colonies, 73 (70%) had “closed” by the end of the third season (i.e., never seen active again, including the 2002-2003 activity season; Fig. 2).

As expected, along the three seasons the number of active colonies was lowest in spring, increased towards December, peaked in February, and decreased again by autumn (Fig. 1B) though with variable numbers in each activity season (Fig. 1C). Assuming complete censuses during each activity season (i.e., after its four sampling periods), many colonies appeared (48 in Season 2 and 22 in Season 3) or disappeared (10 during Season 1, 35 in Season 2 and 28 in Season 3) each year. In fact, only 7 of the 35 colonies detected during the first year showed external activity during Season 3 (5 P. mendozanus and 2 P. rastratus); none of them were inactive in Season 2.

# descriptivos de nidos activos
grid.arrange(
  # A. nidos activos por muestreo
  ggplot(subset(Dest_nidos_long, Status == "a" | Status == "d"), aes(x = Fecha, fill = Sp)) +
    geom_histogram(stat = "count") +
    labs(y = "Active colonies", x = "Date", tag = "A") +
    scale_x_discrete(limits = c("Oct99", "Dic99", "Feb00", "Abr00", "", "Oct00", "Dic00", "Feb01", "Abr01", "", "Oct01", "Dic01", "Feb02", "Abr02")) +
    geom_vline(xintercept = c(5, 10), linetype = 2) +
    scale_fill_manual(values = c("darkgreen", "darkblue", "darkred")) +
    tema + theme(axis.text.x = element_text(size = 9)),
  # B. nidos activos por mes (integrando las 3 temporadas)
  ggplot(subset(Dest_nidos_long, Status == "a" | Status == "d"), aes(x = Month, fill = Sp)) +
    geom_histogram(stat = "count", aes(y = after_stat(count / 3))) +
    labs(y = "(average per season)", x = "Month (1999–2002)", tag = "B") +
    scale_fill_manual(values = c("darkgreen", "darkblue", "darkred")) +
    tema,
  # C. nidos activos por temporada (integrando las 4 fechas de muestreo)
  ggplot(subset(Dest_nidos_long, Status == "a" | Status == "d"), aes(x = Season, fill = Sp)) +
    geom_histogram(stat = "count", aes(y = after_stat(count / 3))) +
    labs(y = "(average per date)", x = "Activity Season (Oct–Apr)", tag = "C") +
    scale_fill_manual(values = c("darkgreen", "darkblue", "darkred")) +
    tema,
  layout_matrix = rbind(c(1, 1, 2), c(1, 1, 3))
)

Figure 1: Number of active Pogonomyrmex spp. colonies detected in each of four sampling dates along three consecutive activity seasons (Oct–Apr 1999–2002) in ~15 ha of algarrobal.

# nidos cerrados por fecha
ggplot(subset(Dest_nidos_long, Status == "c"), aes(x = Fecha, fill = Sp)) +
  geom_histogram(stat = "count") +
  labs(y = "'Closed' colonies", x = "First date inactive") +
  scale_x_discrete(limits = c("Oct99", "Dic99", "Feb00", "Abr00", "", "Oct00", "Dic00", "Feb01", "Abr01", "", "Oct01", "Dic01", "Feb02", "Abr02")) +
  geom_vline(xintercept = c(5, 10), linetype = 2) + 
  scale_fill_manual(values = c("darkgreen", "darkblue", "darkred")) +
  tema + theme(axis.text.x = element_text(size = 8))

Figure 2: Number of detected Pogonomyrmex spp. colonies in ~15 ha of algarrobal that closed (i.e., showed no further external activity or moved its entrance >5 m from its original place and was considered a different colony) along three consecutive activity seasons (Oct–Apr 1999–2002) and a status recheck during the 2002–2003 activity season

Nearly half of the colonies detected along the three years within this area of algarrobal were P. rastratus (n = 52; 49.5%), a third were P. mendozanus (35; 33.3%)) and the rest, P. inermis (18; 17.1%), which increased its proportional abundance of active colonies with time (from no detected colonies in Season 1 to 28% by Season 3; Fig. 1C).

Colonies frequently moved their entrances between sampling periods. We recorded 0.12 – 4.04 m (median = 1 m) displacements of the entrance of the colonies in 18.3% of the observations on previously marked active colonies, particularly by the end of the activity seasons (52.4% of colony displacements occurred between February and April; (Fig. 3). The distributions of colony displacement distances differed among species (Kruskal-Wallis test: 6.386, df = 2, p = 0.041, with shorter movements by P. rastratus (median = 0.68 m, n = 18) than P. mendozanus (1.5 m, n = 17).

grid.arrange(
  # A. nidos desplazados por muestreo
  ggplot(Desplazamiento_nidos, aes(x = Fecha, fill = Sp)) +
    geom_bar() +
    labs(y = "Nest displacements", x = "Date", tag = "A") +
    scale_x_discrete(limits = c("Oct99", "Dec99", "Feb00", "Apr00", "", "Oct00", "Dec00", "Feb01", "Apr01", "", "Oct01", "Dec01", "Feb02", "Apr02")) +
    geom_vline(xintercept = c(5, 10), linetype = 2) +
    scale_fill_manual(values = c("darkgreen", "darkblue", "darkred")) +
    tema + theme(axis.text.x = element_text(size = 9)),
  # B. distancia de desplazamiento por especie
  ggplot(Desplazamiento_nidos, aes(x = Sp, y = Mov_dist, fill = Sp, colour = Sp)) +
    geom_boxplot(alpha = 0.5) +
    geom_point() +
    labs(y = "Distance (m)", x = "Species", tag = "B") +
    scale_fill_manual(values = c("darkgreen", "darkblue", "darkred")) +
    scale_colour_manual(values = c("darkgreen", "darkblue", "darkred")) +
    tema,
  # C. distancia de desplazamiento por mes
  ggplot(Desplazamiento_nidos, aes(x = Month, y = Mov_dist)) +
    geom_boxplot(alpha = 0.5) +
    geom_point() +
    labs(y = "Distance (m)", x = "Month (1999–2002)", tag = "C") +
    scale_fill_manual(values = c("darkgreen", "darkblue", "darkred")) +
    scale_colour_manual(values = c("darkgreen", "darkblue", "darkred")) +
    tema,
  layout_matrix = rbind(c(1, 1, 2), c(1, 1, 3))
)

Figure 3: Number of detected Pogonomyrmex spp. colonies in ~15 ha of algarrobal that moved (i.e., displaced their nest entrances) during three consecutive activity seasons (Oct–Apr 1999–2002), with their measured distance of displacement (m) by species and month.

Spatial patterns by species

Pogonomyrmex colonies were not distributed homogeneously in the algarrobal: we were not able to find any active colonies in approximately half of the sampled area along the three consecutive years (Fig. 4). Within the restricted area with colonies, P. inermis and P. mendozanus appeared mostly aggregated within 2–3 patches, while P. rastratus was not only more abundant but more broadly distributed.

# mapa - figura
ggplot(Extra_mapa_nidos[c(1:6, 1),], aes(x = x, y = y)) + geom_polygon(alpha = 0.2) + geom_path(data = Extra_mapa_nidos[7:9,]) + geom_path(data = Extra_mapa_nidos[10:11,]) + geom_point(data = Mapa_nidos, aes(colour = Sp, size = Active, shape = as.factor(Closed))) + tema + theme(panel.border = element_blank(), axis.title.x = element_blank(), axis.title.y = element_blank()) + coord_fixed() + scale_shape_manual(values=c(19, 21)) + scale_colour_manual(values = c("darkgreen", "darkblue", "darkred"))

Figure 4: Active colonies of Pogonomyrmex mendozanus (blue), P. rastratus (red) and P. inermis (green) in the sampled area of algarrobal during three consecutive activity seasons (spring-autumn 1999–2002). Colonies located on the dirt roads (black lines) crossing the area were not considered. Dot size is proportional to the number of sampling dates (1–12) in which colonies shown external activity. Empty circles are colonies that ‘closed’ (i.e., showed no further external activity or moved its entrance >5 m from its original place and was considered a different colony) during the sampling period (including a status recheck during the 2002–2003 activity season)

Optimal bandwith values of the gaussian kernel estimators (sigma) for a smoothed description of their densities varies within the 14–49 m range (sometimes with similar likelihoods for neighbouring values; Fig. 5, above), so we choose σ = 25 m as a compromise value. Considering all the positions where active colonies were detected sometime along the sampling period, the estimated intensities of their spatial processes appear notoriously heterogeneous (Fig. 5, below), with two sectors close to the dirt roads concentrating most ant colonies.

#determinación de bandwidth óptimos (bw.ppl y bw.diggle) para los kernels gaussianos de los mapas de densidad (curvas de error en función del bandwidth). bw.ppl assumes an heterogeneous Poisson process
#todos los nidos
par(mar = c(2, 2, 3, .5))
layout(matrix(c(1:4, rep.int(5, 4)), nrow = 2, byrow = TRUE), heights=c(1,4))

plot(bw.ppl(xyNidos, sigma = 1:50), xlim = c(5, 50), main = "bw.ppl(Pogos spp)")
plot(bw.ppl(xyNidos[xyNidos$marks == "PM"], sigma = 1:50), col = "darkblue", xlim = c(5, 50), main = "bw.ppl(PM)")
plot(bw.ppl(xyNidos[xyNidos$marks == "PI"], sigma = 1:50), col = "darkgreen", xlim = c(5, 50), main = "bw.ppl(PI)")
plot(bw.ppl(xyNidos[xyNidos$marks == "PR"], sigma = 1:50), col = "darkred", xlim = c(5, 50), main = "bw.ppl(PR)")

#idem with bw.diggle, que asume un modelo contagioso
#optimal values => All: 15.4, PM: 14.4, PI: 8.2, PR: 16.3
#plot(bw.diggle(xyNidos), xlim = c(5, 25), main = "bw.diggle(Pogos spp)")
#plot(bw.diggle(xyNidos[xyNidos$marks == "PM"]), col = "darkblue", xlim = c(5, 25), main = "bw.diggle(PM)")
#plot(bw.diggle(xyNidos[xyNidos$marks == "PI"]), col = "darkgreen", xlim = c(5, 25), main = "bw.diggle(PI)")
#plot(bw.diggle(xyNidos[xyNidos$marks == "PR"]), col = "darkred", xlim = c(5, 25), main = "bw.diggle(PR)")

# mapas de densidad (kernel gaussiano)
plot(density(xyNidos[xyNidos$marks == "PM"], sigma = 25, diggle = TRUE, eps = 1), col = colorRampPalette(c(alpha("darkblue", 0.025), "darkblue"), alpha = TRUE)(256), main = "Pogonomyrmex densities (sigma = 25 m)", ribbon = FALSE)
plot(density(xyNidos[xyNidos$marks == "PR"], sigma = 25, diggle = TRUE), col = colorRampPalette(c(alpha("darkred", 0.025), "darkred"), alpha = TRUE)(256), add = TRUE)
plot(density(xyNidos[xyNidos$marks == "PI"], sigma = 25, diggle = TRUE), col = colorRampPalette(c(alpha("darkgreen", 0.025), "darkgreen"), alpha = TRUE)(256), add = TRUE)
lines(Extra_mapa_nidos$x[7:9], Extra_mapa_nidos$y[7:9], lwd = 4)
lines(Extra_mapa_nidos$x[10:11], Extra_mapa_nidos$y[10:11], lwd = 4)

Figure 5: Density of Pogonomyrmex colonies in ~15 ha of algarrobal, by species. Above: optimal bandwith for gaussian kernel estimation for all colonies and foreach species separately. Below: estimated densities of colony positions (sigma = 25) aggregating the three consecutive activity seasons (spring-autumn 1999–2002) they were sampled.

#segregation test with fixed sigma [similar result with free sigma but lots of warnings]
xyNidos.segr.test <- segregation.test(xyNidos, nsim = 199, sigma = 25, verbose = FALSE)

In fact, summary functions such as Ripley’s K and pair correlation are much higher for every species than expected under CSRI processes, indicating intra-specific spatial aggregation (Fig. 6).

Moreover, the three species were not independent but spatially segregated (Monte Carlo segregation test: T = 27.464, p = 0.005), particularly because of P. inermis occupying almost monospecific patches (Fig. 7).

par(mar = c(2, 2, 1, .5))
layout(matrix(1:12, 4, 3, byrow = TRUE), widths = c(1.5, 1, 1))

#all spp together
  plot(density(xyNidos, sigma = 25, diggle = TRUE, eps = 1), main = "density All", col = colorRampPalette(c(alpha("black", 0.1), "black"), alpha = TRUE)(256), legend = FALSE, cex.axis = 0.8)
  text(250, y = 480, paste("p = ", round(quadrat.test(xyNidos, nx = 4, method = "MonteCarlo", CR = 0)$p.value, 3)), col = "black", adj = c(1, 1))
  
  temp <- envelope(xyNidos, Kest, nsim = 199, nrank = 5, r = seq(0, 100, 5), verbose = FALSE, fix.n = TRUE, savefuns = TRUE, savepatterns = TRUE)
  plot(temp, cbind(obs - theo, lo = lo - theo, hi = hi - theo) ~ r, main = "K-r ~ r", legend = FALSE, cex.axis = 0.8, col = "black", shadecol = alpha("black", 0.3))
  text(100, y = 15000, paste("p = ", round(mad.test(temp, use.theory = TRUE)$p.value, 3)), adj = c(1, 1))
    abline(h = 0, lty = 2, col = "black", lwd = 0.5)
    
  temp2 <- envelope(xyNidos, pcf, nsim = 199, nrank = 5, r = seq(0, 100, 5), verbose = FALSE, simulate = temp,  savefuns = TRUE)
  plot(temp2, cbind(obs - mmean + 1, lo = lo - mmean + 1, hi = hi - mmean + 1) ~ r, main = "pcf ~ r", legend = FALSE, cex.axis = 0.8, col = "black", shadecol = alpha("black", 0.3))
  text(100, y = 4, paste("p = ", round(mad.test(temp2, use.theory = TRUE)$p.value, 3)), adj = c(1, 1))
  abline(h = 1, lty = 2, col = "black", lwd = 0.5)


#loop species with vector of values for vertical text position of test results within plots

#temp dataframe with vectors of values for titles and vertical text position of test results within plots and counter to select each of them
loop.texts <- data.frame(Main = c("density PI", "density PM", "density PR"), K = c(35000, 50000, 25000), PCF = c(17, 15, 4.5), Col = c("darkgreen", "darkblue", "darkred"))
loop.count <- 0

for(i in split(xyNidos)) {
  loop.count <- loop.count + 1
  
  plot(density(i, sigma = 25, diggle = TRUE, eps = 1), main = loop.texts$Main[loop.count], col = colorRampPalette(c(alpha(loop.texts$Col[loop.count], 0.1), loop.texts$Col[loop.count]), alpha = TRUE)(256), legend = FALSE, cex.axis = 0.8)
  text(250, y = 480, paste("p = ", round(quadrat.test(i, nx = 4, method = "MonteCarlo", CR = 0)$p.value, 3)), col = loop.texts$Col[loop.count], adj = c(1, 1))
  
  temp <- envelope(i, Kest, nsim = 199, nrank = 5, r = seq(0, 100, 5), verbose = FALSE, fix.n = TRUE, savefuns = TRUE, savepatterns = TRUE)
  plot(temp, cbind(obs - theo, lo = lo - theo, hi = hi - theo) ~ r, main = "K-r ~ r", legend = FALSE, cex.axis = 0.8, col = loop.texts$Col[loop.count], shadecol = alpha(loop.texts$Col[loop.count], 0.3))
  text(100, y = loop.texts$K[loop.count], paste("p = ", round(mad.test(temp, use.theory = TRUE)$p.value, 3)), adj = c(1, 1))
    abline(h = 0, lty = 2, col = loop.texts$Col[loop.count], lwd = 0.5)
    
  temp2 <- envelope(i, pcf, nsim = 199, nrank = 5, r = seq(0, 100, 5), verbose = FALSE, simulate = temp,  savefuns = TRUE)
  plot(temp2, cbind(obs - mmean + 1, lo = lo - mmean + 1, hi = hi - mmean + 1) ~ r, main = "pcf ~ r", legend = FALSE, cex.axis = 0.8, col = loop.texts$Col[loop.count], shadecol = alpha(loop.texts$Col[loop.count], 0.3))
  text(100, y = loop.texts$PCF[loop.count], paste("p = ", round(mad.test(temp2, use.theory = TRUE)$p.value, 3)), adj = c(1, 1))
  abline(h = 1, lty = 2, col = loop.texts$Col[loop.count], lwd = 0.5)
}

Figure 6: Estimated densities (left) and summary functions K (center) and pair correlation (right) as functions of the interaction radius (r) for the colonies of all and each of the Pogonomyrmex species. Envelopes show the 95% critical points for Monte Carlo pointwise tests estimated after 199 iterations of CSR spatial processes. P values are the results of MonteCarlo quadrat tests (left, 1999 iterations) and Maximum Absolute Deviation or ‘global’ tests (center, right) against CSR.

# figs of K/pcfcross
par(mar = c(2, 2, 1, .5))
layout(matrix(1:6, 2, 3, byrow = FALSE))

loop.texts <- data.frame(Main = c("PM / PI", "PM / PR", "PI / PR"), K = c(-30, -25, 40), PCF = c(1.3, 0.85, 0.85))
loop.count <- 0

for(i in list(c("PM", "PI"), c("PM", "PR"), c("PI", "PR"))) {
  loop.count <- loop.count + 1
  #test of independence through random labelling (= conditional on the global inhomogeneity)
  #ad-hoc function averages linearKcross ij and ji
  #toroidal-shift is not available because non-rectangular window, and other edge managements are problematic in this case
  temp <- envelope(xyNidos, Kcross.avrg, i[1], i[2], sigma = 25, r = seq(0, 100, 5), nsim = 199, nrank = 5, simulate = expression(rlabel(xyNidos)), verbose = FALSE, savefuns = TRUE, savepatterns = TRUE)
  plot(temp, cbind(obs - mmean, lo = lo - mmean, hi = hi - mmean) ~ r, main = paste0(loop.texts$Main[loop.count], ": Kcross"), legend = FALSE, cex.axis = .8)
  text(100, y = loop.texts$K[loop.count], paste("p = ", round(mad.test(temp, use.theory = FALSE)$p.value, 3)), adj = c(1, 0))
  abline(h = 0, lty = 2, col = "red", lwd = 0.5)
    
  temp2 <- envelope(xyNidos, pcfcross.avrg, i[1], i[2], sigma = 25, r = seq(0, 100, 5), simulate = temp, nsim = 199, nrank = 5, verbose = FALSE, savefuns = TRUE)
  plot(temp2, cbind(obs - mmean + 1, lo = lo - mmean + 1, hi = hi - mmean + 1) ~ r, main = paste0(loop.texts$Main[loop.count], ": pcfcross"), legend = FALSE, cex.axis = 0.8)
  text(100, y = loop.texts$PCF[loop.count], paste("p = ", round(mad.test(temp2, use.theory = FALSE)$p.value, 3)), adj = c(1, 0))
  abline(h = 1, lty = 2, col = "red", lwd = 0.5)
}

Figure 7: Averaged multitype K (above) and pair correlation (below) summary functions between each pair of Pogonomyrmex species (cross-type) as function of the interaction radius (r). Envelopes show the 95% critical points for Monte Carlo pointwise tests estimated after 199 iterations of CSRI spatial processes. Observed values below the envelopes indicate spatial segregation (fewer colonies of one species at r distance from a typical colony of the other species than expected by chance). P values are the results of MonteCarlo Maximum Absolute Deviation or ‘global’ tests against CSRI.

#segregation tests per Season [no PI in Season 1]
NidosS2.segr.test <- segregation.test(NidosS2, nsim = 199, sigma = 25, verbose = FALSE)
NidosS3.segr.test <- segregation.test(NidosS3, nsim = 199, sigma = 25, verbose = FALSE)

Despite the intense temporal dynamics of their colonies (appearance/disappearance, activation/deactivation, nest displacements), the areas occupied and spatial patterns of each Pogonomyrmex species mostly repeated among consecutive seasons, always showing aggregation (compared to CSRI) and occupying the same patches (Fig. 8). The three species proved also spatially segregated when only considering those colonies that were active within the same season (Monte Carlo segregation tests; Season 2: T = 18.422, p = 0.005; Season 3: T = 15.655, p = 0.005; no active P. inermis colonies were detected during Season 1).

par(mar = c(2, 2, 1, .5))
layout(matrix(1:9, 3, 3, byrow = FALSE), widths = c(1.5, 1, 1))

loop.graphs <- list(Main = c("Density 1999–2002: PI", "Density 1999–2002: PM", "Density 1999–2002: PR"), K = list(c(-15000, 75000), c(-15000, 75000), c(-15000, 60000)), PCF = list(c(0, 10), c(0, 10), c(0, 10)), Q.texts = list(c(480, 450, 420), c(480, 450, 420), c(480, 450, 420)), K.texts = list(c(75000, 67500, 60000), c(75000, 67500, 60000), c(60000, 52500, 45000)), PCF.texts = list(c(10, 9, 8), c(10, 9, 8), c(10, 9, 8)), Sp = c("PI", "PM", "PR"), Cols = c("#6600FF", "#009933", "#FF6600"))

#intensity graphs
for(i in 1:length(loop.graphs$Sp)) { #especies (plot)
  for(j in 1:length(NidosSeasons)) { #años (curvas)
    #condition to hack lack of PI colonies in Season1
    if(sum(NidosSeasons[[j]]$marks == loop.graphs$Sp[i]) < 1) {
      plot(density(unmark(subset(NidosSeasons[[j]], marks == loop.graphs$Sp[i], drop = TRUE)), sigma = 25, diggle = TRUE, eps = 1), main = loop.graphs$Main[i], col = alpha(loop.graphs$Cols[j], 0.025), ribbon = FALSE, cex.axis = 0.8, add = j > 1)
      next
    }
    plot(density(unmark(subset(NidosSeasons[[j]], marks == loop.graphs$Sp[i], drop = TRUE)), sigma = 25, diggle = TRUE, eps = 1), main = loop.graphs$Main[i], col = colorRampPalette(c(alpha(loop.graphs$Cols[j], 0.025), loop.graphs$Cols[j]), alpha = TRUE)(256), ribbon = FALSE, cex.axis = 0.8, add = j > 1)
    text(x = 250, y = loop.graphs$Q.texts[[i]][j], paste("p = ", round(quadrat.test(unmark(subset(NidosSeasons[[j]], marks == loop.graphs$Sp[i], drop = TRUE)), nx = 4, method = "MonteCarlo", CR = 0)$p.value, 3)), col = loop.graphs$Cols[j], adj = c(1, 1))
  }
}

#K graphs
for(i in 1:length(loop.graphs$Sp)) { #especies (plot)
  for(j in 1:length(NidosSeasons)) { #años (curvas)
    #condition to hack lack of PI colonies in Season1
    if(sum(NidosSeasons[[j]]$marks == loop.graphs$Sp[i]) < 1) {
      plot(x = 0:100, y = rep.int(0, 101), type = "n", ylim = loop.graphs$K[[i]], main = "K-r ~ r", cex.axis = 0.8)
      next
    }
    temp <- envelope(subset(NidosSeasons[[j]], marks == loop.graphs$Sp[i], drop = TRUE), Kest, nsim = 199, nrank = 5, r = seq(0, 100, 5), verbose = FALSE, fix.n = TRUE, savefuns = TRUE)
    plot(temp, cbind(obs - theo, lo = lo - theo, hi = hi - theo) ~ r, main = "K-r ~ r", legend = FALSE, cex.axis = 0.8, col = loop.graphs$Cols[j], shadecol = alpha(loop.graphs$Cols[j], 0.25), ylim = loop.graphs$K[[i]], add = j > 1)
    text(25, y = loop.graphs$K.texts[[i]][j], paste("p = ", round(mad.test(temp, use.theory = TRUE)$p.value, 3)), adj = c(1, 1), col = loop.graphs$Cols[j])
      abline(h = 0, lty = 2, col = "red", lwd = 0.5)
  }
}

#PCF graphs
for(i in 1:length(loop.graphs$Sp)) { #especies (plot)
  for(j in 1:length(NidosSeasons)) { #años (curvas)
    #condition to hack lack of PI colonies in Season1
    if(sum(NidosSeasons[[j]]$marks == loop.graphs$Sp[i]) < 1) {
      plot(x = 0:100, y = rep.int(0, 101), type = "n", ylim = loop.graphs$PCF[[i]], main = "pcf ~ r", cex.axis = 0.8)
      next
    }
    temp <- envelope(subset(NidosSeasons[[j]], marks == loop.graphs$Sp[i], drop = TRUE), pcf, nsim = 199, nrank = 5, r = seq(0, 100, 5), verbose = FALSE, fix.n = TRUE, savefuns = TRUE)
    plot(temp, main = "pcf ~ r", legend = FALSE, cex.axis = 0.8, col = loop.graphs$Cols[j], shadecol = alpha(loop.graphs$Cols[j], 0.25), ylim = loop.graphs$PCF[[i]], add = j > 1)
    text(100, y = loop.graphs$PCF.texts[[i]][j], paste("p = ", round(mad.test(temp, use.theory = TRUE)$p.value, 3)), adj = c(1, 1), col = loop.graphs$Cols[j])
      abline(h = 1, lty = 2, col = "red", lwd = 0.5)
  }
}

Figure 8: Same as Fig. 6, with a curve per activity season: violet: 1999-2000; green: 2001-2001; orange: 2001-2002.

Spatial association with dirt roads

Even when colonies on the dirt roads themselves were not part of this dataset (but see [Analysis of the other dataset]), all but three nests (97%) were <110 m from the closest dirt road (Table 1, Fig. 9A). Pogonomyrmex spp. nest density increased up to the first 50 m and then decreased markedly away from the roads, although the persistence of their external activity along the seasons was not so clearly related with distance (with a slightly negative but poorly explanatory tendency mostly associated with P. mendozanus; Fig. 9B).

#nidos por categoría de distancia
knitr::kable(table(Destino_nidos$Sp, Destino_nidos$Dist), caption = "Number of colonies of each *Pogonomyrmex* species at different distance categories from the closest dirt road. Distance intervals were defined for each to have a similar sampled area (~3 ha, <2% maximum difference) so the same number of colonies is expected in each interval under a completely random distribution.")

Table 1: Number of colonies of each Pogonomyrmex species at different distance categories from the closest dirt road. Distance intervals were defined for each to have a similar sampled area (~3 ha, <2% maximum difference) so the same number of colonies is expected in each interval under a completely random distribution.

	0–25 m	25–40 m	40–90 m	90–140 m	>140 m
PI	6	7	3	1	1
PM	6	9	10	2	0
PR	9	18	14	6	2

grid.arrange(
  #A. distancias nidos a picadas
  ggplot(Mapa_nidos, aes(x = distpic, fill = Sp)) +
    geom_histogram(bins = 20) +
    scale_fill_manual(name = "Species", values = c("darkgreen", "darkblue", "darkred")) +
    labs(x = "Distance to dirt roads (m)", y = "Number of colonies", tag = "A") +
    tema + theme(legend.position.inside = c(.85,.85)),
  #B. distancias nidos a picadas
  ggplot(Mapa_nidos, aes(x = distpic, y = Active, colour = Sp)) +
    geom_jitter(height = 0.2) +
    geom_smooth(method = "lm", formula = y ~ polynom(x, 2), se = FALSE) +
    geom_smooth(aes(colour = 1), colour = "black", method = "lm", formula = y ~ polynom(x, 2), se = FALSE) +
    scale_colour_manual(values = c("darkgreen", "darkblue", "darkred")) +
    labs(x = "Distance to dirt roads (m)", y = "Periods active", tag = "B") +
    tema,
  nrow = 1)

Figure 9: Pogonomyrmex colonies vs. distance to dirt roads. A. Number of detected Pogonomyrmex spp. colonies along three consecutive activity seasons (Oct–Apr 1999–2002) vs. distance to the closest of the two dirt roads in the area. B. Number of sampling periods (out of 12) in which each detected colony showed external activity vs distance to dirt roads (with best quadratic linear fit by species and in general [in black]).

Considering both their positions and the time they were active, the density of Pogonomyrmex colonies (irrespective of species) peaked at some short distance (< 50 m) from the dirt roads. When split by species, the decline of colonies with distance to dirt roads seems strong in P. inermis, starts farther away for P. mendozanus, with the distribution of P. rastratus colonies more uniform within the first 100 m from the roads and strongly declining farther (Fig. 10).

par(mar = c(4, 4, 2, 1))
#intensidad en función de distancia (x sp), pesado por la actividad externa de las colonias
plot(rhohat(xyNidos, DistPicadas, bw = 25, weights = Mapa_nidos$Active), col = "black", shadecol = alpha("black", 0.25), ylim = c(0, 0.004), ylab = "Lambda", xlab = "Distance (m)", legend = FALSE, main = "Distance to dirt roads (weighted by activity)", cex.axis = .8)
text(200, y = 0.004, paste("p = ", round(cdf.test(unmark(xyNidos), DistPicadas, test = "ks", jitter = F)$p.value, 3)), col = "black", adj = c(0, 1))
plot(rhohat(xyNidos[xyNidos$marks == "PR"], DistPicadas, bw = 25, weights = Mapa_nidos$Active[Mapa_nidos$Sp == "PR"]), col = "darkred", shadecol = alpha("darkred", 0.25), cex.axis = .8, add = TRUE)
text(200, y = 0.0038, paste("p = ", round(cdf.test(unmark(xyNidos[xyNidos$marks == "PR"]), DistPicadas, test = "ks", jitter = F)$p.value, 3)), col = "darkred", adj = c(0, 1))
plot(rhohat(xyNidos[xyNidos$marks == "PM"], DistPicadas, bw = 25, weights = Mapa_nidos$Active[Mapa_nidos$Sp == "PM"]), col = "darkblue", shadecol = alpha("darkblue", 0.25), add = TRUE)
text(200, y = 0.0036, paste("p = ", round(cdf.test(unmark(xyNidos[xyNidos$marks == "PM"]), DistPicadas, test = "ks", jitter = F)$p.value, 3)), col = "darkblue", adj = c(0, 1))
plot(rhohat(xyNidos[xyNidos$marks == "PI"], DistPicadas, bw = 25, weights = Mapa_nidos$Active[Mapa_nidos$Sp == "PI"]), col = "darkgreen", shadecol = alpha("darkgreen", 0.25), add = TRUE)
text(200, y = 0.0034, paste("p = ", round(cdf.test(unmark(xyNidos[xyNidos$marks == "PI"]), DistPicadas, test = "ks", jitter = F)$p.value, 3)), col = "darkgreen", adj = c(0, 1))

Figure 10: Nonparametric estimate of the intensity of Poisson spatial point processes for all Pogonomyrmex colonies (black) and for colonies of each species (green: PI; blue: PM; red: PR), weighted by the number of sampling seasons they were active (0–12), as function of the distance to the closest dirt road (m). P values are for Kolmogorov Smirnov tests of goodness-of-fit of observed vs predicted (CSR = uniform) distributions of distances. Marks on the horizontal axis show the actual distance to the nearest road of each colony independently of species.

In fact, distance to dirt roads is a strong predictor of the density of Pogonomyrmex colonies within the algarrobal in log-quadratic Poisson spatial models, with the shape of its influence varying per species. Again, P. inermis decreases with distance, while P. mendozanus and P. rastratus increase in the first 50 m from the closest dirt road and then abruptly decline (Fig. 11).

# Homogeneous Poisson
xyNidos.Poisson <- ppm(xyNidos ~ marks)

# Poisson model with common cuadratic curve or a curve per species
xyNidos.PoisDist <- ppm(xyNidos ~ marks + polynom(DistPicadas, 2))
xyNidos.PoissonDist <- ppm(xyNidos ~ marks * polynom(DistPicadas, 2))

# comparing shapes [cubic is significative but only allows to fit the single far away PI nest, overfitting?]
anova(xyNidos.Poisson,
      ppm(xyNidos ~ marks + DistPicadas),
      xyNidos.PoisDist,
      test = "LR")

## Analysis of Deviance Table
## 
## Model 1: ~marks   Poisson
## Model 2: ~marks + DistPicadas     Poisson
## Model 3: ~marks + (DistPicadas + I(DistPicadas^2))    Poisson
##   Npar Df Deviance   Pr(>Chi)    
## 1    3                           
## 2    4  1   23.545 0.00000122 ***
## 3    5  1   16.186 0.00005741 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#comparing no curve, single curve vs curve per species
anova(xyNidos.Poisson,
      xyNidos.PoisDist,
      xyNidos.PoissonDist,
      test = "LR")

## Analysis of Deviance Table
## 
## Model 1: ~marks   Poisson
## Model 2: ~marks + (DistPicadas + I(DistPicadas^2))    Poisson
## Model 3: ~marks * (DistPicadas + I(DistPicadas^2))    Poisson
##   Npar Df Deviance       Pr(>Chi)    
## 1    3                               
## 2    5  2   39.731 0.000000002357 ***
## 3    9  4   14.625       0.005547 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

layout(matrix(1))
plot(effectfun(xyNidos.PoissonDist, "DistPicadas", marks = "PI", se.fit = TRUE), col = "darkgreen", shadecol = alpha("darkgreen", 0.2), legend = FALSE, main = "Model: log(λ) ~ Sp * polynom(Dist, 2)", ylim = c(0, 0.00085), xlab = "Distance to dirt road (m)", ylab = "Intensity (λ)")
plot(effectfun(xyNidos.PoissonDist, "DistPicadas", marks = "PM", se.fit = TRUE), col = "darkblue", shadecol = alpha("darkblue", 0.2), legend = FALSE, add = TRUE)
plot(effectfun(xyNidos.PoissonDist, "DistPicadas", marks = "PR", se.fit = TRUE), col = "darkred", shadecol = alpha("darkred", 0.2), legend = FALSE, add = TRUE)

Figure 11: Fitted species-specific relationship of colony density with distance to closest dirt road in the log-cuadratic Poisson spatial model y ~ spp * polynom(Distance, 2).

However the residuals of this Poisson model are still strongly aggregated (see Fig. 14 below), indicating that this single predictor cannot fully explain the heterogeneous density of Pogonomyrmex colonies in the area, or that the assumption of spatial independence among colonies/species is not correct. In fact, maps (Figs. 4 and 5) clearly show that colonies of the same species aggregate within particular sectors close to dirt roads, leaving others (equally close) unoccupied by one or more species. The same happens with Poisson models based on only those colonies active in the same season (for Seasons 2 and 3 when all species were present). Evidently there are other factors involved that result in the aggregation or absence of Pogonomyrmex colonies (e.g., along the dirt roads).

Models with intercolony interactions

Lacking information on other covariables, one way to model those aggregations is allowing for positive/negative interactions among the colonies of the same or different species. Gibbs point pattern processes incorporate negative (leading to inhibitory, overdispersed or repulsive patterns) or positive (moderately aggregated patterns) interactions that affect the probability of an event given the distance to other, neighbouring events.

The interaction profiles among neighbouring colonies show a minimum distance of ~4 m between closest neighbours and positive interactions (leading to aggregation) up to a around 15 m (either under CSR or the Poisson model that includes distance to dirt roads as covariable; Fig. 12.

par(mar = c(2, 2, 1, .5))
layout(matrix(1:8, 2, 4, byrow = TRUE))
plot(fitin(ppm(xyNidos ~ marks * polynom(DistPicadas, 2), PairPiece(seq(4, 40, 2)), correction = "isotropic")), main = "Pogonomyrmex spp (full model)", legend = FALSE)
plot(fitin(ppm(unmark(xyNidos[xyNidos$marks == "PI"]) ~ polynom(DistPicadas, 2), PairPiece(seq(4, 40, 2)), correction = "isotropic")), main = "PI ~ DistPicadas", col = "darkgreen", legend = FALSE)
plot(fitin(ppm(unmark(xyNidos[xyNidos$marks == "PM"]) ~ polynom(DistPicadas, 2), PairPiece(seq(4, 40, 2)), correction = "isotropic")), main = "PM ~ DistPicadas", col = "darkblue", legend = FALSE)
plot(fitin(ppm(unmark(xyNidos[xyNidos$marks == "PR"]) ~ polynom(DistPicadas, 2), PairPiece(seq(4, 40, 2)), correction = "isotropic")), main = "PR ~ DistPicadas", col = "darkred", legend = FALSE)

plot(profilepl(data.frame(r = seq(4, 50, 2), hc = 3.99), StraussHard, xyNidos ~ marks * polynom(DistPicadas, 2), correction = "isotropic", AIC = TRUE, fast = FALSE, verbose = FALSE), main = "Strauss+Hard, full model")
plot(profilepl(data.frame(r = seq(4, 50, 2), hc = 3.99), StraussHard, unmark(xyNidos[xyNidos$marks == "PI"]) ~ polynom(DistPicadas, 2), correction = "isotropic", AIC = TRUE, fast = FALSE, verbose = FALSE), main = "Strauss+Hard, PI", col = "darkgreen")
plot(profilepl(data.frame(r = seq(4, 50, 2), hc = 3.99), StraussHard, unmark(xyNidos[xyNidos$marks == "PM"]) ~ polynom(DistPicadas, 2), correction = "isotropic", AIC = TRUE, fast = FALSE, verbose = FALSE), main = "Strauss+Hard, PM", col = "darkblue")
plot(profilepl(data.frame(r = seq(4, 50, 2), hc = 3.99), StraussHard, unmark(xyNidos[xyNidos$marks == "PR"]) ~ polynom(DistPicadas, 2), correction = "isotropic", AIC = TRUE, fast = FALSE, verbose = FALSE), main = "Strauss+Hard, PR", col = "darkred")

Figure 12: Above: Between-colonies interactions (for all colonies or separated by species) as function of their interpoint distance in m (discounting the log-quadratic Poisson model fitted). Values >1 indicate atraction/aggregation, with more colonies separated by that distance interval than expected by chance in a Poisson model; <1 indicates repulsion/segregation. Below: Optimal values for the maximum limit of interpoint interactions in Strauss+Hard models (with hard limit fixed at hc = 3.99 m).

A StraussHard model requiring at least 4 m between Pogonomyrmex colonies and favouring their co-ocurrence up to 15 m (irrespective of species identity) proves a much better fit to the observed spatial pattern than the equivalent Poisson model assuming no interactions, without changing the general shape of the species-specific effect of distance to dirt roads (except to allow for the single but very influential P. inermis colony appearing in Season 3 216.73 m from the closest dirt road; Fig. 13).

#single-parameter StraussHard model
(xyNidos.SHDist <- ppm(xyNidos ~ marks * polynom(DistPicadas, 2), interaction = StraussHard(h = 4, r = 15)))

## Nonstationary Strauss - hard core process
## Fitted to point pattern dataset 'xyNidos'
## 
## Possible marks: 'PI', 'PM' and 'PR'
## 
## Log trend:  ~marks * (DistPicadas + I(DistPicadas^2))
## 
## Fitted trend coefficients:
##              (Intercept)                  marksPM                  marksPR 
##            -8.1801352844            -1.0688098500            -1.0795768372 
##              DistPicadas         I(DistPicadas^2)      marksPM:DistPicadas 
##            -0.0379573857             0.0001576698             0.0631461019 
##      marksPR:DistPicadas marksPM:I(DistPicadas^2) marksPR:I(DistPicadas^2) 
##             0.0754926267            -0.0004802846            -0.0004731118 
## 
## Interaction distance:    15
## Hard core distance:  4
## Fitted interaction parameter gamma:  2.1747596
## 
## Relevant coefficients:
## Interaction 
##   0.7769181 
## 
## For standard errors, type coef(summary(x))

#comparing Poisson vs single-parameter StraussHard
anova(xyNidos.PoissonDist,
      xyNidos.SHDist,
      test = "LR")

## Analysis of Deviance Table
## 
## Model 1: ~marks * (DistPicadas + I(DistPicadas^2))    Poisson
## Model 2: ~marks * (DistPicadas + I(DistPicadas^2))    StraussHard
##   Npar Df AdjDeviance  Pr(>Chi)    
## 1  236                             
## 2  237  1      176.06 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

layout(matrix(1))
plot(effectfun(xyNidos.SHDist, "DistPicadas", marks = "PI", se.fit = TRUE), col = "darkgreen", shadecol = alpha("darkgreen", 0.2), legend = FALSE, main = "Model: log(λ) ~ Sp * polynom(Dist, 2) + StraussHard(single r)", ylim = c(0, 0.00045))
plot(effectfun(xyNidos.SHDist, "DistPicadas", marks = "PM", se.fit = TRUE), col = "darkblue", shadecol = alpha("darkblue", 0.2), legend = FALSE, add = TRUE)
plot(effectfun(xyNidos.SHDist, "DistPicadas", marks = "PR", se.fit = TRUE), col = "darkred", shadecol = alpha("darkred", 0.2), legend = FALSE, add = TRUE)

Figure 13: Fitted species-specific relationship of colony density with distance to closest dirt road in a log-cuadratic Strauss+Hard spatial model: y ~ spp * polynom(Distance, 2) with hc = 4 and r = 15.

This Strauss+Hard model, that forbids neighbouring colonies <4 m apart but doubles the chance of colonies co-ocurring at 4–15 m distances, also substantially corrects for the spatial aggregation of the model residuals (at least up to distances ~60 m; Fig. 14).

layout(matrix(1:4, 1, 4))
plot(Kres(xyNidos.Poisson, r = seq(0, 100, 2)), ires ~ r, shade = c("ihi", "ilo"), main = "Kres (Poisson - null)", legend = FALSE)
plot(Kres(xyNidos.PoisDist, r = seq(0, 100, 2)), ires ~ r, shade = c("ihi", "ilo"), main = "Kres (Poisson - spp + Dist)", legend = FALSE)
plot(Kres(xyNidos.PoissonDist, r = seq(0, 100, 2)), ires ~ r, shade = c("ihi", "ilo"), main = "Kres (Poisson - spp * Dist)", legend = FALSE)
plot(Kres(xyNidos.SHDist, r = seq(0, 100, 2)), ires ~ r, shade = c("ihi", "ilo"), main = "Kres (single StraussHard - spp * Dist)", legend = FALSE)

Figure 14: Residual K function of a null Poisson model (left, assuming homogeneous intensity and no interpoint interactions), a log-quadratic Poisson model dependent on distance to the closest dirt road (center, no interactions) and a similar Strauss+Hard model allowing for aggregation between Pogonomyrmex colonies (over that expected by heterogeneous intensity of the point process associated with distance to the dirt roads).

When allowing for different (species-specific) ranges of interaction, maximum aggregation associates with smaller, more compact clusters for P. inermis than for P. rastratus and P. mendozanus (Fig. 12 and results below).

A Strauss+Hard model allowing for different radii of intraspecific colony interactions (assuming independence among species) is a better fit than both the single-parameter model (above) and a model with only interspecific interactions (assuming intraspecific independence). The more complex model allowing for the estimation of both intra- and interspecific interaction parameters fits the data only slightly better than the only-intraspecific model.

In spite of the different combinations, intraspecific interactions are always estimated positive (aggregation), and interspecific interactions depend on species identities: they are positive when involving PM, but negative (repulsion) between PI and PR (Fig. 15).

# profile para determinar los radios de interacción intra de c/sp; cerca de los valores vistos arriba para cada sp por separado
profilepl(expand.grid(r1 = seq(12, 16, 1), r2 = seq(36, 46, 1), r3 = seq(30, 34, 1)), function(r1, r2, r3) {MultiStraussHard(iradii = diag(c(r1, r2, r3)), hradii = diag(c(3.99, 3.99, 3.99)))}, xyNidos ~ marks * polynom(DistPicadas, 2), correction = "isotropic", AIC = TRUE, fast = FALSE, verbose = FALSE)

## profile log pseudolikelihood for model:
## ppm(xyNidos ~ marks * polynom(DistPicadas,  2),  correction = "isotropic",      
##  interaction = function(r1,  r2,  r3) {         MultiStraussHard(iradii = 
## diag(c(r1,  r2,  r3)),  hradii = diag(c(3.99,               3.99,  3.99)))     
## },  AIC = TRUE,  fast = FALSE)
## interaction: Multitype Strauss Hardcore process
## irregular parameters:
##  r1 in [12, 16]
## r2 in [36, 46]
## r3 in [30, 34]
## optimum values of irregular parameters:
##  r1 = 15, r2 = 44 and r3 = 32

# only-intraspecific model
(xyNidos.SHDist.intrasp <- ppm(xyNidos ~ marks * polynom(DistPicadas, 2), interaction = MultiStraussHard(iradii = diag(c(15, 44, 32)), hradii = diag(c(4, 4, 4)))))

## Nonstationary Multitype Strauss Hardcore process
## Fitted to point pattern dataset 'xyNidos'
## 
## Possible marks: 'PI', 'PM' and 'PR'
## 
## Log trend:  ~marks * (DistPicadas + I(DistPicadas^2))
## 
## Fitted trend coefficients:
##              (Intercept)                  marksPM                  marksPR 
##         -10.956420058103          -3.938047262285          -0.036729808593 
##              DistPicadas         I(DistPicadas^2)      marksPM:DistPicadas 
##           0.009583113998          -0.000004941612           0.111773618678 
##      marksPR:DistPicadas marksPM:I(DistPicadas^2) marksPR:I(DistPicadas^2) 
##           0.045162849665          -0.000868466765          -0.000322987544 
## 
## 3 types of points
## Possible types:
## [1] PI PM PR
## Interaction radii:
##    PI PM PR
## PI 15 NA NA
## PM NA 44 NA
## PR NA NA 32
## Hardcore radii:
##    PI PM PR
## PI  4 NA NA
## PM NA  4 NA
## PR NA NA  4
## Fitted interaction parameters gamma_ij
##          PI       PM       PR
## PI 5.680125       NA       NA
## PM       NA 1.550826       NA
## PR       NA       NA 1.475766
## 
## Relevant coefficients:
## markPIxPI markPMxPM markPRxPR 
## 1.7369733 0.4387877 0.3891770 
## 
## For standard errors, type coef(summary(x))

#profile for the 3 interspecies parameters assuming no intraspecific interactions [after several iterative runs to precise the range and step of candidate values]
profilepl(expand.grid(rIM = seq(12, 18, 1), rIR = seq(20, 24, 1), rMR = seq(22, 26, 1)), function(rIM, rIR, rMR) {MultiStraussHard(iradii = matrix(c(0, rIM, rIR, rIM, 0, rMR, rIR, rMR, 0), nrow = 3), hradii = matrix(rep(3.99, 9), nrow = 3))}, xyNidos ~ marks * polynom(DistPicadas, 2), correction = "isotropic", AIC = TRUE, fast = FALSE, verbose = FALSE)

## profile log pseudolikelihood for model:
## ppm(xyNidos ~ marks * polynom(DistPicadas,  2),  correction = "isotropic",      
##  interaction = function(rIM,  rIR,  rMR) {         MultiStraussHard(iradii = 
## matrix(c(0,  rIM,  rIR,  rIM,               0,  rMR,  rIR,  rMR,  0),  nrow = 
## 3),  hradii = matrix(rep(3.99,               9),  nrow = 3))     },  AIC = 
## TRUE,  fast = FALSE)
## interaction: Multitype Strauss Hardcore process
## irregular parameters:
##  rIM in [12, 18]
## rIR in [20, 24]
## rMR in [22, 26]
## optimum values of irregular parameters:
##  rIM = 16, rIR = 22 and rMR = 24

# only-interspecific model
(xyNidos.SHDist.intersp <- ppm(xyNidos ~ marks * polynom(DistPicadas, 2), interaction = MultiStraussHard(iradii = matrix(c(0, 16, 22, 16, 0, 24, 22, 24, 0), nrow = 3), hradii = matrix(rep(4, 9), nrow = 3))))

## Nonstationary Multitype Strauss Hardcore process
## Fitted to point pattern dataset 'xyNidos'
## 
## Possible marks: 'PI', 'PM' and 'PR'
## 
## Log trend:  ~marks * (DistPicadas + I(DistPicadas^2))
## 
## Fitted trend coefficients:
##              (Intercept)                  marksPM                  marksPR 
##            -8.1626810671            -3.9244506505            -0.5887593151 
##              DistPicadas         I(DistPicadas^2)      marksPM:DistPicadas 
##            -0.0268373693             0.0001089175             0.1479890176 
##      marksPR:DistPicadas marksPM:I(DistPicadas^2) marksPR:I(DistPicadas^2) 
##             0.0729293150            -0.0010970499            -0.0005139456 
## 
## 3 types of points
## Possible types:
## [1] PI PM PR
## Interaction radii:
##    PI PM PR
## PI NA 16 22
## PM 16 NA 24
## PR 22 24 NA
## Hardcore radii:
##    PI PM PR
## PI  4  4  4
## PM  4  4  4
## PR  4  4  4
## Fitted interaction parameters gamma_ij
##           PI       PM        PR
## PI        NA 3.118408 0.2438349
## PM 3.1184077       NA 1.5011648
## PR 0.2438349 1.501165        NA
## 
## Relevant coefficients:
##  markPIxPM  markPIxPR  markPMxPR 
##  1.1373225 -1.4112641  0.4062413 
## 
## For standard errors, type coef(summary(x))

#profile exploring all 6 parameters 3 units up and down around their previously selected values [to check if they move or behave as independent; increasing candidate values simultaneously soon becomes prohibitive (n^6 models)!] -> apparently they behave as independent [matching previously selected radii]
profilepl(expand.grid(r1 = seq(12, 18, 3), r2 = seq(41, 47, 3), r3 = seq(29, 35, 3), rIM = seq(12, 18, 3), rIR = seq(19, 25, 3), rMR = seq(21, 27, 3)), function(r1, r2, r3, rIM, rIR, rMR) {MultiStraussHard(iradii = matrix(c(r1, rIM, rIR, rIM, r2, rMR, rIR, rMR, r3), nrow = 3), hradii = matrix(rep(3.99, 9), nrow = 3))}, xyNidos ~ marks * polynom(DistPicadas, 2), correction = "isotropic", AIC = TRUE, fast = FALSE, verbose = FALSE)

## profile log pseudolikelihood for model:
## ppm(xyNidos ~ marks * polynom(DistPicadas,  2),  correction = "isotropic",      
##  interaction = function(r1,  r2,  r3,  rIM,  rIR,  rMR) {         
## MultiStraussHard(iradii = matrix(c(r1,  rIM,  rIR,  rIM,               r2,  
## rMR,  rIR,  rMR,  r3),  nrow = 3),  hradii = matrix(rep(3.99,               9), 
##  nrow = 3))     },  AIC = TRUE,  fast = FALSE)
## interaction: Multitype Strauss Hardcore process
## irregular parameters:
##  r1 in [12, 18]
## r2 in [41, 47]
## r3 in [29, 35]
## rIM in [12, 18]
## rIR in [19, 25]
## rMR in [21, 27]
## optimum values of irregular parameters:
##  r1 = 15, r2 = 44, r3 = 32, rIM = 15, rIR = 22 and rMR = 24

# both intra- and interspecific model (full)
(xyNidos.SHDist.fullsp <- ppm(xyNidos ~ marks * polynom(DistPicadas, 2), interaction = MultiStraussHard(iradii = matrix(c(15, 15, 22, 15, 44, 24, 22, 24, 32), nrow = 3), hradii = matrix(rep(4, 9), nrow = 3))))

## Nonstationary Multitype Strauss Hardcore process
## Fitted to point pattern dataset 'xyNidos'
## 
## Possible marks: 'PI', 'PM' and 'PR'
## 
## Log trend:  ~marks * (DistPicadas + I(DistPicadas^2))
## 
## Fitted trend coefficients:
##              (Intercept)                  marksPM                  marksPR 
##          -10.29431590039           -2.65755802360           -1.12502453847 
##              DistPicadas         I(DistPicadas^2)      marksPM:DistPicadas 
##           -0.00402045744            0.00005024147            0.05790604694 
##      marksPR:DistPicadas marksPM:I(DistPicadas^2) marksPR:I(DistPicadas^2) 
##            0.06201512594           -0.00044601716           -0.00037323648 
## 
## 3 types of points
## Possible types:
## [1] PI PM PR
## Interaction radii:
##    PI PM PR
## PI 15 15 22
## PM 15 44 24
## PR 22 24 32
## Hardcore radii:
##    PI PM PR
## PI  4  4  4
## PM  4  4  4
## PR  4  4  4
## Fitted interaction parameters gamma_ij
##          PI       PM       PR
## PI 4.970279 2.483445 0.384900
## PM 2.483445 1.465238 1.367639
## PR 0.384900 1.367639 1.471754
## 
## Relevant coefficients:
##  markPIxPI  markPIxPM  markPMxPM  markPIxPR  markPMxPR  markPRxPR 
##  1.6034760  0.9096469  0.3820176 -0.9547716  0.3130857  0.3864547 
## 
## For standard errors, type coef(summary(x))

#SH.fullsp wo. DistPicadas as predictor
xyNidos.SH.fullsp <- ppm(xyNidos ~ marks, interaction = MultiStraussHard(iradii = matrix(c(15, 15, 22, 15, 44, 24, 22, 24, 32), nrow = 3), hradii = matrix(rep(4, 9), nrow = 3)))

#comparing Model 1: single-parameter, 2: interspp, 3: intraspp, 4: both (full) StraussHard
anova(xyNidos.SHDist,
      xyNidos.SHDist.intersp,
      xyNidos.SHDist.intrasp,
      xyNidos.SHDist.fullsp,
      test = "LR")

## Analysis of Deviance Table
## 
## Model 1: ~marks * (DistPicadas + I(DistPicadas^2))    StraussHard
## Model 2: ~marks * (DistPicadas + I(DistPicadas^2))    MultiStraussHard
## Model 3: ~marks * (DistPicadas + I(DistPicadas^2))    MultiStraussHard
## Model 4: ~marks * (DistPicadas + I(DistPicadas^2))    MultiStraussHard
##   Npar Df AdjDeviance       Pr(>Chi)    
## 1   10                                  
## 2   12  2     -37.528 0.000000007093 ***
## 3   12  0      62.844      < 2.2e-16 ***
## 4   15  3       8.670        0.03402 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

par(mar = c(3.7, 4, .5, .5))
plot(fitin(xyNidos.SHDist.fullsp), title = "Strauss+Hard model w/intra+interspp interactions")

Figure 15: Intra and interspecific interactions as estimated with a log-quadratic Strauss+Hard model with fixed hard limit (hc = 4 m) but different interaction radii (r). Values of h > 1 imply aggregation (more frequent than expected by a Poisson model), and h < 1 repulsion. Estimated values for model with only-intraspp interactions (i.e., only values on diagonal) or with only-interspp interactions (i.e., only off-diagonal) were very similar (see fitted models above).

# similar estimated values in only-intra and only-interspp models [not shown]
# plot(fitin(xyNidos.SHDist.intrasp), title = "only intraspp")
# plot(fitin(xyNidos.SHDist.intersp), xlim = c(0, 55), title = "only interspp")

#SH.fullsp wo. DistPicadas as predictor
xyNidos.SH.fullsp <- ppm(xyNidos ~ marks, interaction = MultiStraussHard(iradii = matrix(c(15, 15, 22, 15, 44, 24, 22, 24, 32), nrow = 3), hradii = matrix(rep(4, 9), nrow = 3)))

The residuals of both models including intraspecific aggregation are no longer spatially aggregated (Fig. 16), so they seem enough to model the observed spatial pattern of active colonies along three consecutive seasons.

layout(matrix(1:6, 2, 3, byrow = TRUE))
plot(Kres(xyNidos.PoissonDist, r = seq(0, 100, 2)), ires ~ r, shade = c("ihi", "ilo"), main = "Kres (Poisson: no interactions)", legend = FALSE)
plot(Kres(xyNidos.SHDist, r = seq(0, 100, 2)), ires ~ r, shade = c("ihi", "ilo"), main = "Kres (S+H: any Pogo)", legend = FALSE)
plot(Kres(xyNidos.SHDist.intrasp, r = seq(0, 100, 2)), ires ~ r, shade = c("ihi", "ilo"), main = "Kres (S+H: only intraspp)", legend = FALSE)
plot(Kres(xyNidos.SHDist.intersp, r = seq(0, 100, 2)), ires ~ r, shade = c("ihi", "ilo"), main = "Kres (S+H: only interspp)", legend = FALSE)
plot(Kres(xyNidos.SHDist.fullsp, r = seq(0, 100, 2)), ires ~ r, shade = c("ihi", "ilo"), main = "Kres (S+H: both interactions)", legend = FALSE)
plot(Kres(xyNidos.SH.fullsp, r = seq(0, 100, 2)), ires ~ r, shade = c("ihi", "ilo"), main = "Kres (S+H: homogeneous, both interactions)", legend = FALSE)

Figure 16: Residual K function of log-quadraticPoisson model and Strauss+Hard models with a single interaction parameter for any Pogonomyrmex colonies, estimating only intraspecific interactions, only interspecific interactions, both kinds of interactions, and with both interactions but homogeneous intensity (cf. Fig. 14).

However, when the intraspecific aggregation of colonies is estimated explicitly in these Strauss+Hard models, the (additional) predictive power of the distance to dirt roads is strongly reduced.

This means that there is a strong intra-specific aggregation of the positions of Pogonomyrmex active colonies detected during three consecutive activity seasons that relates but cannot be fully explained by Poisson models with distance to dirt roads as the single predictor. Other environmental or population-level processes result in (additional) intra-specific colony aggregation.

#comparing SH models with/out DistPicadas as predictor
anova(xyNidos.SH.fullsp, 
      update(xyNidos.SHDist.fullsp, xyNidos ~ marks + polynom(DistPicadas, 2)),
      xyNidos.SHDist.fullsp, 
      test = "LR")

## Analysis of Deviance Table
## 
## Model 1: ~marks   MultiStraussHard
## Model 2: ~marks + DistPicadas + I(DistPicadas^2)      MultiStraussHard
## Model 3: ~marks * (DistPicadas + I(DistPicadas^2))    MultiStraussHard
##   Npar Df AdjDeviance Pr(>Chi)
## 1    9                        
## 2   11  2     0.36415   0.8335
## 3   15  4     2.45260   0.6531

A summary of the main Poisson and Gibbs fitted models follows (“Residuals”: number of signs (1-5) show the degree of spatial aggregation (+) or repulsion (-) of model residuals, with “0” indicating no significant spatial autocorrelation according to Ripley’s K function; see Figs. 14 and 16):

Model	Structure	Interactions	Intensities/Trends	Residuals	AIC
PS-NULL	spp	—	βPI: 0.00012 βPM: 0.00019 βPR: 0.00034	+++++	1781.83
PS-DIST	spp+dist2	—	βPI: 0.00011+1.02333m+0.99977m² βPM: 0.00017+1.02333m+0.99977m² βPR: 0.00031+1.02333m+0.99977m²	++++	1746.1
PS-SP×D	spp*dist2	—	βPI: 0.0003+0.98002m+1.00006m² βPM: 0.00006+1.0793m+0.99928m² βPR: 0.00021+1.03783m+0.99967m²	++++	1739.48
SH-SP×D-POGO	spp*dist2	ɣ: 2.175 [15]	βPI: 0.00028+0.96275m+1.00016m² βPM: 0.0001+1.02551m+0.99968m² βPR: 0.0001+1.03825m+0.99968m²	+	1418.18
SH-SP×D-INTRA	spp*dist2	ɣPI: 5.68 [15] ɣPM: 1.551 [44] ɣPR: 1.476 [32]	βPI: 0.00002+1.00963m+1m² βPM: 0+1.12903m+0.99913m² βPR: 0.00002+1.05627m+0.99967m²	0	755.42
SH-SP×D-INTER	spp*dist2	ɣPI-PM: 3.118 [16] ɣPI-PR: 0.244 [22] ɣPM-PR: 1.501 [24]	βPI: 0.00029+0.97352m+1.00011m² βPM: 0.00001+1.1288m+0.99901m² βPR: 0.00016+1.04717m+0.9996m²	++	1251.66
SH-SP×D-BOTH	spp*dist2	ɣPI: 4.97 [15] ɣPM: 1.465 [44] ɣPR: 1.472 [32] ɣPI-PM: 2.483 [15] ɣPI-PR: 0.385 [22] ɣPM-PR: 1.368 [24]	βPI: 0.00003+0.99599m+1.00005m² βPM: 0+1.05536m+0.9996m² βPR: 0.00001+1.05971m+0.99968m²	0	731.13
SH-NULL-BOTH	spp	ɣPI: 4.507 [15] ɣPM: 1.519 [44] ɣPR: 1.501 [32] ɣPI-PM: 2.458 [15] ɣPI-PR: 0.393 [22] ɣPM-PR: 1.386 [24]	βPI: 0.00004 βPM: 0.00001 βPR: 0.00009	0	723.7

Conclusions

• The number of active Pogonomyrmex colonies within the algarrobal varies within and between consecutive activity seasons, appearing and disappearing quite frequently. Moreover, they also move their nest entrances frequently.
  
• There is considerable differences in the abundance of each species: P. rastratus accounted for almost half of the colonies of harvester ants in this ~15 ha area.
  
• There are very heterogeneously distributed: almost half the area sampled had no active colonies during those three years.
  
• All three species are spatially aggregated, but while two of them (P. mendozanus and P. inermis) occurred in patches, P. rastratus was a bit more extended. Species are spatially segregated, particularly P. rastratus with P. inermis.
  
• The spatial patterns of the accumulated positions where active colonies occurred along three activity seasons also appear when only those colonies active in the same season were analysed.
  
• Distance to dirt roads is a strong predictor of the heterogeneous density of Pogonomyrmex colonies. Even when colonies on the roads themselves were not considered here, sectors close (<50 m) to dirt roads concentrate almost all colonies.
  
• Still, there are other environmental (e.g., soil properties, seed abundance) or population-level processes (e.g., restricted dispersal) that result in aggregation of the active colonies of the same species, and segregation between some species, that is stronger than what can be explained only by distance to dirt roads.