Pogos en picadas

---

Pogonomyrmex spp. on picadas (dirt roads) [Ñacuñán 2004–2010]

Summary of objectives and methods

Census in three summers (2004, 2007, 2010) of active nests of the three Pogonomyrmex species along 1.5 km segments of two dirt roads (La Fadina and Doble), plus a single census in a third dirt road (2006: de la Torre) [colonies were not marked permanently so each search is “random” every year; there may be minor variations in position for the same colonies even if they didn’t really move, particularly in Doble which had an arbitrarily marked starting point] PR = Pogonomyrmex rastratus (now P. propinqua [red]
PM = Pogonomyrmex mendozanus (ex P. pronotalis) [blue]
PI = Pogonomyrmex inermis [green]

Quick field measurements of permeability of the soil surface (infiltration time) and topography (msnm) along both segments (2010)

Soil texture (percentages of sand, silt and clay) in soil samples (0-30 cm) associated with sectors (“spatial clusters”) of active colonies of two of the species (PI and PM) along the roads Goals for 2004 census: ¿is there any inter/intraspecific spatial pattern? ¿are they associated with edaphic factors?

Goal 2006: spatial consistency: ¿do detected patterns appear elsewhere (in a different dirt road)?

Goals 2007-2010: temporal consistency and drought effect (dry period between 2007 and 2010) -> changes in abundance? (active colonies, species-specific effects), do spatial patterns and environmental associations persist through environmental change? (random loss? spatial aggregation towards core areas? density-dependent thinning? changes in intra/interspecific interactions/associations?)

Packages, databases and extra functions

#Esta es una versión de "Nidos en picadas.Rmd" para analizar los patrones y asociaciones espaciales con spatstat (usando redes de un solo segmento) en lugar de las funciones propias desarrolladas en "Ripley 1D.R"
#Algunas líneas se marcan temporalmente como comentarios para evaluar su necesidad/conveniencia al modificar el código

library(dplyr)
library(tidyr)
library(ggplot2)
library(gridExtra)
library(scales)
library(ggtern) # for ternary plot (soil texture)
library(glmmTMB) # for mixed models with beta distribution (soil texture)
library(DHARMa) # to test beta mixed model assumptions
library(spatstat) # documented in "Spatial Point Patterns: Methodology and Applications with R" by Baddeley et al (2016)
library(knitr) #for kable function -> nice tables in html

# tema general para ggplot en B&N sin grilla ni legend
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 y manejo de datos
# importación de planilla de nidos en picadas
Pogos_picadas <- read.delim("Picadas.csv", colClasses = c("factor", "factor", "numeric", "factor", "logical"), comment.char = "#")
## algunas categorías combinadas
#Pogos_picadas$Cats <- as.factor(paste0(Pogos_picadas$Picada, Pogos_picadas$Fecha, Pogos_picadas$Sp))
#Pogos_picadas$Sample <- as.factor(paste0(Pogos_picadas$Picada, Pogos_picadas$Fecha))

# importación de planilla de topografía
Alt_picadas <- read.delim("Altura_picadas.csv", colClasses = c("factor", "numeric", "integer", "numeric", "numeric", "numeric"), comment.char = "#")

# importación de planilla de permeabilidad
Perm_picadas <- read.delim("Permeabilidad.csv", colClasses = c("factor", "numeric", "factor", "integer"), comment.char = "#", na.strings = "")

# importacion de datos de textura de suelo
Textura <- read.delim("Textura_suelo.csv", colClasses = c("character", "factor", "factor", rep("numeric", 3), "factor"), comment.char = "#", na.strings = "")
levels(Textura$Especie) <- c("Random", "PI", "PM")

Textura_long <- Textura %>%
  pivot_longer(cols = c("Arena", "Limo", "Arcilla"), names_to = "Tamaño", values_to = "Percentage")

TexturaXCluster <- Textura %>%
  group_by(Cluster) %>%
  summarize(Especie = unique(Especie), Arena = mean(Arena), Limo = mean(Limo), Arcilla = mean(Arcilla))

#variable para definir SPI (índice de permeabilidad) invirtiendo el log(Tiempo) restándole log(120 sec) para más o menos centrarlo en valores típicos del algarrobal
Perm_picadas$SPI <- -(log(Perm_picadas$Tiempo) - log(120))

#Generic window to contain any 1500m transect; width x10 (60m) to facilitate graphics
Ventana <- owin(c(0,1500), c(0,40))

#Generic network with single 1500m segment in the middle of the window
Picada <- linnet(ppp(x = c(0, 1500), y = c(30, 30), window = Ventana), m = matrix(c(FALSE, TRUE, TRUE, FALSE), 2, 2))

#Transectas para análisis con spatstat (solo colonias adultas) incluyendo marcas de Fecha y de Especie [por si se pudiera analizar espacio/tiempo junto como en Baddeley et al 2020]. Por ahora, cada transecta*fecha va por separado [abajo]
#LaFad_Pogos <- Pogos_picadas %>%
#                filter(Picada == "LaFad", Adlt_Actv == T) %>%
#                transmute(x = Dist, y = 3, Sp, Year = Fecha) %>%
#                droplevels() %>%
#                lpp(L = Picada)
#Doble_Pogos <- Pogos_picadas %>%
#                filter(Picada == "Doble", Adlt_Actv == T) %>%
#                transmute(x = Dist, y = 3, Sp, Year = Fecha) %>%
#                droplevels() %>%
#                lpp(L = Picada)

#Transectas para análisis con spatstat (solo colonias adultas) con marcas de Especie [una por transecta*año] 
LaFad2004_Pogos <- Pogos_picadas %>%
                    filter(Picada == "LaFad", Fecha == "2004", Adlt_Actv == T) %>%
                    transmute(x = Dist, y = 30, Sp) %>%
                    droplevels() %>%
                    lpp(L = Picada)

LaFad2007_Pogos <- Pogos_picadas %>%
                    filter(Picada == "LaFad", Fecha == "2007", Adlt_Actv == T) %>%
                    transmute(x = Dist, y = 30, Sp) %>%
                    droplevels() %>%
                    lpp(L = Picada)

LaFad2010_Pogos <- Pogos_picadas %>%
                    filter(Picada == "LaFad", Fecha == "2010", Adlt_Actv == T) %>%
                    transmute(x = Dist, y = 30, Sp) %>%
                    droplevels() %>%
                    lpp(L = Picada)

Doble2004_Pogos <- Pogos_picadas %>%
                    filter(Picada == "Doble", Fecha == "2004", Adlt_Actv == T) %>%
                    transmute(x = Dist, y = 30, Sp) %>%
                    droplevels() %>%
                    lpp(L = Picada)

Doble2007_Pogos <- Pogos_picadas %>%
                    filter(Picada == "Doble", Fecha == "2007", Adlt_Actv == T) %>%
                    transmute(x = Dist, y = 30, Sp) %>%
                    droplevels() %>%
                    lpp(L = Picada)

Doble2010_Pogos <- Pogos_picadas %>%
                    filter(Picada == "Doble", Fecha == "2010", Adlt_Actv == T) %>%
                    transmute(x = Dist, y = 30, Sp) %>%
                    droplevels() %>%
                    lpp(L = Picada)

Torre2006_Pogos <- Pogos_picadas %>%
                    filter(Picada == "Torre", Adlt_Actv == T) %>%
                    transmute(x = Dist, y = 30, Sp) %>%
                    lpp(L = Picada)

#Idem alturas
LaFad_Alt <- Alt_picadas %>%
              filter(Picada == "LaFad") %>%
              transmute(x, y = 30, Altitude) %>%
              lpp(L = Picada)

Doble_Alt <- Alt_picadas %>%
              filter(Picada == "Doble") %>%
              transmute(x, y = 30, Altitude) %>%
              lpp(L = Picada)

#Idem permeabilidad
LaFad_Perm <- Perm_picadas %>%
                filter(Lugar == "LaFad") %>%
                transmute(x = Dist, y = 30, SPI) %>%
                lpp(L = Picada)

Doble_Perm <- Perm_picadas %>%
                filter(Lugar == "Doble") %>%
                transmute(x = Dist, y = 30, SPI) %>%
                lpp(L = Picada)

Picadas2004 <- anylist(LaFad2004 = LaFad2004_Pogos, Doble2004 = Doble2004_Pogos)
PicadasAll <- list(LaFad = anylist(LaFad2004_Pogos, LaFad2007_Pogos, LaFad2010_Pogos), Doble = anylist(Doble2004_Pogos, Doble2007_Pogos, Doble2010_Pogos))

# function to shift marked points patterns (categorical marks) on single-lined network (lpp, seg == 1, variable x, fixed y) with toroidal shift to test independence of components [similar to rshift.ppp but for a line; toroidal shift based on upper x limit or endline; x values rounded as data]
linear_shift <- function(LineToShift) {
  for (i in levels(LineToShift$data$marks)) {
    LineToShift$data$tp[LineToShift$data$marks == i] <- LineToShift$data$tp[LineToShift$data$marks == i] + round(runif(1), 4)  
  }
  LineToShift$data$tp[LineToShift$data$tp > 1] <- LineToShift$data$tp[LineToShift$data$tp > 1] - 1
  LineToShift$data$x <- LineToShift$data$tp * LineToShift$domain$lines$ends$x1
  return(LineToShift)
}

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

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

Results

Abundance and spatio-temporal dynamics of Pogonomyrmex colonies on dirt roads

Tens of active Pogonomyrmex nests were detected along each 1.5 km dirt road segment sampled, with all three species present in every single case (Figure 1).

ggplot(data = Pogos_picadas, aes(x = Dist, y = Sp, colour = Sp, shape = Adlt_Actv)) + geom_point(fill = "white", size = .9) + labs(x = "Distance (m)") + tema + theme(axis.title.y = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank()) + facet_grid(Picada + Fecha ~ .) + scale_shape_manual(values = c(42, 19)) + scale_colour_manual(values = c("darkgreen", "darkblue", "darkred")) + scale_y_discrete(limits=c("", "","PI","PM", "PR", "", ""))  

Figure 1: Active Pogonomyrmex colonies of each species (PI: green, PM: blue, PR: red) detected along 1.5 km segments of three dirt roads (LaFad, Torre and Doble) during their peak external activity season (February) of different years (2004–2010). Species are vertically separated to ease visualization, but only the distance along the paths (and not across) was recorded. Adult colonies are represented with big dots, and incipient colonies (only single queens) with small points.

The number of active adult colonies of Pogonomyrmex inermis (PI) and Pogonomyrmex rastratus (PR) was always higher than Pogonomyrmex mendozanus (PM) in the 1.5 km segments of both dirt roads that were censused in 2004 and 2007 (Table 1). The dirt road sampled (only) in 2006 showed the same pattern (PI: 54 nests, PR: 33; PM: 22).

All species reduced their average abundances by 2010, though more strikingly PR (Figure 2): in 2004 it accounted for 46.5% of Pogonomyrmex active nests in La Fadina and Doble, 35.9% in 2007, and only 26.0% in 2010.

kable(with(subset(Pogos_picadas, Adlt_Actv == TRUE), table(Fecha, Sp, exclude = "2006")), format = "html", table.attr = "style='width:70%; margin-right:auto; margin-left:auto;'", caption = "Number of active colonies of each *Pogonomyrmex* species detected in different summers (2004, 2007 and 2010) on the same two 1.5 km segments of dirt roads crossing the algarrobal in Ñacuñán, central Monte desert.")

Table 1: Number of active colonies of each Pogonomyrmex species detected in different summers (2004, 2007 and 2010) on the same two 1.5 km segments of dirt roads crossing the algarrobal in Ñacuñán, central Monte desert.

	PI	PM	PR
2004	48	28	66
2007	59	25	47
2010	38	16	19

#Densidad promedio de nidos activos (adultos) en caminos dentro del algarrobal
ggplot(subset(Pogos_picadas, Adlt_Actv == TRUE & Fecha != "2006"), aes(x = Fecha, colour = Sp, linetype = Picada)) +
  geom_line(aes(group = interaction(Sp, Picada)), stat= "count") +
  geom_point(aes(group = interaction(Sp, Picada)), stat= "count") +
  labs(x = "Year", y = "Active nests (1.5 km)", colour = "Species", linetype = "Dirt road") +
  scale_colour_manual(values = c("darkgreen", "darkblue", "darkred")) +
  tema + theme(legend.position = "right")

Figure 2: Number of adult Pogonomyrmex colonies of each species in the two 1.5 km segments of two dirt roads (LaFad and Doble) sampled repeatedly in 2004, 2007 and 2010.

Some incipient colonies (i.e., with only a single founding queen detected, no workers yet) were also found occasionally, most of them of P. mendozanus concentrated in a ~200 m sector of Doble during February 2010 (Table 2). Their rate of establishment/survival in the area proved rather low, so they were not considered for further analyses.

# colonias no adultas:
kable(with(subset(Pogos_picadas, Adlt_Actv == FALSE), table(Fecha, Sp)), format = "html", table.attr = "style='width:70%; margin-right:auto; margin-left:auto;'", caption = "Number of incipient colonies (only single queens excavating or foraging) of each *Pogonomyrmex* species detected in different summers (2004, 2007 and 2010) on the same two 1.5 km segments of dirt roads crossing the algarrobal in Ñacuñán, central Monte desert.")

Table 2: Number of incipient colonies (only single queens excavating or foraging) of each Pogonomyrmex species detected in different summers (2004, 2007 and 2010) on the same two 1.5 km segments of dirt roads crossing the algarrobal in Ñacuñán, central Monte desert.

	PI	PM	PR
2004	0	0	1
2006	0	17	0
2007	0	6	0
2010	0	104	0

Intra- and interspecific spatial patterns (LaFad & Doble 2004)

The intensity of the spatial process underlying the distribution of active Pogonomyrmex spp. nests (i.e., all three species together) does not change notoriously along the road segments compared to what is expected under a CSR (random) process, both when describing its nonparametric estimate at each distance and with spatial summary functions (K and pair correlation) that assume spatially homogeneous intensities (Figure 3). However, heterogeneous intensities are notorious when analysing linear point patterns by species (comparing them against expectations for CSRI processes, i.e., CSR + Independence among species). Both goodness-of-fit tests between Poisson models of observed patterns and CSR processes and global tests of the K and pair-correlation functions (that test for autocorrelation assuming homogeneity) suggest heterogeneous processes with intra-specific aggregation in several cases (Figure 3).

par(mar = c(2, 2, 1, .5))
#layout to incluge subfigures by genus and by species together for each transect
layout(matrix(c(1:3, 7:15, 4:6, 16:24), 8, 3, byrow = TRUE), widths = c(3, 1, 1))

#Intensities for Pogo colonies (all species together)

#loop of transects with vector of values for vertical text position of tests results within plots
loop.texts <- data.frame(Main = c("Intensity (LaFad 2004)", "Intensity (Doble 2004)"), Lambda = c(0.08, 0.09), K = c(-15, -15), PCF = c(0.8, 0.8))
loop.count <- 0

for(i in Picadas2004) {
  loop.count <- loop.count + 1
  plot(rhohat(unmark(i), "x"), legend = FALSE, main = loop.texts$Main[loop.count], cex.axis = .8)
  text(1525, y = loop.texts$Lambda[loop.count], paste("p = ", round(cdf.test(unmark(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). fix.n = TRUE ver Baddeley et al 2014 Ecol.Monog.
  temp <- envelope(unmark(i), linearK, nsim = 199, nrank = 5, r = seq(0, 350, 10), verbose = FALSE, fix.n = TRUE, savefuns = TRUE, warnWbw = FALSE)
  plot(temp, cbind(obs - theo, lo = lo - theo, hi = hi - theo) ~ r, main = "K-r ~ r", legend = FALSE, cex.axis = 0.8)
  text(350, 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 = "red", lwd = 0.5)
    
  temp <- envelope(unmark(i), linearpcf, nsim = 199, nrank = 5, r = seq(0, 350, 10), verbose = FALSE, fix.n = TRUE, savefuns = TRUE, warnWbw = FALSE)
  plot(temp, main = "pcf ~ r", legend = FALSE, cex.axis = 0.8)
  text(350, y = loop.texts$PCF[loop.count], paste("p = ", round(mad.test(temp, use.theory = TRUE)$p.value, 3)), adj = c(1, 1))
}

#Intensities by species

#double loop transect/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("LaFad: PI", "LaFad: PM", "LaFad: PR", "Doble: PI", "Doble: PM", "Doble: PR"), Lambda = c(0.03, 0.025, 0.065, 0.045, 0.027, 0.047), K = c(-40, -60, -20, -35, -50, -30), PCF = c(2, 2.5, 1.45, 2.5, 1.8, 1.4), Col = c("darkgreen", "darkblue", "darkred", "darkgreen", "darkblue", "darkred"))
loop.count <- 0

for(i in Picadas2004) {
  for(j in split(i)) {
    loop.count <- loop.count + 1
    plot(rhohat(j, "x"), legend = FALSE, main = loop.texts$Main[loop.count], cex.axis = .8, col = loop.texts$Col[loop.count], shadecol = alpha(loop.texts$Col[loop.count], 0.3))
    text(1525, y = loop.texts$Lambda[loop.count], paste("p = ", round(cdf.test(j, "x", test = "ks", jitter = F)$p.value, 3)), adj = c(1, 1))
    
    temp <- envelope(j, linearK, nsim = 199, nrank = 5, r = seq(0, 350, 10), verbose = FALSE, fix.n = TRUE, savefuns = TRUE, savepatterns = TRUE, warnWbw = FALSE)
    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(350, 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(j, linearpcf, nsim = 199, nrank = 5, r = seq(0, 350, 10), verbose = FALSE, simulate = temp,  savefuns = TRUE, warnWbw = FALSE)
    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(350, 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)
  }
}

## Ties were encountered

Figure 3: Left: Nonparametric estimate of the intensity of a Poisson spatial point process for all Pogonomyrmex spp. colonies and for each species separately along each of the two segments of dirt road sampled in 2004. P values result from Kolmogorov-Smirnov goodness-of-fit tests of single-predictor (the x coordinate) Poisson models fitted to the observed linear point process against a CSR model. Center and right: Summary functions K (center) and pair correlation (right) for homogeneous Poisson linear point processes as functions of the distance (r) among colonies. 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 Maximum Absolute Deviation (or ‘global’ tests) against CSR.

Those intraspecific aggregations may result either from Poisson processes that vary their intensity along the dirt roads (with each species responding independently) or from rather homogeneous base conditions where aggregation/segregation results from association/dissociation interactions among species (=positive/negative interactions among their colonies).

Observed patters are compatible with Poisson spatial processes with heterogeneous intensities (Figure 4), i.e., the probability of finding a nest varies along each road segment, with no important evidence to support positive or negative interactions among neighbouring colonies (i.e., of not being a Poisson process, which assumes distance-related independence among colonies).

par(mar = c(2, 2, 1, .5))
layout(matrix(1:9, 3, 3, byrow = TRUE), widths = c(3, 1, 1))
loop.texts <- data.frame(Main = c("LaFad: PI", "LaFad: PM", "LaFad: PR", "Doble: PI", "Doble: PM", "Doble: PR"), K = c(-200, -100, 100, -130, -130, -75), PCF = c(1.8, 2.5, 1.5, 2.0, 2.0, 1.6), Col = c("darkgreen", "darkblue", "darkred", "darkgreen", "darkblue", "darkred"))
loop.count <- 0

#based in lppm Poisson model with density(sigma = 100) as variable lambda estimate
#loop for transects
for(i in Picadas2004) {
  #loop for species
  for(j in split(i)) {
    loop.count <- loop.count + 1
    
    #estimate density with gaussian kernel and fixed sigma in case the model asks for it
    est_density <- densityfun(j, sigma = 100, dx = 5)
    #define model for sp j in transect i
    templppm <- lppm(j ~ log(est_density))
    #plot the model with K-S testing fit of observed against model
    plot(templppm, legend = FALSE, main = loop.texts$Main[loop.count], cex.axis = .8, style = "w", adjust = 30, col = loop.texts$Col[loop.count])
    text(1525, y = 200, 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, sigma = 100, r = seq(0, 350, 10), normpower = 2, nsim = 199, nrank = 5, verbose = FALSE, fix.n = TRUE, savefuns = TRUE, savepatterns = TRUE, warnWbw = FALSE)
    #plot pointwise envelope with (single-staged) MAD test according to linearKinhom [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 = loop.texts$Col[loop.count], shadecol = alpha(loop.texts$Col[loop.count], 0.3))
    text(350, y = loop.texts$K[loop.count], paste("p = ", round(mad.test(temp, use.theory = FALSE)$p.value, 3)), adj = c(1, 0), col = loop.texts$Col[loop.count])
    #(two-staged) MAD global test based on linearKinhom [Baddeley et al 2017]
    text(350, y = loop.texts$K[loop.count], paste("p = ", round(bits.test(templppm, fun = linearKinhom, exponent = Inf, rinterval = c(0, 350), nsim = 19, leaveoneout = FALSE, use.theory = FALSE, verbose = FALSE)$p.value, 3)), adj = c(1, 1.5), col = loop.texts$Col[loop.count], font = 2)
    abline(h = 0, lty = 2, col = loop.texts$Col[loop.count], lwd = 0.5)
    # idem with linearpcfinhom
    temp2 <- envelope(templppm, linearpcfinhom, sigma = 100, r = seq(0, 350, 10), normpower = 2, simulate = temp, nsim = 199, nrank = 5, verbose = FALSE, fix.n = TRUE, savefuns = TRUE, warnWbw = FALSE)
    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$Col[loop.count], shadecol = alpha(loop.texts$Col[loop.count], 0.3))
    text(350, y = loop.texts$PCF[loop.count], paste("p = ", round(mad.test(temp2, use.theory = FALSE)$p.value, 3)), adj = c(1, 0), col = loop.texts$Col[loop.count])
    #(two-staged) MAD global test based on linearpcfinhom [Baddeley et al 2017]
    text(350, y = loop.texts$PCF[loop.count], paste("p = ", round(bits.test(templppm, fun = linearpcfinhom, exponent = Inf, rinterval = c(0, 350), nsim = 19, leaveoneout = FALSE, use.theory = FALSE, verbose = FALSE)$p.value, 3)), adj = c(1, 1.5), col = loop.texts$Col[loop.count], font = 2)
    abline(h = 1, lty = 2, col = loop.texts$Col[loop.count], lwd = 0.5)
  }
}

Figure 4: Left: Intensities of the heterogeneous Poisson linear point process model estimated through gaussian kernels (with fixed smoothing bandwidth: sigma = 100) for the active colonies of each Pogonomyrmex spp. along each of the two segments of dirt roads sampled in 2004. P values result from Kolmogorov-Smirnov goodness-of-fit tests between the observed pattern and the Poisson model fitted. Center and right:Summary functions K (center) and pair correlation (right) for inhomogeneous point patterns as functions of the distance (r) among colonies. Envelopes show the 95% critical points for Monte Carlo pointwise tests estimated after 199 iterations of the fitted heterogeneous Poisson model. P values in normal font result from single-staged (‘global’) tests of Maximum Absolute Deviation (MAD) against the fitted model, which may be too conservative (Baddeley et al. 2017); p values in bold result from a Balanced Independent Two-Stage MonteCarlo global test of goodness-of-fit (with 19 patterns per stage for a total of 380 simulations) that avoids that bias.

Those intraspecific aggregated point patterns are not independent among species, resulting in interspecific spatial segregation. Tests of independence through random labelling, which are conditional on the observed linear distribution of Pogonomyrmex spp. nests, clearly show this non-random distribution (Figure 5): in Doble, PM appears segregated from the other two species (which overlap their distributions); in La Fadina PM segregates from PR. Tests of independence by random toroidal shift, conditional on the three species-specific observed distributions, are much more conservative but also show some evidence compatible with interspecific spatial segregation, particularly in Doble.

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

loop.texts <- data.frame(Main = c("LaFad: PM / PI", "LaFad: PM / PR", "LaFad: PI / PR", "Doble: PM / PI", "Doble: PM / PR", "Doble: PI / PR"), K = c(-30, -25, 40, 20, 15, -25), PCF = c(1.3, 0.85, 0.85, 0.75, 0.8, 1.3))
loop.count <- 0

for(i in Picadas2004) {
  for(j 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)
    #function averages linearKcross ij and ji
    temp <- envelope(i, linearKcross.inhom.avrg, j[1], j[2], sigma = 100, r = seq(0, 350, 10), nsim = 199, nrank = 5, simulate = expression(rlabel(i)), verbose = FALSE, savefuns = TRUE, savepatterns = TRUE)
    # test of independence-of-components con función linear_shift para generar patrón con puntos desplazados (= toroidal-shift)
    temp2 <- envelope(i, linearKcross.inhom.avrg, j[1], j[2], sigma = 100, r = seq(0, 350, 10), nsim = 199, nrank = 5, simulate = expression(linear_shift(i)), verbose = FALSE, savefuns = TRUE, savepatterns = TRUE)
    
    plot(temp2, cbind(obs - mmean, lo = lo - mmean, hi = hi - mmean) ~ r, main = loop.texts$Main[loop.count], legend = FALSE, cex.axis = .8, shadecol = "#00000030")
    plot(temp, cbind(obs - mmean, lo = lo - mmean, hi = hi - mmean) ~ r, main = loop.texts$Main[loop.count], legend = FALSE, cex.axis = .8, col = "#996600", shadecol = "#99660030", add = TRUE)
    text(350, y = loop.texts$K[loop.count], paste("p = ", round(mad.test(temp, use.theory = FALSE)$p.value, 3)), col = "#996600", adj = c(1, 1.5))
    text(350, y = loop.texts$K[loop.count], paste("p = ", round(mad.test(temp2, use.theory = FALSE)$p.value, 3)), adj = c(1, 0))
    abline(h = 0, lty = 2, col = "red", lwd = 0.5)
    
    temp3 <- envelope(i, linearpcfcross.inhom.avrg, j[1], j[2], sigma = 100, r = seq(0, 350, 10), simulate = temp, nsim = 199, nrank = 5, verbose = FALSE, savefuns = TRUE)
    temp4 <- envelope(i, linearpcfcross.inhom.avrg, j[1], j[2], sigma = 100, r = seq(0, 350, 10), simulate = temp2, nsim = 199, nrank = 5, verbose = FALSE, savefuns = TRUE)
    
    plot(temp4, cbind(obs - mmean + 1, lo = lo - mmean + 1, hi = hi - mmean + 1) ~ r, main = "", legend = FALSE, cex.axis = 0.8, shadecol = "#00000030")
    plot(temp3, cbind(obs - mmean + 1, lo = lo - mmean + 1, hi = hi - mmean + 1) ~ r, main = "", legend = FALSE, cex.axis = 0.8, col = "#996600", shadecol = "#99660030", add = TRUE)
    text(350, y = loop.texts$PCF[loop.count], paste("p = ", round(mad.test(temp3, use.theory = FALSE)$p.value, 3)), col = "#996600", adj = c(1, 1.5))
    text(350, y = loop.texts$PCF[loop.count], paste("p = ", round(mad.test(temp4, use.theory = FALSE)$p.value, 3)), adj = c(1, 0))
    abline(h = 1, lty = 2, col = "red", lwd = 0.5)
  }
}

Figure 5: Tests of independence among the distributions of active nests of each Pogonomyrmex species along each of the two segments of dirt roads sampled in 2004, based on the inhomogeneous cross-type linear versions of two summary functions, K (cumulative, above) and pair correlation (instantaneous, below), for point patterns as functions of the distance (r) among colonies. Values lower than 0 or 1, respectively, suggest inhibition between points; values greater than 0 or 1 suggest clustering. Envelopes show the 95% critical points for Monte Carlo pointwise tests estimated after 199 iterations of two different null models: random labelling (conditional on the global inhomogeneity; brown) and toroidal shift (conditional on both inhomogeneous point patterns; grey), with p values from (‘global’) tests of Maximum Absolute Deviation (MAD).

In fact, the probability that an active nest belongs to each species changes along each dirt road, with one particularly notorious case: sectors where PM is more probable are complementary of those where the probability of PI increases (with more intense segregation in Doble; Figure 6).

#Fig intensidades x spp; sigma "óptimo promedio" (bw.relrisk = 40–60)
par(mar = c(2, 2, 1, .5))
layout(matrix(1:2))

loop.spp <- data.frame(Sp = c("PI", "PM", "PR"), Cols = c("darkgreen", "darkblue", "darkred"))

for(i in 1:length(Picadas2004)) { #picadas (panel)
  temp <- relrisk(Picadas2004[[i]], sigma = 50, dx = 5)
  for(k in 1:length(temp)) { #sp (curvas)
    plot(temp[[k]], style = "width", scale = 100, legend = F, main = names(Picadas2004)[i], col = alpha(loop.spp$Cols[k], 0.3), add = k > 1)
    plot(split(Picadas2004[[i]])[[k]], add = T, col = loop.spp$Cols[k], pch = 19, size = 0.5, lty = 0)
  }
}

Figure 6: Spatially-varying probability of a nest being of each of the three Pogonomyrmex species (PI: green; PM: blue; PR: red) along the two dirt road sampled in 2004, estimated with a fixed bandwidth (sigma = 50 m) and represented by the width of the coloured band, with the observed active nests as dots.

Assuming Poisson linear point processes (= no interactions between neighbouring colonies), models both parametric (log-cubic relationship with the coordinate) and non-parametric (smoothed cubic splines) posing species-specific variations along the dirt road segments (i.e., with one parameter per species; Models 2) fit the observed patterns much better than those with common parameters (Models 1):

#con modelo paramétricos (polinomio)
#rescale es para pasar escala a km y evitar error de matriz singular; no cambia nada
cat("Parametric models (cubic polynomials)\n")
cat("La Fadina\n")
anova(lppm(rescale(LaFad2004_Pogos, 1000) ~ polynom(x, 3) + marks), lppm(rescale(LaFad2004_Pogos, 1000) ~ polynom(x, 3) * marks), test = "LRT")
cat("\nDoble\n")
anova(lppm(rescale(Doble2004_Pogos, 1000) ~ polynom(x, 3) + marks), lppm(rescale(Doble2004_Pogos, 1000) ~ polynom(x, 3) * marks), test = "LRT")

#con modelo no paramétrico (cubic spline)
cat("\nNon-parametric models (cubic splines, df = 5)\n")
cat("La Fadina\n")
anova(lppm(LaFad2004_Pogos ~ splines::bs(x, df = 5) + marks), lppm(LaFad2004_Pogos ~ splines::bs(x, df = 5) * marks), test = "LRT")
cat("\nDoble\n")
anova(lppm(Doble2004_Pogos ~ splines::bs(x, df = 5) + marks), lppm(Doble2004_Pogos ~ splines::bs(x, df = 5) * marks), test = "LRT")

## Parametric models (cubic polynomials)
## La Fadina
## Analysis of Deviance Table
## 
## Model 1: ~x + I(x^2) + I(x^3) + marks     Poisson
## Model 2: ~(x + I(x^2) + I(x^3)) * marks   Poisson
##   Npar Df Deviance Pr(>Chi)   
## 1    6                        
## 2   12  6   22.118 0.001153 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Doble
## Analysis of Deviance Table
## 
## Model 1: ~x + I(x^2) + I(x^3) + marks     Poisson
## Model 2: ~(x + I(x^2) + I(x^3)) * marks   Poisson
##   Npar Df Deviance    Pr(>Chi)    
## 1    6                            
## 2   12  6   36.564 0.000002141 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Non-parametric models (cubic splines, df = 5)
## La Fadina
## Analysis of Deviance Table
## 
## Model 1: ~splines::bs(x, df = 5) + marks      Poisson
## Model 2: ~splines::bs(x, df = 5) * marks      Poisson
##   Npar Df Deviance       Pr(>Chi)    
## 1    8                               
## 2   18 10   58.684 0.000000006426 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Doble
## Analysis of Deviance Table
## 
## Model 1: ~splines::bs(x, df = 5) + marks      Poisson
## Model 2: ~splines::bs(x, df = 5) * marks      Poisson
##   Npar Df Deviance        Pr(>Chi)    
## 1    8                                
## 2   18 10    64.68 0.0000000004667 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

This is strong evidence for species-specific density variations along the dirt roads, suggesting they respond differently to the fluctuations of environmental conditions (or, less probable, that they keep exclusive sectors through negative interactions among species but they are not related to the distance among their individual colonies).

In summary, even when Pogonomyrmex spp. nests can be found with almost homogeneous density along the dirt roads, their species identities do not occur randomly. Each species show intra- and interspecific patterns of variation in their densities along the sampled segments, probably caused by species-specific responses to underlying variations in environmental conditions, resulting in 3–7 patches/species in each 1.5 km censused.

¿Does heterogeneous, species-specific spatial patterns repeat elsewhere?: Torre 2006

In the third dirt road segment, sampled two years later and where the nest density was higher, all intra- and inter-specific patters agree with what was described for the two focal dirt roads in 2004: even when the segment was mostly occupied by Pogonomyrmex spp. nests all along (though a bit more variably than in the other two), colonies of two of the three species (PM and PI) showed very aggregated linear distributions, with the other species (PR) occurring more homogeneously (Figure 7).

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

#loop Pogos+species with vector of values for vertical text position of tests results within plots
loop.texts <- data.frame(Main = c("Torre: Pogonomyrmex", "Torre: PI", "Torre: PM", "Torre: PR"), Lambda = c(0.12, 0.13, 0.05, 0.045), K = c(-15, -15, -40, -30), PCF = c(1.15, 2.3, 1.8, 1.35), Col = c("black", "darkgreen", "darkblue", "darkred"))
loop.count <- 0

for(j in c(Pogos = list(unmark(Torre2006_Pogos)), split(Torre2006_Pogos))) {
  loop.count <- loop.count + 1
  plot(rhohat(j, "x"), legend = FALSE, main = loop.texts$Main[loop.count], cex.axis = .8, col = loop.texts$Col[loop.count], shadecol = alpha(loop.texts$Col[loop.count], 0.3))
  text(1525, y = loop.texts$Lambda[loop.count], paste("p = ", round(cdf.test(j, "x", test = "ks", jitter = F)$p.value, 3)), adj = c(1, 1))
    
  temp <- envelope(j, linearK, nsim = 199, nrank = 5, r = seq(0, 350, 10), verbose = FALSE, fix.n = TRUE, savefuns = TRUE, savepatterns = TRUE, warnWbw = FALSE)
  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(350, 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(j, linearpcf, nsim = 199, nrank = 5, r = seq(0, 350, 10), verbose = FALSE, simulate = temp, fix.n = TRUE, savefuns = TRUE, warnWbw = FALSE)
  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(350, y = loop.texts$PCF[loop.count], paste("p = ", round(mad.test(temp, 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 7: Left: Nonparametric estimate of the intensity of a Poisson spatial point process for all Pogonomyrmex spp. colonies (above) and for each species (PI: green; PM: blue; PR: red) along the 1.5 km segment of dirt road sampled in 2006. P values result from Kolmogorov-Smirnov goodness-of-fit tests of single-predictor (the x coordinate) Poisson models fitted to the observed linear point process against a CSR model. Center and right: Summary functions K (center) and pair correlation (right) for homogeneous Poisson linear point processes as functions of the distance (r) among colonies. 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 Maximum Absolute Deviation (or ‘global’ tests) against CSR.

These aggregated patterns are again compatible with Poisson point process models with inhomogeneous intensities along the segment (i.e., with no evidence of distance-related atraction/repulsion interactions between individual colonies; Fig. 8).

par(mar = c(2, 2, 1, .5))
layout(matrix(1:12, 4, 3, byrow = TRUE), widths = c(3, 1, 1))
loop.texts <- data.frame(Main = c("Torre: Pogonomyrmex", "Torre: PI", "Torre: PM", "Torre: PR"), K = c(5, -200, -150, -150), PCF = c(1.15, 1.9, 1.8, 1.75), Col = c("black", "darkgreen", "darkblue", "darkred"))
loop.count <- 0

#based in density() como estimador de lambda; simula las iteraciones siguiendo esa heterogeneidad (con simulate = expression(rpoislpp(lambda = density(j, sigma = 100))))
for(j in c(Pogos = list(unmark(Torre2006_Pogos)), split(Torre2006_Pogos))) {
    loop.count <- loop.count + 1
    
    #estimate density with gaussian kernel and fixed sigma in case the model asks for it
    est_density <- densityfun(j, sigma = 100, dx = 5)
    #define model for sp j in transect i
    templppm <- lppm(j ~ log(est_density))
    #plot the model with K-S testing fit of observed against model
    plot(templppm, legend = FALSE, main = loop.texts$Main[loop.count], cex.axis = .8, style = "w", adjust = 30, col = loop.texts$Col[loop.count])
    text(1525, y = 150, 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, sigma = 100, r = seq(0, 350, 10), normpower = 2, nsim = 199, nrank = 5, verbose = FALSE, fix.n = TRUE, savefuns = TRUE, savepatterns = TRUE, warnWbw = FALSE)
    #plot pointwise envelope with (single-staged) MAD test according to linearKinhom [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 = loop.texts$Col[loop.count], shadecol = alpha(loop.texts$Col[loop.count], 0.3))
    text(350, y = loop.texts$K[loop.count], paste("p = ", round(mad.test(temp, use.theory = FALSE)$p.value, 3)), adj = c(1, 0), col = loop.texts$Col[loop.count])
    #(two-staged) MAD global test based on linearKinhom [Baddeley et al 2017]
    text(350, y = loop.texts$K[loop.count], paste("p = ", round(bits.test(templppm, fun = linearKinhom, exponent = Inf, rinterval = c(0, 350), nsim = 19, leaveoneout = FALSE, use.theory = FALSE, verbose = FALSE, warnWbw = FALSE)$p.value, 3)), adj = c(1, 1.5), col = loop.texts$Col[loop.count], font = 2)
    abline(h = 0, lty = 2, col = loop.texts$Col[loop.count], lwd = 0.5)
    # idem with linearpcfinhom
    temp2 <- envelope(templppm, linearpcfinhom, sigma = 100, r = seq(0, 350, 10), normpower = 2, simulate = temp, nsim = 199, nrank = 5, verbose = FALSE, fix.n = TRUE, savefuns = TRUE, warnWbw = FALSE)
    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$Col[loop.count], shadecol = alpha(loop.texts$Col[loop.count], 0.3))
    text(350, y = loop.texts$PCF[loop.count], paste("p = ", round(mad.test(temp2, use.theory = FALSE)$p.value, 3)), adj = c(1, 0), col = loop.texts$Col[loop.count])
    #(two-staged) MAD global test based on linearpcfinhom [Baddeley et al 2017]
    text(350, y = loop.texts$PCF[loop.count], paste("p = ", round(bits.test(templppm, fun = linearpcfinhom, exponent = Inf, rinterval = c(0, 350), nsim = 19, leaveoneout = FALSE, use.theory = FALSE, verbose = FALSE, warnWbw = FALSE)$p.value, 3)), adj = c(1, 1.5), col = loop.texts$Col[loop.count], font = 2)
    abline(h = 1, lty = 2, col = loop.texts$Col[loop.count], lwd = 0.5)
}

Figure 8: Same as Fig. 4 but for all Pogonomyrmex nests (black), and then split by species, in the 1.5 km segment of Torre dirt road censused in 2006.

The intensity of the spatial processes of each species vary along the road, but also the probability of an active nest belonging to each species changes as before: sectors where PM predominates appear complementary to those where the chance of PI increases (Fig. 9).

#Fig intensidades x spp; sigma "óptimo promedio" (bw.relrisk = 40–60)
par(mar = c(2, 2, 1, .5))
layout(matrix(1:2))

loop.spp <- data.frame(Sp = c("PI", "PM", "PR"), Cols = c("darkgreen", "darkblue", "darkred"))

temp <- relrisk(Torre2006_Pogos, sigma = 50, dx = 5)
for(k in 1:length(temp)) { #sp (curvas)
  plot(temp[[k]], style = "width", scale = 80, legend = F, main = "Torre2006", col = alpha(loop.spp$Cols[k], 0.3), add = k > 1)
  plot(split(Torre2006_Pogos)[[k]], add = T, col = loop.spp$Cols[k], pch = 19, size = 0.5, lty = 0)
}

Figure 9: Same as Fig. 6 but along the segment of dirt road sampled in 2006.

Accordingly, the two species with strongly heterogeneous patterns (PM and PI) are again spatially segregated (Fig. 10).

# linearKcross: tests of random labelling and independence-of-components (random toroidal shift)
par(mar = c(2, 2, 1, .5))
layout(matrix(1:6, 2, 3, byrow = FALSE))

loop.texts <- data.frame(Main = c("Torre: PM / PI", "Torre: PM / PR", "Torre: PI / PR"), K = c(20, -20, -20), PCF = c(1.1, 0.9, 0.9))
loop.count <- 0

for(i in list(Torre2006_Pogos)) {
  for(j 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)
    #function averages linearKcross ij and ji
    temp <- envelope(i, linearKcross.inhom.avrg, j[1], j[2], sigma = 100, r = seq(0, 350, 10), nsim = 199, nrank = 5, simulate = expression(rlabel(i)), verbose = FALSE, savefuns = TRUE, savepatterns = TRUE)
    # test of independence-of-components con función linear_shift para generar patrón con puntos desplazados (= toroidal-shift)
    temp2 <- envelope(i, linearKcross.inhom.avrg, j[1], j[2], sigma = 100, r = seq(0, 350, 10), nsim = 199, nrank = 5, simulate = expression(linear_shift(i)), verbose = FALSE, savefuns = TRUE, savepatterns = TRUE)
    
    plot(temp2, cbind(obs - mmean, lo = lo - mmean, hi = hi - mmean) ~ r, main = loop.texts$Main[loop.count], legend = FALSE, cex.axis = .8, shadecol = "#00000030")
    plot(temp, cbind(obs - mmean, lo = lo - mmean, hi = hi - mmean) ~ r, main = loop.texts$Main[loop.count], legend = FALSE, cex.axis = .8, col = "#996600", shadecol = "#99660030", add = TRUE)
    text(350, y = loop.texts$K[loop.count], paste("p = ", round(mad.test(temp, use.theory = FALSE)$p.value, 3)), col = "#996600", adj = c(1, 1.5))
    text(350, y = loop.texts$K[loop.count], paste("p = ", round(mad.test(temp2, use.theory = FALSE)$p.value, 3)), adj = c(1, 0))
    abline(h = 0, lty = 2, col = "red", lwd = 0.5)
    
    temp3 <- envelope(i, linearpcfcross.inhom.avrg, j[1], j[2], sigma = 100, r = seq(0, 350, 10), simulate = temp, nsim = 199, nrank = 5, verbose = FALSE, savefuns = TRUE)
    temp4 <- envelope(i, linearpcfcross.inhom.avrg, j[1], j[2], sigma = 100, r = seq(0, 350, 10), simulate = temp2, nsim = 199, nrank = 5, verbose = FALSE, savefuns = TRUE)
    
    plot(temp4, cbind(obs - mmean + 1, lo = lo - mmean + 1, hi = hi - mmean + 1) ~ r, main = "", legend = FALSE, cex.axis = 0.8, shadecol = "#00000030")
    plot(temp3, cbind(obs - mmean + 1, lo = lo - mmean + 1, hi = hi - mmean + 1) ~ r, main = "", legend = FALSE, cex.axis = 0.8, col = "#996600", shadecol = "#99660030", add = TRUE)
    text(350, y = loop.texts$PCF[loop.count], paste("p = ", round(mad.test(temp3, use.theory = FALSE)$p.value, 3)), col = "#996600", adj = c(1, 1.5))
    text(350, y = loop.texts$PCF[loop.count], paste("p = ", round(mad.test(temp4, use.theory = FALSE)$p.value, 3)), adj = c(1, 0))
    abline(h = 1, lty = 2, col = "red", lwd = 0.5)
  }
}

Figure 10: Same as Fig. 5 but along the segment of dirt road sampled in 2006.

Both parametric and non-parametric Poisson spatial models with inhomogeneous intensities fit the observed patterns much better when they pose species-specific parameters of variation along the road (Models 2) than those with a common parameter defining the change of intensity (Models 1):

#con modelo paramétricos (polinomio)
#rescale es para pasar escala a km y evitar error de matriz singular; no cambia nada
cat("Parametric models (cubic polynomials)\n")
anova(lppm(rescale(Torre2006_Pogos, 1000) ~ polynom(x, 3) + marks), lppm(rescale(Torre2006_Pogos, 1000) ~ polynom(x, 3) * marks), test = "LRT")

#con modelo no paramétrico (cubic spline)
cat("/nNon-parametric models (cubic splines, df = 5)\n")
anova(lppm(Torre2006_Pogos ~ splines::bs(x, df = 5) + marks), lppm(Torre2006_Pogos ~ splines::bs(x, df = 5) * marks), test = "LRT")

## Parametric models (cubic polynomials)
## Analysis of Deviance Table
## 
## Model 1: ~x + I(x^2) + I(x^3) + marks     Poisson
## Model 2: ~(x + I(x^2) + I(x^3)) * marks   Poisson
##   Npar Df Deviance       Pr(>Chi)    
## 1    6                               
## 2   12  6   50.752 0.000000003321 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## /nNon-parametric models (cubic splines, df = 5)
## Analysis of Deviance Table
## 
## Model 1: ~splines::bs(x, df = 5) + marks      Poisson
## Model 2: ~splines::bs(x, df = 5) * marks      Poisson
##   Npar Df Deviance  Pr(>Chi)    
## 1    8                          
## 2   18 10   88.827 9.146e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Association between Pogonomyrmex spp. and soil characteristics

Soil permeability

The quick field technique used to roughly estimate soil permeability through the time it takes a volume of water to infiltrate through the soil surface, seems able to differentiate the most typical plant communities in Ñacuñán Reserve: algarrobal (an open woodland of scattered algarrobos within a shrub matrix), jarillal (shrubland) and medanal (vegetated sand dunes)(Fig. 11). Although our soil sampling was concentrated on a single spot per habitat (20 samples 10–15 m apart), infiltration times were distinct enough to be used as diagnostic (threshold midway between medanal and algarrobal means: ~60 sec; between algarrobal and jarillal: ~284 sec).

We defined a Soil Permeability Index (SPI) as an inverse function of the logarithm of the infiltration time measured at field (i.e., high infiltration time = low permeability, and viceversa) after substracting log(2 min) so it is centered around typical values for the algarrobal (mean SPI: algarrobal = -0.36, medanal = 1.75, jarillal = -1.37). The proposed values to differentiate among habitats, then, correspond to SPI = -0.86 and SPI = 0.7.

grid.arrange(
  # A. tiempo de infiltración por ambiente
  ggplot(subset(Perm_picadas, Lugar != "LaFad" & Lugar != "Doble"), aes(y = Lugar, x = log(Tiempo), fill = Lugar)) +
  geom_boxplot(coef = 3) +
  scale_fill_manual(values = c("brown4", "darkgreen", "orange")) +
  labs(x = "log(Infiltration Time [sec])", y = "Habitat") +
  scale_y_discrete(limits = c("Jarillal","Algarrobal","Medanal")) +
  geom_vline(xintercept = c(4.09, 5.65), colour = "blue", linewidth = 0.5, linetype = 2) +
  tema,
  # B. frecuencia de muestras según tiempo de infiltración
  ggplot(subset(Perm_picadas, Lugar != "LaFad" & Lugar != "Doble"), aes(log(Tiempo), fill = Lugar)) +
  geom_histogram(position = "identity", bins = 15, alpha = 0.8) +
  scale_fill_manual(values = c("brown4", "darkgreen", "orange")) +
  labs(x = "log(Infiltration Time [sec])", y = "Samples") +
  ylim(c(0, 11)) +
  geom_vline(xintercept = c(4.09, 5.65), colour = "blue", linewidth = 0.5, linetype = 2) +
  tema,
layout_matrix = rbind(c(1, 1, 1, 1, 1, 2, 2, 2, 2))
)

Figure 11: Soil permeability, estimated by the infiltration time of a constant volume of water through the soil surface, in three typical habitats of the central Monte desert: algarrobal, jarillal and medanal (20 samples each). Vertical dotted lines show possible threshold values, midway between the estimated habitat means.

Infiltration times along the roads showed values within those obtained in the sampled natural habitats (Fig. 12), though they showed important variations in short distances, even between samples taken ~4 m apart towards both sides of the road centers. The bold measuring technique employed probably adds to the naturally varying factors, the somewhat artificial height profiles of the roads and the soil being altered and transported by transit and heavy machinery during periods of road maintenance. The relationship between the two measurements taken at the same road mark was closer in Doble that in La Fadina (lineal regression: R²: 0.33 and 0.15, respectively), where the side of the road farther from a parallel wire fence (left) showed higher infiltration times (in average, 46 sec higher) than the other (Fig. 12).

In spite of variation between close samples, predicted permeability indices (obtained by fitting cubic smoothing splines to the observed values at both sides; df = 35) changed notably along both dirt roads, with sectors that reach soil permeability values close to those of jarillal (lower) or to medanal (higher) (Fig. 12), particularly in Doble.

#función para elegir el nivel de suavizado óptimo (best AIC) de permeabilidad según su capacidad predictiva de la intensidad del lppm (modelo perm * marks)
BestSmoothDF <- function(LppY, LppX) {
  Npoints <- npoints(LppX)/2
  Best <- 1:Npoints
  for(j in 2:Npoints) {
    Smoothed <- funxy(function(x, y){predict(with(as.data.frame(LppX), smooth.spline(marks ~ x, df = j)), x)$y}, W = Window(LppY))
    Best[j] <- AIC(lppm(LppY ~ marks * Smoothed))
  }
  Best[1] <- max(Best)
  plot(1:Npoints, Best, type = "l", main = paste(match.call()$LppY, "~ marks *", match.call()$LppX))
  abline(h = min(Best), lty = 2, col = "red", lwd = 0.5)
  abline(h = min(Best) + 2, lty = 2, col = "red", lwd = 0.5)
  text(x = 20, y = max(Best) - 5, paste(which(Best == min(Best)), "[", paste(range(which(Best <= min(Best) + 2)), collapse = "–"), "]"))
}

# Dado que los valores se tomaron a distancias particulares de cada segmento muestreado, es necesario hacer una interpolación para estimar los valores de permeabilidad en todas las distancias intermedias. El nivel de suavizado óptimo de acuerdo a su valor predictivo de la variación de la densidad de cada especie a lo largo de las transectas en los modelos Poisson posteriores se puede ver en los siguientes gráficos (df óptimo en la función _spline_ de suavizado, con el rango indistinguible ΔAIC < 2 entre corchetes):

par(mar = c(2, 2, 3, .5))
layout(matrix(1:6, 3, 2, byrow = FALSE))
BestSmoothDF(LaFad2004_Pogos, LaFad_Perm)
BestSmoothDF(LaFad2007_Pogos, LaFad_Perm)
BestSmoothDF(LaFad2010_Pogos, LaFad_Perm)
BestSmoothDF(Doble2004_Pogos, Doble_Perm)
BestSmoothDF(Doble2007_Pogos, Doble_Perm)
BestSmoothDF(Doble2010_Pogos, Doble_Perm)

#se decide usar df=35, que es un valor cercano al óptimo (y siempre dentro del rango indistinguible del óptimo) en todos los casos

grid.arrange(
  # A. LaFad: valores de permeabilidad obtenidos y suavizado + interpolación (smooth.spline con df elegido mediante BestSmoothDF según gráfico anterior)
  Perm_picadas %>%
    filter(Lugar == "LaFad") %>%
    mutate(Sm_all = predict(smooth.spline(SPI ~ Dist, df = 35, keep.data = FALSE), Dist)$y,
           Sm_der = predict(smooth.spline(SPI[Lado == "der"] ~ Dist[Lado == "der"], df = 35, keep.data = FALSE), Dist)$y,
           Sm_izq = predict(smooth.spline(SPI[Lado == "izq"] ~ Dist[Lado == "izq"], df = 35, keep.data = FALSE), Dist)$y) %>%
    ggplot(aes(x = Dist, y = SPI, colour = Lado)) +
      geom_hline(yintercept = c(0.7, -0.86), linetype = 2, colour = "blue", linewidth = 0.5) +
      geom_point() +
      geom_line(aes(y = Sm_der), colour = "green4") +
      geom_line(aes(y = Sm_izq), colour = "red") +
      geom_line(aes(y = Sm_all), colour = "black", linewidth = 1) +
      scale_colour_manual(values = c("green4", "red")) +
      labs(y = "SPI", title = "LaFad") +
      tema,
  # LaFad: relación entre valores medidos a derecha e izquierda
  Perm_picadas %>%
  filter(Lugar == "LaFad") %>%
  pivot_wider(id_cols = Dist, names_from = Lado, values_from = SPI) %>%
  ggplot(aes(x = izq, y = der)) +
    geom_hline(yintercept = c(0.7, -0.86), linetype = 2, colour = "blue", linewidth = 0.5) +
    geom_vline(xintercept = c(0.7, -0.86), linetype = 2, colour = "blue", linewidth = 0.5) +
    geom_point() +
    geom_abline(aes(intercept = 0, slope = 1), linetype = "dashed") +
    labs(x = "SPI (left)", y = "SPI (right)") +
    scale_x_continuous(limits = c(-2.3, 2)) +
    scale_y_continuous(limits = c(-2.3, 2)) +
    tema,
  # B. Doble: valores de permeabilidad obtenidos y suavizado + interpolación (smooth.spline con df elegido mediante BestSmoothDF según gráfico anterior)
  Perm_picadas %>%
    filter(Lugar == "Doble") %>%
    mutate(Sm_all = predict(smooth.spline(SPI ~ Dist, df = 35, keep.data = FALSE), Dist)$y,
           Sm_der = predict(smooth.spline(SPI[Lado == "der"] ~ Dist[Lado == "der"], df = 35, keep.data = FALSE), Dist)$y,
           Sm_izq = predict(smooth.spline(SPI[Lado == "izq"] ~ Dist[Lado == "izq"], df = 35, keep.data = FALSE), Dist)$y) %>%
    ggplot(aes(x = Dist, y = SPI, colour = Lado)) +
      geom_hline(yintercept = c(0.7, -0.86), linetype = 2, colour = "blue", linewidth = 0.5) +
      geom_point() +
      geom_line(aes(y = Sm_der), colour = "green4") +
      geom_line(aes(y = Sm_izq), colour = "red") +
      geom_line(aes(y = Sm_all), colour = "black", linewidth = 1) +
      scale_colour_manual(values = c("green4", "red")) +
      labs(y = "SPI", title = "Doble") +
      tema,
  # Doble: relación entre valores medidos a derecha e izquierda
  Perm_picadas %>%
  filter(Lugar == "Doble") %>%
  pivot_wider(id_cols = Dist, names_from = Lado, values_from = SPI) %>%
  ggplot(aes(x = izq, y = der)) +
    geom_hline(yintercept = c(0.7, -0.86), linetype = 2, colour = "blue", linewidth = 0.5) +
    geom_vline(xintercept = c(0.7, -0.86), linetype = 2, colour = "blue", linewidth = 0.5) +
    geom_point() +
    geom_abline(aes(intercept = 0, slope = 1), linetype = "dashed") +
    labs(x = "SPI (left)", y = "SPI (right)") +
    scale_x_continuous(limits = c(-2.3, 2)) +
    scale_y_continuous(limits = c(-2.3, 2)) +
    tema,
layout_matrix = rbind(c(1, 1, 1, 2), c(3, 3, 3, 4))
)

Figure 12: Soil permeability along two dirt roads (La Fadina, above, and Doble, below) crossing the algarrobal in Ñacuñán Reserve, central Monte desert. Permeability was estimated by the time (log(seconds)) it took a fixed amount of water to completely infiltrate through the soil surface (see details in Methods) in two spots ~4 m apart (red: left; green: right) every 25 m along each 1.5 km segment of road. Continuous lines show the smoothed profile (fitting cubic splines with df = 35) of log-infiltration times (red: left side samples; green: right side samples; black: all samples); discontinuous blue lines show the proposed thresholds between characteristic infiltration times in three habitats (see Fig. 11).

# pasado a funcion funxy para la predicción en cualquier "x" 
LaFad_Perm_Smooth <- funxy(function(x, y){predict(with(as.data.frame(LaFad_Perm), smooth.spline(marks ~ x, df = 35)), x)$y}, W = Ventana)
Doble_Perm_Smooth <- funxy(function(x, y){predict(with(as.data.frame(Doble_Perm), smooth.spline(marks ~ x, df = 35)), x)$y}, W = Ventana)
Both_Perm_Smooth <- list(LaFad_Perm_Smooth, Doble_Perm_Smooth)

The estimated intensities of heterogeneous Poisson spatial processes for each Pogonomyrmex species show a strong relationship with soil permeability: PM associates with high permeability (< 135 seg) in both segments, PI with low permeability (> 180 seg), particularly in Doble, and PR showed more variability, though avoiding the lowest values in La Fadina (Fig. 13).

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

#estimación de intensidad de cada sp según permeabilidad
#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("LaFad 2004: PI", "LaFad 2004: PM", "LaFad 2004: PR", "Doble 2004: PI", "Doble 2004: PM", "Doble 2004: PR"), Lambda = c(0.15, 0.04, 0.06, 0.15, 0.05, 0.04))
loop.count <- 0

for(i in 1:length(Picadas2004)) {
  for(j in split(Picadas2004[[i]])) {
    loop.count <- loop.count + 1
    plot(rhohat(j, Both_Perm_Smooth[[i]]), legend = FALSE, main = loop.texts$Main[loop.count], cex.axis = .8)
    text(0.65, y = loop.texts$Lambda[loop.count], paste("p = ", round(cdf.test(j, Both_Perm_Smooth[[i]], test = "ks", jitter = F)$p.value, 3)), adj = c(1, 1))
    abline(v = 0.7, lty = 2, col = "blue", lwd = 0.5)
    abline(v = -0.86, lty = 2, col = "blue", lwd = 0.5)
  }
}

Figure 13: Resource selection functions, or intensities of heterogeneous Poisson linear point process models for the active colonies of each Pogonomyrmex spp. in 2004 (black line: fixed-bandwith kernel density estimator; grey areas: pointwise confidence intervals) as functions of (log) infiltration times [sec] along the two segments of dirt roads (left: La Fadina; right: Doble). Discontinuous vertical lines show the proposed infiltration time thresholds between three typical habitats (see Fig. 11). P values result from Kolmogorov-Smirnov goodness-of-fit tests against (homogeneous = no relationship) CSR models with same mean intensities.

Similar conclusions can be derived from observing the estimated soil permeability at every position along the dirt roads occupied by active nests of each Pogonomyrmex spp.: PI associated with lower permeability (particularly in Doble), PM with the highest, and PR, more variable, with middle values (Fig. 14).

#datos de permeabilidad estimado para cada colonia por especie para las dos picadas
data.frame(Perm = c(LaFad_Perm_Smooth(x = LaFad2004_Pogos$data$x, y = LaFad2004_Pogos$data$y),
                    Doble_Perm_Smooth(x = Doble2004_Pogos$data$x, y = Doble2004_Pogos$data$y)),
             Sp = c(LaFad2004_Pogos$data$marks, 
                    Doble2004_Pogos$data$marks), 
         Picada = c(rep("LaFad", dim(LaFad2004_Pogos$data)[1]), 
                    rep("Doble", dim(Doble2004_Pogos$data)[1]))) %>%
  ggplot(aes(x = Sp, y = Perm, fill = Sp, colour = Sp)) +
    geom_boxplot() +
    geom_point() +
    facet_wrap(vars(Picada)) +
    scale_colour_manual(values = c("darkgreen", "darkblue", "darkred")) +
    scale_fill_manual(values = alpha(c("darkgreen", "darkblue", "darkred"), 0.5)) +
    labs(x = "Species", y = "Soil permeability (SPI)") +
    geom_hline(yintercept = c(0.7, -0.86), lty = 2, col = "blue", lwd = 0.5) +
    tema

Figure 14: Distributions (median and quartiles) of infiltration times (log(sec)) associated to each of the active nests (dots) of each Pogonomyrmex species detected in 2004 in the 1.5 km segment sampled in each dirt road (La Fadina and Doble) crossing the algarrobal in Ñacuñán, central Monte desert.

Topography

The general slope of these roads, as in most of the area, is very gentle (La Fadina: 0.26%; Doble: 0.22%). Still, both of them (particularly Doble) show several sectors of relative small peaks and valleys of a few decameters above and below the general slope (Fig. 15).

#Para calcular la altura relativa (o TPI) de un punto de la transecta es necesario determinar la altura de todos los puntos sobre la transecta (interpolando, dado que la altura sobre el nivel del mar se midió en distancias particulares) y la altura promedio de sus entornos. Lo primero se puede obtener con una función de suavizado que resulte "cercana" a cada punto medido (i.e., con d.f. elevado; en este caso, de la mitad del número de datos). La estimación de la altura promedio exige determinar un tamaño de ventana para calcular una media móvil o un suavizado con menos grados de libertad. El nivel de suavizado óptimo para esa altura promedio, de acuerdo a su valor predictivo de la variación de la densidad de cada especie a lo largo de las transectas en los modelos Poisson posteriores, se puede ver en los gráficos siguientes (df óptimo en la función _spline_ de suavizado para la altura promedio, rango indistinguible ΔAIC < 2 entre corchetes)

#función para elegir los niveles de suavizado óptimos (best AIC) para obtener la altura relativa según su capacidad predictiva de la intensidad del lpp (modelo perm * marks)
BestSmoothedTPI <- function(LppY, LppX) {
  Npoints <- npoints(LppX)
  Best <- 1:(Npoints/2)
  for(j in 2:(Npoints/2)) {
    SmoothedTPI <- funxy(function(x, y){predict(with(as.data.frame(LppX), smooth.spline(marks ~ x, df = Npoints/2)), x)$y - predict(with(as.data.frame(LppX), smooth.spline(marks ~ x, df = j)), x)$y}, W = Window(LppY))
      Best[j] <- AIC(lppm(LppY ~ marks * SmoothedTPI))
  }
  Best[1] <- max(Best, na.rm = TRUE)
  plot(1:(Npoints/2), Best, type = "l", main = paste(match.call()$LppY, "~ marks *", match.call()$LppX))
  abline(h = min(Best), lty = 2, col = "red", lwd = 0.5)
  abline(h = min(Best) + 2, lty = 2, col = "red", lwd = 0.5)
  text(x = Npoints/4, y = max(Best) - 5, paste(which(Best == min(Best)), "[", paste(range(which(Best <= min(Best) + 2)), collapse = "–"), "]"))
}

par(mar = c(2, 2, 3, .5))
layout(matrix(1:6, 3, 2, byrow = FALSE))
BestSmoothedTPI(LaFad2004_Pogos, LaFad_Alt)
BestSmoothedTPI(LaFad2007_Pogos, LaFad_Alt)
BestSmoothedTPI(LaFad2010_Pogos, LaFad_Alt)
BestSmoothedTPI(Doble2004_Pogos, Doble_Alt)
BestSmoothedTPI(Doble2007_Pogos, Doble_Alt)
BestSmoothedTPI(Doble2010_Pogos, Doble_Alt)

# se elige df = (muestras/2) para describir la altura en cada punto
# se elige df = 3 como mejor valor para describir la pendiente general en ambos casos

#valores de altura obtenidos y suavizado/interpolación para generar el DEM [Digital Elevation Model] con predicciones de un smooth.spline con df = length(x)/2; precisión horizontal = 1 m
LaFad_Alt_Smoothed <- funxy(function(x , y) {predict(with(as.data.frame(LaFad_Alt), smooth.spline(marks ~ x, df = length(x)/2)), x)$y}, W = Ventana)
Doble_Alt_Smoothed <- funxy(function(x , y) {predict(with(as.data.frame(Doble_Alt), smooth.spline(marks ~ x, df = length(x)/2)), x)$y}, W = Ventana)

# si hacen falta indicadores que usen pendientes se podrían obtener con los valores de la primera derivada directo de smooth.spline
#[tipo CTI (Compound Topographic Index): measure of both slope and the upstream contributing area, o TWI tendency of an area to accumulate water; higher = wetter) TWI=ln(SCA/tanφ) SCA: specific catchment area and φ is the slope angle (radians), ¿pero cómo calcular la cuenca de cada punto?)] 
#LaFad_Alt_Slope <- funxy(function(x , y) {predict(with(as.data.frame(LaFad_Alt), smooth.spline(marks ~ x, df = length(x)/2)), x, deriv = 1)$y}, W = Ventana)
#Doble_Alt_Slope <- funxy(function(x , y) {predict(with(as.data.frame(Doble_Alt), smooth.spline(marks ~ x, df = length(x)/2)), x, deriv = 1)$y}, W = Ventana)

#tpi: diferencia entre Alt[DEM] (smooth con df = x/2) y ~rollingmean(Alt[DEM]) (smooth con df = 3)
LaFad_Alt_TPI <- funxy(function(x , y) {predict(with(as.data.frame(LaFad_Alt), smooth.spline(marks ~ x, df = length(x)/2)), x)$y - predict(with(as.data.frame(LaFad_Alt), smooth.spline(marks ~ x, df = 3)), x)$y}, W = Ventana)
Doble_Alt_TPI <- funxy(function(x , y) {predict(with(as.data.frame(Doble_Alt), smooth.spline(marks ~ x, df = length(x)/2)), x)$y - predict(with(as.data.frame(Doble_Alt), smooth.spline(marks ~ x, df = 3)), x)$y}, W = Ventana)
Both_Alt_TPI <- list(LaFad_Alt_TPI, Doble_Alt_TPI)

par(mar = c(1, 4, 1.5, 0.5))
layout(matrix(1:2, 2, 1), heights = c(2, 1))
#LaFad
with(as.data.frame(LaFad_Alt), plot(marks ~ x, ylab = "Altitude (masl)", main = "La Fadina"))
lines(0:1500, LaFad_Alt_Smoothed(0:1500, 0), col ="red")
lines(0:1500, predict(with(as.data.frame(LaFad_Alt), smooth.spline(marks ~ x, df = 3)), 0:1500)$y, col = "blue")
with(as.data.frame(LaFad_Alt), plot(LaFad_Alt_TPI(x, 0) * 100 ~ x, type = "l", ylab = "TPI (cm)"))
abline(h = 0, lty = 2, lwd = 0.5)
#with(as.data.frame(LaFad_Alt), plot(LaFad_Alt_Slope(x, 0) ~ x, ylab = "Slope (%)"))
#lines(0:1500, LaFad_Alt_Slope(0:1500, 0))
#abline(h = 0, lty = 2, lwd = 0.5)

#Doble
with(as.data.frame(Doble_Alt), plot(marks ~ x, ylab = "Altitude (masl)", main = "Doble"))
lines(0:1500, Doble_Alt_Smoothed(0:1500, 0), col = "red")
lines(0:1500, predict(with(as.data.frame(Doble_Alt), smooth.spline(marks ~ x, df = 3)), 0:1500)$y, col = "blue")
with(as.data.frame(Doble_Alt), plot(Doble_Alt_TPI(x, 0) * 100 ~ x, type = "l", ylab = "TPI (cm)"))

abline(h = 0, lty = 2, lwd = 0.5)
#with(as.data.frame(Doble_Alt), plot(Doble_Alt_Slope(x, 0) ~ x, ylab = "Slope (%)"))
#lines(0:1500, Doble_Alt_Slope(0:1500, 0))
#abline(h = 0, lty = 2, lwd = 0.5)

Figure 15: Altitude (meters above sea level) and TPI (topographic position index) along each of the sampled 1.5 km segments of dirt roads (La Fadina, above, and Doble, below) crossing the algarrobal in Ñacuñán Reserve, central Monte desert. Dots: sampled altitude values; red line: estimation of local altitude (cubic smoothing spline with df = number of samples/2); blue line: moving average altitude (spline with df = 3); black line: estimated topographic index, or the difference between local altitude and the average altitude in the surroundings.

The estimated intensities of the Poisson processes associated with each Pogonomyrmex species active colonies show a similar relationship with the peaks and valleys of local topography than was detected with soil permeability: PM associates with relative high areas, PI with low areas (though only markedly in Doble) and PR with both average and relatively low areas (Fig. 16).

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

#estimación de intensidad de cada sp según altura relativa

#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("LaFad 2004: PI", "LaFad 2004: PM", "LaFad 2004: PR", "Doble 2004: PI", "Doble 2004: PM", "Doble 2004: PR"), Lambda = c(0.07, 0.045, 0.10, 0.06, 0.035, 0.045))
loop.count <- 0

for(i in 1:length(Picadas2004)) {
  for(j in split(Picadas2004[[i]])) {
    loop.count <- loop.count + 1
    plot(rhohat(j, Both_Alt_TPI[[i]]), legend = FALSE, main = loop.texts$Main[loop.count], cex.axis = .8)
    text(0.17 * (1 + floor(loop.count/4)), y = loop.texts$Lambda[loop.count], paste("p = ", round(cdf.test(j, Both_Alt_TPI[[i]], test = "ks", jitter = F)$p.value, 3)), adj = c(1, 1))
    abline(v = 0, lty = 2, col = "blue", lwd = 0.5)
  }
}

Figure 16: Same as Fig. 13 but as function of the estimated TPI along both dirt road segments.

Again, something similar can be observed when plotting the TPI associated with each active nest detected: PI associates with relative depressions in Doble, PM with the highest areas, and PR is rather variable (Fig. 17).

#datos de altitud relativa estimados para cada colonia por especie para las dos picadas
data.frame(AltTPI = c(LaFad_Alt_TPI(x = LaFad2004_Pogos$data$x, y = LaFad2004_Pogos$data$y),
                    Doble_Alt_TPI(x = Doble2004_Pogos$data$x, y = Doble2004_Pogos$data$y)),
             Sp = c(LaFad2004_Pogos$data$marks, 
                    Doble2004_Pogos$data$marks), 
         Picada = c(rep("LaFad", dim(LaFad2004_Pogos$data)[1]), 
                    rep("Doble", dim(Doble2004_Pogos$data)[1]))) %>%
  ggplot(aes(x = Sp, y = AltTPI, fill = Sp, colour = Sp)) +
  geom_boxplot() +
  geom_point() +
  scale_colour_manual(values = c("darkgreen", "darkblue", "darkred")) +
  scale_fill_manual(values = alpha(c("darkgreen", "darkblue", "darkred"), 0.5)) +
  labs(y = "TPI [m]", x = "Species") +
  geom_hline(yintercept = 0, colour = "blue", linewidth = 0.5, linetype = 2) +
  facet_wrap(vars(Picada)) +
  tema

Figure 17: Distributions (median and quartiles) of estimated TPI values associated to each of the active nests (dots) of each Pogonomyrmex species detected in 2004 in the 1.5 km segment sampled in each dirt road (La Fadina and Doble) crossing the algarrobal in Ñacuñán, central Monte desert.

Both soil characteristics

Both environmental variables measured associated with the variation in observed colony density of Pogonomyrmex species along the two sampled dirt roads (more strictly, they proved relevant single predictors of the intensity of Poisson point pattern processes that influence the spatial locations of the colonies of one or more of the species).

However, although not strongly, both measured variables were correlated (Fig. 18). In Doble, as could have been expected: permeability is lower in relative depressions (where fine-particle, clayish soil tends to accumulate) while highest sectors have the highest permeability values (sandy soil). In La Fadina, instead, the linear correlation is almost none, maybe as consequence of the end of the segment (1400–1500 m) appearing as a depression with high permeability. This may be an edge artifact because of a lack of altitude information out of the segment to properly estimate the average altitude around; the linear tendency becomes negative, though still rather flat, when avoiding the last 150 m; Fig. 18).

layout(matrix(1))
#Perm vs TPI
data.frame(Perm = c(LaFad_Perm_Smooth(seq(0, 1500, 25), 0), Doble_Perm_Smooth(seq(0, 1500, 25), 0)), TPI = c(LaFad_Alt_TPI(seq(0, 1500, 25), 0), Doble_Alt_TPI(seq(0, 1500, 25), 0)), Picada = rep(c("LaFad", "Doble"), each = 1525/25)) %>%
  ggplot(aes(x = TPI, y = Perm, colour = Picada)) +
  geom_smooth(data = data.frame(Perm = LaFad_Perm_Smooth(seq(0, 1350, 25), 0), TPI = LaFad_Alt_TPI(seq(0, 1350, 25), 0)), colour = "deepskyblue", fill = "lightgrey", alpha = 0.5, linetype = 2, method = "lm", formula = y ~ x) +
  geom_point() +
  geom_smooth(method = "lm", formula = y ~ x) +
  labs(y = "Soil permeability (SPI)", x = "TPI [m]") +
  tema

Figure 18: Linear relationship between the two edaphic variables estimated along the two segments of dirt roads (blue: La Fadina; red: Doble) crossing the algarrobal in Ñacuñán, central Monte desert. Dots are estimated values every 25 m of the 1.5 km segments. The dashed line is the linear relationship estimated for the first 1350 m of La Fadina.

As a consequence, both factors seems to be involved in a better separation of the environmental conditions associated with active colonies of each of the Pogonomyrmex species (Fig. 19).

#datos de altitud relativa y permeabilidad estimados para cada colonia por especie para las dos picadas
temp.data <- data.frame(Alt = c(LaFad_Alt_TPI(x = LaFad2004_Pogos$data$x, y = LaFad2004_Pogos$data$y),
                    Doble_Alt_TPI(x = Doble2004_Pogos$data$x, y = Doble2004_Pogos$data$y)),
           Perm = c(LaFad_Perm_Smooth(x = LaFad2004_Pogos$data$x, y = LaFad2004_Pogos$data$y),
                    Doble_Perm_Smooth(x = Doble2004_Pogos$data$x, y = Doble2004_Pogos$data$y)),
             Sp = c(LaFad2004_Pogos$data$marks, 
                    Doble2004_Pogos$data$marks), 
         Picada = c(rep("LaFad", dim(LaFad2004_Pogos$data)[1]), 
                    rep("Doble", dim(Doble2004_Pogos$data)[1])))
grid.arrange(
  ggplot(data.frame(Alt = c(LaFad_Alt_TPI(0:1500, 0), Doble_Alt_TPI(0:1500, 0))), aes (x = Alt)) +
    geom_histogram(bins = 25) +
    coord_cartesian(xlim = c(-0.51, 0.51)) +
    labs(y = "\n") +
    tema + theme(axis.text = element_blank(), axis.title.x = element_blank(), rect = element_blank(), line = element_blank()),
  ggplot() + geom_blank() + theme_void(),
  ggplot(temp.data, aes(x = Alt, y = Perm, fill = Sp, colour = Sp, shape = Picada)) +
    geom_point(size = 2) +
    scale_colour_manual(values = c("darkgreen", "darkblue", "darkred")) +
    scale_fill_manual(values = alpha(c("darkgreen", "darkblue", "darkred"), 0.2)) +
    labs(y = "Soil permeability (SPI)", x = "Topography (TPI [m])") +
    coord_cartesian(xlim = c(-0.51, 0.51), ylim = c(-1.6, 1.2)) +
    geom_hline(yintercept = c(0.7, -0.86), lty = 2, col = "blue", lwd = 0.5) +
    geom_vline(xintercept = 0, colour = "blue", linewidth = 0.5, linetype = 2) +
    geom_rug(sides = "rt") +
    stat_ellipse(aes(group = Sp), type = "norm", level = 0.9) +
    scale_shape_manual(values = c(21, 19)) +
    tema,
  ggplot(data.frame(Perm = c(LaFad_Perm_Smooth(0:1500, 0), Doble_Perm_Smooth(0:1500, 0))), aes(y = Perm)) +
    geom_histogram(bins = 25) +
    coord_cartesian(ylim = c(-1.6, 1.2)) +
    labs(x = "\n") +
    tema + theme(axis.text = element_blank(), axis.title.y = element_blank(), rect = element_blank(), line = element_blank()),
  nrow = 2, widths = c(4, 1), heights = c(1, 4))

Figure 19: Relative topography (TPI) and soil permeability (SPI) values associated with active nests of each Pogonomyrmex species (PI: green, PM: blue, PR: red) detected in 2004 in the 1.5 km segment sampled in each dirt road (La Fadina: full dots; Doble: open dots) crossing the algarrobal in Ñacuñán, central Monte desert, with data ellipsoids for each species at the 90% levels assuming multivariate normal distributions. Histograms (above and right) show the distribution of values of each variable estimated every meter along both road segments.

Poisson point pattern models

— La Fadina

Through the following series of comparisons among increasingly complex (=forward) nested Poisson models, we see that in the segment of La Fadina dirt road, [1] the three Pogonomyrmex species differ in their average densities, [2] that the edaphic variables are not very relevant single linear or quadratic predictors of the variation of Pogonomyrmex spp. active colonies as a whole (though soil permeability is a bit better) but [3] both are strong predictors when their effect is allowed to be species-specific [= interaction] (again, with soil permeability a bit more relevant than topography).

In spite of both edaphic variables being partially related (Fig. 18), [4] topography adds to the predictive power of the model already containing soil permeability, though [5] the interaction among them is not that important.

#comparación de modelos Poisson anidados
#(1) comparación modelo homogeneo por especie (cada sp con su intensidad constante en el espacio) vs. intensidad única (misma intensidad homogénea de las 3 spp)
cat("[1] Homogeneous models with a single (Model 1) or species-specific (Model 2) intensities\n")
anova(lppm(LaFad2004_Pogos ~ 1),
      lppm(LaFad2004_Pogos ~ marks),
      test = "LR")

#(2) agregando a la permeabilidad o a la altura como predictores lineales o cuadráticos independientes de la variación espacial de la intensidad  (i.e., afectando a las 3 spp por igual) 
cat("\n[2] Model with species-specific spatially homogeneous intensities (Models 1) or with their intensities varying (in the same way) linearly (Models 2) or quadratically (Models 3) with either of the edaphic variables\n")
anova(lppm(LaFad2004_Pogos ~ marks),
      lppm(LaFad2004_Pogos ~ marks + LaFad_Perm_Smooth),
      lppm(LaFad2004_Pogos ~ marks + polynom(LaFad_Perm_Smooth, 2)),
      test = "LR")
anova(lppm(LaFad2004_Pogos ~ marks),
      lppm(LaFad2004_Pogos ~ marks + LaFad_Alt_TPI),
      lppm(LaFad2004_Pogos ~ marks + polynom(LaFad_Alt_TPI, 2)),
      test = "LR")

#(3) agregando a la permeabilidad o a la altura como predictores lineales de la variación espacial de la intensidad  (i.e., afectando a las 3 spp por igual) o con interacción entre edáficas y especies [= test de interacción], de modo lineal o cuadrático
cat("\n[3] Model with heterogeneous intensities varying in the same way linearly with each edaphic factor (Models 1) or with species-specific linear (Models 2) or quadratic (Models 3) relationships\n")
anova(lppm(LaFad2004_Pogos ~ marks + LaFad_Perm_Smooth),
      lppm(LaFad2004_Pogos ~ marks * LaFad_Perm_Smooth),
      lppm(LaFad2004_Pogos ~ marks * polynom(LaFad_Perm_Smooth, 2)),
      test = "LR")
anova(lppm(LaFad2004_Pogos ~ marks + LaFad_Alt_TPI),      
      lppm(LaFad2004_Pogos ~ marks * LaFad_Alt_TPI),
      lppm(LaFad2004_Pogos ~ marks * polynom(LaFad_Alt_TPI, 2)),
      test = "LR")

#(4) agregando a la topografía sobre el modelo que ya tiene a permeabilidad * especie
cat("\n[4] Model with species-specific quadratic effects of soil permeability (Model 1), with an additional independent quadratic effect of topography (Model 2), or with both species-specific effects interacting (Model 3)\n")
anova(lppm(LaFad2004_Pogos ~ marks * polynom(LaFad_Perm_Smooth, 2)),
      lppm(LaFad2004_Pogos ~ marks * polynom(LaFad_Perm_Smooth, 2) + marks * polynom(LaFad_Alt_TPI, 2)),
      lppm(LaFad2004_Pogos ~ marks * polynom(LaFad_Perm_Smooth, 2) * polynom(LaFad_Alt_TPI, 2)),
      test = "LR")

#(5) cuadráticas independientes o lineales interactuando
cat("\n[5] Independent quadratic relationships or interacting linear effects?\n")
AIC(lppm(LaFad2004_Pogos ~ marks * LaFad_Perm_Smooth * LaFad_Alt_TPI), lppm(LaFad2004_Pogos ~ marks * polynom(LaFad_Perm_Smooth, 2) + marks * polynom(LaFad_Alt_TPI, 2)))

# single predictor, best lineal model and best quadratic model
PoissonPermLaFad2004 <- lppm(LaFad2004_Pogos ~ marks * polynom(LaFad_Perm_Smooth, 2))
PoissonAltLaFad2004 <- lppm(LaFad2004_Pogos ~ marks * polynom(LaFad_Alt_TPI, 2))
PoissonLaFad2004 <- lppm(LaFad2004_Pogos ~ marks * LaFad_Perm_Smooth * LaFad_Alt_TPI)
Poisson2LaFad2004 <- lppm(LaFad2004_Pogos ~ marks * polynom(LaFad_Perm_Smooth, 2) + marks * polynom(LaFad_Alt_TPI, 2))

## [1] Homogeneous models with a single (Model 1) or species-specific (Model 2) intensities
## Analysis of Deviance Table
## 
## Model 1: ~1   Poisson
## Model 2: ~marks   Poisson
##   Npar Df Deviance Pr(>Chi)   
## 1    1                        
## 2    3  2   11.098 0.003892 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## [2] Model with species-specific spatially homogeneous intensities (Models 1) or with their intensities varying (in the same way) linearly (Models 2) or quadratically (Models 3) with either of the edaphic variables
## Analysis of Deviance Table
## 
## Model 1: ~marks   Poisson
## Model 2: ~marks + LaFad_Perm_Smooth   Poisson
## Model 3: ~marks + (LaFad_Perm_Smooth + I(LaFad_Perm_Smooth^2))    Poisson
##   Npar Df Deviance Pr(>Chi)  
## 1    3                       
## 2    4  1   4.5565  0.03279 *
## 3    5  1   1.0867  0.29720  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Analysis of Deviance Table
## 
## Model 1: ~marks   Poisson
## Model 2: ~marks + LaFad_Alt_TPI   Poisson
## Model 3: ~marks + (LaFad_Alt_TPI + I(LaFad_Alt_TPI^2))    Poisson
##   Npar Df Deviance Pr(>Chi)
## 1    3                     
## 2    4  1  2.42849   0.1191
## 3    5  1  0.93034   0.3348
## 
## [3] Model with heterogeneous intensities varying in the same way linearly with each edaphic factor (Models 1) or with species-specific linear (Models 2) or quadratic (Models 3) relationships
## Analysis of Deviance Table
## 
## Model 1: ~marks + LaFad_Perm_Smooth   Poisson
## Model 2: ~marks * LaFad_Perm_Smooth   Poisson
## Model 3: ~marks * (LaFad_Perm_Smooth + I(LaFad_Perm_Smooth^2))    Poisson
##   Npar Df Deviance   Pr(>Chi)    
## 1    4                           
## 2    6  2   20.410 0.00003699 ***
## 3    9  3   13.772   0.003232 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Analysis of Deviance Table
## 
## Model 1: ~marks + LaFad_Alt_TPI   Poisson
## Model 2: ~marks * LaFad_Alt_TPI   Poisson
## Model 3: ~marks * (LaFad_Alt_TPI + I(LaFad_Alt_TPI^2))    Poisson
##   Npar Df Deviance  Pr(>Chi)    
## 1    4                          
## 2    6  2   18.515 0.0000954 ***
## 3    9  3   11.303    0.0102 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## [4] Model with species-specific quadratic effects of soil permeability (Model 1), with an additional independent quadratic effect of topography (Model 2), or with both species-specific effects interacting (Model 3)
## Analysis of Deviance Table
## 
## Model 1: ~marks * (LaFad_Perm_Smooth + I(LaFad_Perm_Smooth^2))    Poisson
## Model 2: ~marks * (LaFad_Perm_Smooth + I(LaFad_Perm_Smooth^2)) + marks * (LaFad_Alt_TPI + I(LaFad_Alt_TPI^2))     Poisson
## Model 3: ~marks * (LaFad_Perm_Smooth + I(LaFad_Perm_Smooth^2)) * (LaFad_Alt_TPI + I(LaFad_Alt_TPI^2))     Poisson
##   Npar Df Deviance   Pr(>Chi)    
## 1    9                           
## 2   15  6   28.388 0.00007941 ***
## 3   27 12    5.427     0.9422    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## [5] Independent quadratic relationships or interacting linear effects?
##                                                                                                   df
## lppm(LaFad2004_Pogos ~ marks * LaFad_Perm_Smooth * LaFad_Alt_TPI)                                 12
## lppm(LaFad2004_Pogos ~ marks * polynom(LaFad_Perm_Smooth, 2) + marks * polynom(LaFad_Alt_TPI, 2)) 15
##                                                                                                        AIC
## lppm(LaFad2004_Pogos ~ marks * LaFad_Perm_Smooth * LaFad_Alt_TPI)                                 707.1753
## lppm(LaFad2004_Pogos ~ marks * polynom(LaFad_Perm_Smooth, 2) + marks * polynom(LaFad_Alt_TPI, 2)) 699.4964

Predicted nest positions according to all these Poisson spatial model are generally good, but the selected model that includes both edaphic variables seems to improve the prediction of ant colonies in some particular sectors, especially for the two species that show stronger intraspecific aggregation (PI and PM; Fig. 20). Notice that predictions with both variables included are similar but anyway not perfect, particularly for PR (e.g., there are active nest where predicted intensity is well below the average); there is also a sector (near the mark of 1000 m) where PM should be present according to this model but no active nests were detected.

# predicciones vs colonias
par(mar = c(2, 2, 1, .5))
layout(matrix(1:3, 3, 1))

loop.spp <- data.frame(Sp = c("PI", "PM", "PR"), Points = c(0.08, 0.03, 0.06), Cols = c("darkgreen", "darkblue", "darkred"))
graph.count <- 0

for(i in loop.spp$Sp) {
  plot(x = seq(0, 1500, 5), y = predict(Poisson2LaFad2004, locations = data.frame(x = seq(0, 1500, 5), y = 30, marks = i)), type = "l", lwd = 2, main = paste("LaFad2004: ", i))
  lines(x = seq(0, 1500, 5), y = predict(PoissonLaFad2004, locations = data.frame(x = seq(0, 1500, 5), y = 30, marks = i)), col = "darkgrey")
  lines(x = seq(0, 1500, 5), y = predict(PoissonPermLaFad2004, locations = data.frame(x = seq(0, 1500, 5), y = 30, marks = i)), col = "lightblue3")
  lines(x = seq(0, 1500, 5), y = predict(PoissonAltLaFad2004, locations = data.frame(x = seq(0, 1500, 5), y = 30, marks = i)), col = "sienna")
  points(x = LaFad2004_Pogos$data$x[marks(LaFad2004_Pogos) == i], y = rep(loop.spp$Points[loop.spp$Sp == i], length(LaFad2004_Pogos$data$x[marks(LaFad2004_Pogos) == i])), col = loop.spp$Cols[loop.spp$Sp == i])
  abline(h = intensity(LaFad2004_Pogos)[i], lty = 2, col = "red")
}

Figure 20: Intensity predicted by Poisson spatial models based on species-specific parameters for the effects of soil permeability and topographic position [black line: both variables as independent quadratic effects; grey: both variables as interacting linear effects; light blue: only permeability; brown: only topography] on the location of active nests of each Pogonomyrmex ant species in 2004 (circles) along a 1.5 km segment of La Fadina dirt road (Ñacuñán Reserve, central Monte desert). The discontinuous red line shows the average intensity for each species (i.e. corresponding to a CSR model).

The intensity of the spatial process predicted by the best model (i.e., both edaphic variables with per-species quadratic parameters) reflects both the correlation between variables and their partial effects on each Pogonomyrmex species, as previously described (Figs. 21 and 22):

• PI density increases with lower permeability and topographic positions, and decreases with more permeability and higher relative altitudes
  
• In opposition, PM density increases with permeability and relative elevation and is absent from sectors with low permability and relative elevation
  
• The density of PR peaks at mid permeabilities and decreases towards high relative altitudes.

# predicción bivariada en la escala del modelo; para matriz de n×n (=50) valores del rango observado de cada variable en la transecta
temp.data <- data.frame(Alt = rep(seq(min(LaFad_Alt_TPI(0:1500, 0)), max(LaFad_Alt_TPI(0:1500, 0)), length.out = 50), 50), Perm = rep(seq(min(LaFad_Perm_Smooth(0:1500, 0)), max(LaFad_Perm_Smooth(0:1500, 0)), length.out = 50), each = 50)) %>%
  mutate(PI = predict(Poisson2LaFad2004, locations = data.frame(marks = "PI", x = rep(0, length(Alt)), y = rep(0, length(Alt))), covariates = data.frame(LaFad_Alt_TPI = Alt, LaFad_Perm_Smooth = Perm))) %>%
  mutate(PM = predict(Poisson2LaFad2004, locations = data.frame(marks = "PM", x = rep(0, length(Alt)), y = rep(0, length(Alt))), covariates = data.frame(LaFad_Alt_TPI = Alt, LaFad_Perm_Smooth = Perm))) %>%
  mutate(PR = predict(Poisson2LaFad2004, locations = data.frame(marks = "PR", x = rep(0, length(Alt)), y = rep(0, length(Alt))), covariates = data.frame(LaFad_Alt_TPI = Alt, LaFad_Perm_Smooth = Perm))) %>% pivot_longer(., cols = c("PI", "PM", "PR"), names_to = "Sp", values_to = "Lambda") # para un solo ggplot

ggplot(temp.data, aes(x = Alt, y = Perm, fill = log(Lambda))) +
  geom_raster() + facet_wrap(.~Sp) +
  scale_fill_gradient(name = "log(λ)", high = "black", low = "#FFFFFF00") +
  labs(x = "TPI [m]", y = "Soil Permeability (SPI)") +
  geom_hline(yintercept = c(0.7, -0.86), colour = "blue", linewidth = 0.5, linetype = 2) +
  geom_vline(xintercept = 0, colour = "blue", linewidth = 0.5, linetype = 2) +
  tema +
  theme(legend.position = "right")

Figure 21: Log-intensity (i.e. in the scale of the model) predicted by the best Poisson spatial model based on species-specific parameters for the effects of soil permeability and topographic position on the location of active nests of each Pogonomyrmex ant species in 2004 along a 1.5 km segment of La Fadina dirt road (Ñacuñán Reserve, central Monte desert). Discontinuous blue lines show the proposed infiltration time thresholds between three typical habitats (see Fig. 11) and the neutral topography position (TPI = 0, where local altitude matches the average altitude in the surroundings.

#para hacer gráficos separados eliminar pivot_longer y llamar a cada variable en un gráfico
#grid.arrange(ggplot(temp.data, aes(x = Alt, y = Perm, fill = log(PI), alpha = log(PI))) + geom_raster() + scale_fill_gradient(low = "transparent", high = "green"),
#             ggplot(temp.data, aes(x = Alt, y = Perm, fill = log(PM), alpha = log(PM))) + geom_raster() + scale_fill_gradient(low = "transparent", high = "black"),
#             ggplot(temp.data, aes(x = Alt, y = Perm, fill = log(PR), alpha = log(PR))) + geom_raster() + scale_fill_gradient(low = "transparent", high = "red"),
#             nrow = 1)

# predicciones en la escala de los datos
par(mar = c(2, 2, 1, .5))
layout(matrix(1:6, 2, 3, byrow = TRUE))
# efectos: perm (con 5 valores de altTPI = {min, 25%, median, 75%, max})
loop.spp <- data.frame(Sp = c("PI", "PM", "PR"), Y = c(0.4, 0.08, 0.2), Cols = c("darkgreen", "darkblue", "darkred"))
temp.data <- quantile(LaFad_Alt_TPI(0:1500, 0), c(0, .25, .5, .75, 1))

for(i in loop.spp$Sp) {
  for(j in 1:length(temp.data)) {
    plot(effectfun(Poisson2LaFad2004, "LaFad_Perm_Smooth", marks = i, LaFad_Alt_TPI = temp.data[j]), col = loop.spp$Cols[loop.spp$Sp == i], main = paste("LaFad2004: ", i), ylim = c(0, loop.spp$Y[loop.spp$Sp == i]), lwd = .5 * j, add = j > 1)
    abline(v = c(0.7, -0.86), col = "blue", lwd = 0.5, lty = 2)
  }
}
# efectos: alt (con perm = {min, 25%, median, 75%, max})
temp.data <- quantile(LaFad_Perm_Smooth(0:1500, 0), c(0, .25, .5, .75, 1))

for(i in loop.spp$Sp) {
  for(j in 1:length(temp.data)) {
    plot(effectfun(Poisson2LaFad2004, "LaFad_Alt_TPI", marks = i, LaFad_Perm_Smooth = temp.data[j]), col = loop.spp$Cols[loop.spp$Sp == i], main = paste("LaFad2004: ", i), ylim = c(0, loop.spp$Y[loop.spp$Sp == i]), lwd = .5 * j, add = j > 1)
    abline(v = 0, col = "blue", lwd = 0.5, lty = 2)
  }
}

Figure 22: Partial effects of soil permeability (above) and topographic position (below), according to the selected model, for different values of the other variable (in order of increasing line width: minimum, 25% quantile, median, 75% quantile, and maximum) on the intensity of the point pattern spatial process determining the location of active nests of each Pogonomyrmex ant species in 2004 along a 1.5 km segment of La Fadina dirt road (Ñacuñán Reserve, central Monte desert).

— Doble

In the other dirt road, Doble, the analyses of nested models of increasing complexity show almost the same results, with a few subtleties: [1] the difference among species densities is smaller and [2] the contribution of edaphic variables to explain the variation in Pogonomyrmex spp. density along the road are a bit higher but [3], again, both soil permeability (as a linear predictor in this case) and topography are stronger predictors when allowing for different per-species parameters. Then, [4] both are again important predictors in a combined model, [5] best as independent linear/quadratic effects than as interacting linear effects.

#comparación de modelos Poisson anidados
#(1) comparación modelos homogéneos con intensidad por especie vs. misma intensidad para cada una
cat("[1] Homogeneous models with a single (Model 1) or species-specific (Model 2) intensities\n")
anova(lppm(Doble2004_Pogos ~ 1),
      lppm(Doble2004_Pogos ~ marks),
      test = "LR")

#(2) agregando a la permeabilidad o a la altura como predictores lineales o cuadráticos independientes de la variación espacial de la intensidad  (i.e., afectando a las 3 spp por igual) 
cat("\n[2] Model with species-specific spatially homogeneous intensities (Models 1) or with their intensities varying (in the same way) linearly (Models 2) or quadratically (Models 3) with either of the edaphic variables\n")
anova(lppm(Doble2004_Pogos ~ marks),
      lppm(Doble2004_Pogos ~ marks + Doble_Perm_Smooth),
      lppm(Doble2004_Pogos ~ marks + polynom(Doble_Perm_Smooth, 2)),
      test = "LR")
anova(lppm(Doble2004_Pogos ~ marks),
      lppm(Doble2004_Pogos ~ marks + Doble_Alt_TPI),
      lppm(Doble2004_Pogos ~ marks + polynom(Doble_Alt_TPI, 2)),
      test = "LR")

#(3) agregando a la permeabilidad o a la altura como predictores lineales de la variación espacial de la intensidad  (i.e., afectando a las 3 spp por igual) o con interacción entre edáficas y especies [= test de interacción], de modo lineal o cuadrático
cat("\n[3] Model with heterogeneous intensities varying in the same way linearly with each edaphic factor (Models 1) or with species-specific linear (Models 2) or quadratic (Models 3) relationships\n")
anova(lppm(Doble2004_Pogos ~ marks + Doble_Perm_Smooth),
      lppm(Doble2004_Pogos ~ marks * Doble_Perm_Smooth),
      lppm(Doble2004_Pogos ~ marks * polynom(Doble_Perm_Smooth, 2)),
      test = "LR")
anova(lppm(Doble2004_Pogos ~ marks + Doble_Alt_TPI),      
      lppm(Doble2004_Pogos ~ marks * Doble_Alt_TPI),
      lppm(Doble2004_Pogos ~ marks * polynom(Doble_Alt_TPI, 2)),
      test = "LR")

#(4) agregando a la topografía sobre el modelo que ya tiene a permeabilidad * especie
cat("\n[4] Model with species-specific quadratic effects of soil permeability (Model 1), with an additional independent quadratic effect of topography (Model 2), or with both species-specific effects interacting (Model 3)\n")
anova(lppm(Doble2004_Pogos ~ marks * Doble_Perm_Smooth),
      lppm(Doble2004_Pogos ~ marks * Doble_Perm_Smooth + marks * polynom(Doble_Alt_TPI, 2)),
      lppm(Doble2004_Pogos ~ marks * Doble_Perm_Smooth * polynom(Doble_Alt_TPI, 2)),
      test = "LR")

#(5) lineal+cuadrática independientes, cuadráticas independientes, o lineales interactuando
cat("\n[5] Independent quadratic relationships or interacting linear effects?\n")
AIC(lppm(Doble2004_Pogos ~ marks * Doble_Perm_Smooth * Doble_Alt_TPI), lppm(Doble2004_Pogos ~ marks * Doble_Perm_Smooth + marks * polynom(Doble_Alt_TPI, 2)), lppm(Doble2004_Pogos ~ marks * polynom(Doble_Perm_Smooth, 2) + marks * polynom(Doble_Alt_TPI, 2)))

# single predictor, best lineal model and best quadratic model
PoissonPermDoble2004 <- lppm(Doble2004_Pogos ~ marks * Doble_Perm_Smooth)
PoissonAltDoble2004 <- lppm(Doble2004_Pogos ~ marks * polynom(Doble_Alt_TPI, 2))
PoissonDoble2004 <- lppm(Doble2004_Pogos ~ marks * Doble_Perm_Smooth * Doble_Alt_TPI)
Poisson2Doble2004 <- lppm(Doble2004_Pogos ~ marks * Doble_Perm_Smooth + marks * polynom(Doble_Alt_TPI, 2))

## [1] Homogeneous models with a single (Model 1) or species-specific (Model 2) intensities
## Analysis of Deviance Table
## 
## Model 1: ~1   Poisson
## Model 2: ~marks   Poisson
##   Npar Df Deviance Pr(>Chi)  
## 1    1                       
## 2    3  2   5.9196  0.05183 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## [2] Model with species-specific spatially homogeneous intensities (Models 1) or with their intensities varying (in the same way) linearly (Models 2) or quadratically (Models 3) with either of the edaphic variables
## Analysis of Deviance Table
## 
## Model 1: ~marks   Poisson
## Model 2: ~marks + Doble_Perm_Smooth   Poisson
## Model 3: ~marks + (Doble_Perm_Smooth + I(Doble_Perm_Smooth^2))    Poisson
##   Npar Df Deviance Pr(>Chi)  
## 1    3                       
## 2    4  1   6.2468  0.01244 *
## 3    5  1   0.6882  0.40679  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Analysis of Deviance Table
## 
## Model 1: ~marks   Poisson
## Model 2: ~marks + Doble_Alt_TPI   Poisson
## Model 3: ~marks + (Doble_Alt_TPI + I(Doble_Alt_TPI^2))    Poisson
##   Npar Df Deviance Pr(>Chi)  
## 1    3                       
## 2    4  1   4.5804  0.03234 *
## 3    5  1   3.6888  0.05478 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## [3] Model with heterogeneous intensities varying in the same way linearly with each edaphic factor (Models 1) or with species-specific linear (Models 2) or quadratic (Models 3) relationships
## Analysis of Deviance Table
## 
## Model 1: ~marks + Doble_Perm_Smooth   Poisson
## Model 2: ~marks * Doble_Perm_Smooth   Poisson
## Model 3: ~marks * (Doble_Perm_Smooth + I(Doble_Perm_Smooth^2))    Poisson
##   Npar Df Deviance        Pr(>Chi)    
## 1    4                                
## 2    6  2   41.818 0.0000000008306 ***
## 3    9  3    5.102          0.1645    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Analysis of Deviance Table
## 
## Model 1: ~marks + Doble_Alt_TPI   Poisson
## Model 2: ~marks * Doble_Alt_TPI   Poisson
## Model 3: ~marks * (Doble_Alt_TPI + I(Doble_Alt_TPI^2))    Poisson
##   Npar Df Deviance  Pr(>Chi)    
## 1    4                          
## 2    6  2   21.829 0.0000182 ***
## 3    9  3   18.232  0.000394 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## [4] Model with species-specific quadratic effects of soil permeability (Model 1), with an additional independent quadratic effect of topography (Model 2), or with both species-specific effects interacting (Model 3)
## Analysis of Deviance Table
## 
## Model 1: ~marks * Doble_Perm_Smooth   Poisson
## Model 2: ~marks * Doble_Perm_Smooth + marks * (Doble_Alt_TPI + I(Doble_Alt_TPI^2))    Poisson
## Model 3: ~marks * Doble_Perm_Smooth * (Doble_Alt_TPI + I(Doble_Alt_TPI^2))    Poisson
##   Npar Df Deviance  Pr(>Chi)    
## 1    6                          
## 2   12  6  23.8518 0.0005561 ***
## 3   18  6   9.2635 0.1592900    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## [5] Independent quadratic relationships or interacting linear effects?
##                                                                                                   df
## lppm(Doble2004_Pogos ~ marks * Doble_Perm_Smooth * Doble_Alt_TPI)                                 12
## lppm(Doble2004_Pogos ~ marks * Doble_Perm_Smooth + marks * polynom(Doble_Alt_TPI, 2))             12
## lppm(Doble2004_Pogos ~ marks * polynom(Doble_Perm_Smooth, 2) + marks * polynom(Doble_Alt_TPI, 2)) 15
##                                                                                                        AIC
## lppm(Doble2004_Pogos ~ marks * Doble_Perm_Smooth * Doble_Alt_TPI)                                 664.2115
## lppm(Doble2004_Pogos ~ marks * Doble_Perm_Smooth + marks * polynom(Doble_Alt_TPI, 2))             660.6903
## lppm(Doble2004_Pogos ~ marks * polynom(Doble_Perm_Smooth, 2) + marks * polynom(Doble_Alt_TPI, 2)) 665.0852

The bivariate models, either with quadratic or interacting linear relationships. seems to once again better predict some sectors with the presence of colonies than the univariate models (Fig. 23). The combination of edaphic values associated with the lowest densities of PI and PM are very similar to the previous road, though with some differences on the association between topography and the highest predicted densities of PI (Fig. 24). In this road the agreement between high predicted intensity and presence of active PR colonies is better, and relies heavier on the mid values of relative altitude instead of on soil permeability (Fig. 25).

# predicciones vs colonias
par(mar = c(2, 2, 1, .5))
layout(matrix(1:3, 3, 1))

loop.spp <- data.frame(Sp = c("PI", "PM", "PR"), Points = c(0.09, 0.04, 0.03), Cols = c("darkgreen", "darkblue", "darkred"))
graph.count <- 0

for(i in loop.spp$Sp) {
  plot(x = seq(0, 1500, 5), y = predict(Poisson2Doble2004, locations = data.frame(x = seq(0, 1500, 5), y = 30, marks = i)), type = "l", lwd = 2, main = paste("Doble2004: ", i))
  lines(x = seq(0, 1500, 5), y = predict(PoissonDoble2004, locations = data.frame(x = seq(0, 1500, 5), y = 30, marks = i)), col = "darkgrey")
  lines(x = seq(0, 1500, 5), y = predict(PoissonPermDoble2004, locations = data.frame(x = seq(0, 1500, 5), y = 30, marks = i)), col = "lightblue3")
  lines(x = seq(0, 1500, 5), y = predict(PoissonAltDoble2004, locations = data.frame(x = seq(0, 1500, 5), y = 30, marks = i)), col = "sienna")
  points(x = Doble2004_Pogos$data$x[marks(Doble2004_Pogos) == i], y = rep(loop.spp$Points[loop.spp$Sp == i], length(Doble2004_Pogos$data$x[marks(Doble2004_Pogos) == i])), col = loop.spp$Cols[loop.spp$Sp == i])
  abline(h = intensity(Doble2004_Pogos)[i], lty = 2, col = "red")
}

Figure 23: Same as Fig. 20 but for the active Pogonomyrmex colonies in 2004 along the 1.5km segment of the Doble dirt road.

# predicción bivariada en la escala del modelo; para matriz de n×n (=50) valores del rango observado de cada variable en la transecta
temp.data <- data.frame(Alt = rep(seq(min(Doble_Alt_TPI(0:1500, 0)), max(Doble_Alt_TPI(0:1500, 0)), length.out = 50), 50), Perm = rep(seq(min(Doble_Perm_Smooth(0:1500, 0)), max(Doble_Perm_Smooth(0:1500, 0)), length.out = 50), each = 50)) %>%
  mutate(PI = predict(Poisson2Doble2004, locations = data.frame(marks = "PI", x = rep(0, length(Alt)), y = rep(0, length(Alt))), covariates = data.frame(Doble_Alt_TPI = Alt, Doble_Perm_Smooth = Perm))) %>%
  mutate(PM = predict(Poisson2Doble2004, locations = data.frame(marks = "PM", x = rep(0, length(Alt)), y = rep(0, length(Alt))), covariates = data.frame(Doble_Alt_TPI = Alt, Doble_Perm_Smooth = Perm))) %>%
  mutate(PR = predict(Poisson2Doble2004, locations = data.frame(marks = "PR", x = rep(0, length(Alt)), y = rep(0, length(Alt))), covariates = data.frame(Doble_Alt_TPI = Alt, Doble_Perm_Smooth = Perm))) %>% pivot_longer(., cols = c("PI", "PM", "PR"), names_to = "Sp", values_to = "Lambda") # para un solo ggplot

ggplot(temp.data, aes(x = Alt, y = Perm, fill = log(Lambda))) +
  geom_raster() +
  facet_wrap(.~Sp) +
  scale_fill_gradient(name = "log(λ)", high = "black", low = "#FFFFFF00") +
  labs(x = "TPI [m]", y = "Soil permeability (SPI)") +
  geom_hline(yintercept = c(0.7, -0.86), colour = "blue", linewidth = 0.5, linetype = 2) +
  geom_vline(xintercept = 0, colour = "blue", linewidth = 0.5, linetype = 2) +
  tema +
  theme(legend.position = "right")

Figure 24: Same as Fig. 21 but for the active Pogonomyrmex colonies in 2004 along the 1.5km segment of the Doble dirt road.

#para hacer gráficos separados eliminar pivot_longer y llamar a cada variable en un gráfico
#grid.arrange(ggplot(temp.data, aes(x = Alt, y = Perm, fill = log(PI), alpha = log(PI))) + geom_raster() + scale_fill_gradient(low = "transparent", high = "green"),
#             ggplot(temp.data, aes(x = Alt, y = Perm, fill = log(PM), alpha = log(PM))) + geom_raster() + scale_fill_gradient(low = "transparent", high = "black"),
#             ggplot(temp.data, aes(x = Alt, y = Perm, fill = log(PR), alpha = log(PR))) + geom_raster() + scale_fill_gradient(low = "transparent", high = "red"),
#             nrow = 1)

# predicciones en la escala de los datos
par(mar = c(2, 2, 1, .5))
layout(matrix(1:6, 2, 3, byrow = TRUE))

# efectos: perm (con 5 valores de altTPI = {min, 25%, median, 75%, max})
loop.spp <- data.frame(Sp = c("PI", "PM", "PR"), Y = c(0.3, 0.08, 0.05), Cols = c("darkgreen", "darkblue", "darkred"))
temp.data <- quantile(Doble_Alt_TPI(0:1500, 0), c(0, .25, .5, .75, 1))

for(i in loop.spp$Sp) {
  for(j in 1:length(temp.data)) {
    plot(effectfun(Poisson2Doble2004, "Doble_Perm_Smooth", marks = i, Doble_Alt_TPI = temp.data[j]), col = loop.spp$Cols[loop.spp$Sp == i], main = paste("Doble2004: ", i), ylim = c(0, loop.spp$Y[loop.spp$Sp == i]), lwd = .5 * j, add = j > 1)
    abline(v = c(0.7, -0.86), col = "blue", lwd = 0.5, lty = 2)
  }
}

# efectos: alt (con perm = {min, 25%, median, 75%, max})
temp.data <- quantile(Doble_Perm_Smooth(0:1500, 0), c(0, .25, .5, .75, 1))

for(i in loop.spp$Sp) {
  for(j in 1:length(temp.data)) {
    plot(effectfun(Poisson2Doble2004, "Doble_Alt_TPI", marks = i, Doble_Perm_Smooth = temp.data[j]), col = loop.spp$Cols[loop.spp$Sp == i], main = paste("Doble2004: ", i), ylim = c(0, loop.spp$Y[loop.spp$Sp == i]), lwd = .5 * j, add = j > 1)
    abline(v = 0, col = "blue", lwd = 0.5, lty = 2)
  }
}

Figure 25: Same as Fig. 22 but for the active Pogonomyrmex colonies in 2004 along the 1.5km segment of the Doble dirt road.

— Roads as replicates

Given the highly similar point pattern models fitted to the information of each of the two road segments, they can be analysed as replicates assuming the same spatial process(es). Again, the best predictive model involves each of the three species associated in very different way to both edaphic factors, slightly better in quadratic, non-interacting form than as interacting linear effects:

# ¡OJO! (Baddeley dixit 20220830):
# We have not yet implemented mppm for linear networks. However there is a way to 'cheat' and obtain some of the analysis at least. Assume you have a hyperframe H that contains a column 'X' of 'lpp' objects. First convert these to quadrature schemes via H$X <- with(H, quadscheme(X))
# Then call 'mppm' on this modified hyperframe.
# You can read off the fitted coefficients, perform anova, etc.
# However you can't simulate from the fitted model,
# and I don't guarantee that other things will be correct even if they appear to work.
# Please do not share this information widely because it is a temporary hack that may cease to work later.

# So…
Ambas2004 <- hyperframe(Nidos = solist(quadscheme(LaFad2004_Pogos), quadscheme(Doble2004_Pogos)),
                        Perm = solist(LaFad_Perm_Smooth, Doble_Perm_Smooth),
                        Alt = solist(LaFad_Alt_TPI, Doble_Alt_TPI), 
                        Picada = as.factor(c("LF", "D")))

#(1) agregando a la permeabilidad o a la altura como predictores lineales de la variación espacial de la intensidad  (i.e., afectando a las 3 spp por igual) o con interacción entre edáficas y especies [= test de interacción], de modo lineal o cuadrático
cat("\n[1] Model with heterogeneous intensities varying in the same way linearly with each edaphic factor (Models 1) or with species-specific linear (Models 2) or quadratic (Models 3) relationships\n")

## 
## [1] Model with heterogeneous intensities varying in the same way linearly with each edaphic factor (Models 1) or with species-specific linear (Models 2) or quadratic (Models 3) relationships

anova(mppm(data = Ambas2004, formula = Nidos ~ marks + Perm),
      mppm(data = Ambas2004, formula = Nidos ~ marks * Perm),
      mppm(data = Ambas2004, formula = Nidos ~ marks * polynom(Perm, 2)),
      test = "LR")

## Analysis of Deviance Table
## 
## Model 1: Nidos ~ marks + Perm     Poisson
## Model 2: Nidos ~ marks * Perm     Poisson
## Model 3: Nidos ~ marks * polynom(Perm, 2)     Poisson
##   Npar Df Deviance  Pr(>Chi)    
## 1    4                          
## 2    6  2   60.934 5.867e-14 ***
## 3    9  3   14.845  0.001954 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(mppm(data = Ambas2004, formula = Nidos ~ marks + Alt),
      mppm(data = Ambas2004, formula = Nidos ~ marks * Alt),
      mppm(data = Ambas2004, formula = Nidos ~ marks * polynom(Alt, 2)),
      test = "LR")

## Analysis of Deviance Table
## 
## Model 1: Nidos ~ marks + Alt      Poisson
## Model 2: Nidos ~ marks * Alt      Poisson
## Model 3: Nidos ~ marks * polynom(Alt, 2)      Poisson
##   Npar Df Deviance      Pr(>Chi)    
## 1    4                              
## 2    6  2   35.079 0.00000002414 ***
## 3    9  3   24.008 0.00002488495 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#(2) agregando a la topografía sobre el modelo que ya tiene a permeabilidad * especie
cat("\n[2] Model with species-specific quadratic effects of soil permeability (Model 1), with an additional independent quadratic effect of topography (Model 2), or with both species-specific effects interacting (Model 3)\n")

## 
## [2] Model with species-specific quadratic effects of soil permeability (Model 1), with an additional independent quadratic effect of topography (Model 2), or with both species-specific effects interacting (Model 3)

anova(mppm(data = Ambas2004, formula = Nidos ~ marks * polynom(Perm, 2)),
      mppm(data = Ambas2004, formula = Nidos ~ marks * polynom(Perm, 2) + marks * polynom(Alt, 2)),
      mppm(data = Ambas2004, formula = Nidos ~ marks * polynom(Perm, 2) * polynom(Alt, 2)),
      test = "LR")

## Analysis of Deviance Table
## 
## Model 1: Nidos ~ marks * polynom(Perm, 2)     Poisson
## Model 2: Nidos ~ marks * polynom(Perm, 2) + marks * polynom(Alt, 2)   Poisson
## Model 3: Nidos ~ marks * polynom(Perm, 2) * polynom(Alt, 2)   Poisson
##   Npar Df Deviance   Pr(>Chi)    
## 1    9                           
## 2   15  6   31.407 0.00002119 ***
## 3   27 12   18.805    0.09334 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#(3) cuadráticas independientes o lineales interactuando
cat("\n[3] Independent quadratic relationships or interacting linear effects?\n")

## 
## [3] Independent quadratic relationships or interacting linear effects?

AIC(mppm(data = Ambas2004, formula = Nidos ~ marks * Perm * Alt))

## [1] 2415.724

AIC(mppm(data = Ambas2004, formula = Nidos ~ marks * polynom(Perm, 2) + marks * polynom(Alt, 2)))

## [1] 2406.976

PoissonAmbas2004 <- mppm(data = Ambas2004, formula = Nidos ~ marks * Perm * Alt)
Poisson2Ambas2004 <- mppm(data = Ambas2004, formula = Nidos ~ marks * polynom(Perm, 2) + marks * polynom(Alt, 2))

The inclusion of a term for road identity as explanatory variable either modifying the density of all or of each species (an intercept-model in mixed models jargon) is not justified for any of the two candidate models, supporting the assumption of common spatial processes. Assuming the edaphic effects are the same for both replicates then the (minor) differences observed in the density of Pogonomyrmex spp. active colonies in each road result from the (minor) differences in their characteristics of soil permeability and topography:

#(1) comparación modelos con efecto picada (i.e., simil efecto bloque; proceso independiente de segmento o no)
cat("[1] Best model assuming roads as replicates (Model 1) vs., adding road identity as explanatory variable of road density with a single (Model 2) or species-specific parameters (Model 3)\n")
anova(PoissonAmbas2004, update(PoissonAmbas2004, . ~ . + Picada), update(PoissonAmbas2004, . ~ . + Picada * marks), test = "LRT")
anova(Poisson2Ambas2004, update(Poisson2Ambas2004, . ~ . + Picada), update(Poisson2Ambas2004, . ~ . + Picada * marks), test = "LRT")

## [1] Best model assuming roads as replicates (Model 1) vs., adding road identity as explanatory variable of road density with a single (Model 2) or species-specific parameters (Model 3)
## Analysis of Deviance Table
## 
## Model 1: Nidos ~ marks * Perm * Alt   Poisson
## Model 2: Nidos ~ marks + Perm + Alt + Picada + marks:Perm + marks:Alt + Perm:Alt + marks:Perm:Alt     Poisson
## Model 3: Nidos ~ marks + Perm + Alt + Picada + marks:Perm + marks:Alt + Perm:Alt + marks:Picada + marks:Perm:Alt      Poisson
##   Npar Df Deviance Pr(>Chi)
## 1   12                     
## 2   13  1  1.27268   0.2593
## 3   15  2  0.86093   0.6502
## Analysis of Deviance Table
## 
## Model 1: Nidos ~ marks * polynom(Perm, 2) + marks * polynom(Alt, 2)   Poisson
## Model 2: Nidos ~ marks + polynom(Perm, 2) + polynom(Alt, 2) + Picada + marks:polynom(Perm, 2) + marks:polynom(Alt, 2)     Poisson
## Model 3: Nidos ~ marks + polynom(Perm, 2) + polynom(Alt, 2) + Picada + marks:polynom(Perm, 2) + marks:polynom(Alt, 2) + marks:Picada      Poisson
##   Npar Df Deviance Pr(>Chi)
## 1   15                     
## 2   16  1   1.2182   0.2697
## 3   18  2   1.8617   0.3942

Predictions from this best model (Figs. 26–29) agree with those already mentioned for each of the roads separately, meaning:

• the density of PI decreases strongly with soil permeability except when in the lowest sectors, and decreases slightly with relative height except when in the most impermeable soils
• the density of PM increases with soil permeability (with a stronger change when in relative depressions) and increases with relative altitude
• the density of PR increases in sectors where relative altitude and soil permeability compensate: it decreases with relative altitude in permeable soils but increases with height when in the most impermeable soils (or, from the other point of view, increases with permeability in depressions but with impermeability in relative peaks).

par(mar = c(2, 2, 1, .5))
layout(matrix(1:3, 3, 1))
# predict y effectfun no funcionan con mppm sobre objetos lpp ('cause hack…)
# con subfits(mppm)[[1]] obtengo ppm (y no lppm) ¿y entonces predict estima mal los lamdba (los relativos están ok, ¿pero las magnitudes son areales en vez de lineales?
# So… "a mano":
# predicción bivariada en la escala del modelo multiplicando coeficientes a partir de coef(mppm)
# matriz de n×n (=50) valores dentro del rango observado de cada variable en ambas transectas
temp.data <- data.frame(Alt = rep(seq(min(c(LaFad_Alt_TPI(0:1500, 0), Doble_Alt_TPI(0:1500, 0))), max(c(LaFad_Alt_TPI(0:1500, 0), Doble_Alt_TPI(0:1500, 0))), length.out = 50), 50), 
                        Perm = rep(seq(min(c(LaFad_Perm_Smooth(0:1500, 0), Doble_Perm_Smooth(0:1500, 0))), max(c(LaFad_Perm_Smooth(0:1500, 0), Doble_Perm_Smooth(0:1500, 0))), length.out = 50), each = 50)) %>%
    bind_rows("PI" = ., "PM" = ., "PR" = ., .id = "Sp")
# "matrices de diseño" construidas a partir de coef(mppm)
# para modelo de lineales interactuando
coef.data <- coef(PoissonAmbas2004)
temp.matrix <- temp.data %>%
  mutate(V01 = TRUE, V02 = Sp == "PM", V03 = Sp == "PR", V04 = Perm, V05 = Alt, V06 = V02 * V04, V07 = V03 * V04, V08 = V02 * V05, V09 = V03 * V05, V10 = V04 * V05, V11 = V02 * V10, V12 = V03 * V10, .keep = "none")  
temp.data$Lambda1 <- exp(as.matrix(temp.matrix) %*% coef.data)[,1]
# para modelo de cuadráticas independientes
coef.data <- coef(Poisson2Ambas2004)
temp.matrix <- temp.data %>%
  mutate(V01 = TRUE, V02 = Sp == "PM", V03 = Sp == "PR", V04 = Perm, V05 = Perm ^ 2, V06 = Alt, V07 = Alt ^ 2, V08 = V02 * V04, V09 = V03 * V04, V10 = V02 * V05, V11 = V03 * V05, V12 = V02 * V06, V13 = V03 * V06, V14 = V02 * V07, V15 = V03 * V07, .keep = "none")
temp.data$Lambda2 <- exp(as.matrix(temp.matrix) %*% coef.data)[,1]

grid.arrange(ggplot(temp.data, aes(x = Alt, y = Perm, fill = log(Lambda1))) +
                geom_raster() +
                facet_wrap(. ~ Sp) +
                scale_fill_gradient(name = "log(λ)", high = "black", low = "#FFFFFF00") +
                labs(x = NULL, y = NULL, caption = "Linear interacting") +
                geom_hline(yintercept = c(0.7, -0.86), colour = "blue", linewidth = 0.5, linetype = 2) +
                geom_vline(xintercept = 0, colour = "blue", linewidth = 0.5, linetype = 2) +
                tema +
                theme(legend.position = "right"),
              ggplot(temp.data, aes(x = Alt, y = Perm, fill = log(Lambda2))) +
                geom_raster() +
                facet_wrap(. ~ Sp) +
                scale_fill_gradient(name = "log(λ)", high = "black", low = "#FFFFFF00") +
                labs(x = NULL, y = NULL, caption = "Quadratic independent") +
                geom_hline(yintercept = c(0.7, -0.86), colour = "blue", linewidth = 0.5, linetype = 2) +
                geom_vline(xintercept = 0, colour = "blue", linewidth = 0.5, linetype = 2) +
                tema +
                theme(legend.position = "right"),
              nrow = 2, bottom = "TPI [m]", left = "Soil permeability (SPI)")

Figure 26: Same as Fig. 21 but for the active Pogonomyrmex colonies in 2004 along both 1.5 km segments of dirt roads. Abpve: Model with interacting linear effects. Below: Model with independent quadratic effects.

#para hacer gráficos separados (diferentes escalas y colores) eliminar pivot_longer y llamar a cada variable en un ggplot diferente

grid.arrange(ggplot(subset(temp.data, Sp == "PI"), aes(x = Alt, y = Perm, fill = log(Lambda1))) +
               geom_raster() +
               scale_fill_gradient(low = "transparent", high = "darkgreen") +
               labs(x = NULL, y = NULL, subtitle = "PI", caption = "") +
               geom_hline(yintercept = c(0.7, -0.86), colour = "blue", linewidth = 0.5, linetype = 2) +
               geom_vline(xintercept = 0, colour = "blue", linewidth = 0.5, linetype = 2) +
               tema,
             ggplot(subset(temp.data, Sp == "PM"), aes(x = Alt, y = Perm, fill = log(Lambda1))) +
               geom_raster() +
               scale_fill_gradient(low = "transparent", high = "darkblue") +
               labs(x = NULL, y = NULL, subtitle = "PM", caption = "") +
               geom_hline(yintercept = c(0.7, -0.86), colour = "blue", linewidth = 0.5, linetype = 2) +
               geom_vline(xintercept = 0, colour = "blue", linewidth = 0.5, linetype = 2) +
               tema,
             ggplot(subset(temp.data, Sp == "PR"), aes(x = Alt, y = Perm, fill = log(Lambda1))) +
               geom_raster() +
               scale_fill_gradient(low = "transparent", high = "darkred") +
               labs(x = NULL, y = NULL, subtitle = "PR", caption = "Linear interacting") +
               geom_hline(yintercept = c(0.7, -0.86), colour = "blue", linewidth = 0.5, linetype = 2) +
               geom_vline(xintercept = 0, colour = "blue", linewidth = 0.5, linetype = 2) +
               tema,
             ggplot(subset(temp.data, Sp == "PI"), aes(x = Alt, y = Perm, fill = log(Lambda2))) +
               geom_raster() +
               scale_fill_gradient(low = "transparent", high = "darkgreen") +
               labs(x = NULL, y = NULL, subtitle = "", caption = "") +
               geom_hline(yintercept = c(0.7, -0.86), colour = "blue", linewidth = 0.5, linetype = 2) +
               geom_vline(xintercept = 0, colour = "blue", linewidth = 0.5, linetype = 2) +
               tema,
             ggplot(subset(temp.data, Sp == "PM"), aes(x = Alt, y = Perm, fill = log(Lambda2))) +
               geom_raster() +
               scale_fill_gradient(low = "transparent", high = "darkblue") +
               labs(x = NULL, y = NULL, subtitle = "", caption = "") +
               geom_hline(yintercept = c(0.7, -0.86), colour = "blue", linewidth = 0.5, linetype = 2) +
               geom_vline(xintercept = 0, colour = "blue", linewidth = 0.5, linetype = 2) +
               tema,
             ggplot(subset(temp.data, Sp == "PR"), aes(x = Alt, y = Perm, fill = log(Lambda2))) +
               geom_raster() +
               scale_fill_gradient(low = "transparent", high = "darkred") +
               labs(x = NULL, y = NULL, subtitle = "", caption = "Quadratic independent") +
               geom_hline(yintercept = c(0.7, -0.86), colour = "blue", linewidth = 0.5, linetype = 2) +
               geom_vline(xintercept = 0, colour = "blue", linewidth = 0.5, linetype = 2) +
               tema,
             nrow = 2, bottom = "TPI [m]", left = "Soil permeability (SPI)")

Figure 27: Same as Fig. 26 but with colour transparency defined for each species range of intensity values (= different scales).

# predicciones en la escala de los datos (simil effectfun)
# efectos: perm (con alt = {min, 25%, median, 75%, max})
temp.data <- data.frame(Alt = rep(fivenum(c(LaFad_Alt_TPI(0:1500, 0), Doble_Alt_TPI(0:1500, 0))), each = 50), 
                        Perm = rep(seq(min(c(LaFad_Perm_Smooth(0:1500, 0), Doble_Perm_Smooth(0:1500, 0))), max(c(LaFad_Perm_Smooth(0:1500, 0), Doble_Perm_Smooth(0:1500, 0))), length.out = 50), 5),
                        Lwidth = rep(1:5, each = 50)) %>%
                        bind_rows("PI" = ., "PM" = ., "PR" = ., .id = "Sp")

coef.data <- coef(PoissonAmbas2004)
temp.matrix <- temp.data %>%
  mutate(V01 = TRUE, V02 = Sp == "PM", V03 = Sp == "PR", V04 = Perm, V05 = Alt, V06 = V02 * V04, V07 = V03 * V04, V08 = V02 * V05, V09 = V03 * V05, V10 = V04 * V05, V11 = V02 * V10, V12 = V03 * V10, .keep = "none")  
temp.data$Lambda <- exp(as.matrix(temp.matrix) %*% coef.data)[,1]

grid.arrange(ggplot(subset(temp.data, Sp == "PI"), aes(x = Perm, y = Lambda, group = as.factor(Alt), linewidth = Lwidth)) +
               geom_line(colour = "darkgreen", lineend = "round") +
               scale_linewidth(breaks = 1.5:5.5, range = c(0.3, 1.3)) +
               geom_vline(xintercept = c(0.7, -0.86), colour = "blue", linewidth = 0.5, linetype = 2) +
               coord_cartesian(ylim = c(0, 0.007)) +
               labs(subtitle = "PI", y  = NULL, x = NULL) +
               tema,
             ggplot(subset(temp.data, Sp == "PM"), aes(x = Perm, y = Lambda, group = as.factor(Alt), linewidth = Lwidth)) +
               geom_line(colour = "darkblue", lineend = "round") +
               scale_linewidth(breaks = 1.5:5.5, range = c(0.3, 1.3)) +
               coord_cartesian(ylim = c(0, 0.0015)) +
               geom_vline(xintercept = c(0.7, -0.86), colour = "blue", linewidth = 0.5, linetype = 2) +
               labs(subtitle = "PM", y  = NULL, x = NULL) +
               tema,
             ggplot(subset(temp.data, Sp == "PR"), aes(x = Perm, y = Lambda, group = as.factor(Alt), linewidth = Lwidth)) +
               geom_line(colour = "darkred", lineend = "round") +
               scale_linewidth(breaks = 1.5:5.5, range = c(0.3, 1.3)) +
               coord_cartesian(ylim = c(0, 0.0035)) +
               geom_vline(xintercept = c(0.7, -0.86), colour = "blue", linewidth = 0.5, linetype = 2) +
               labs(subtitle = "PR", y  = NULL, x = NULL) +
               tema,
             nrow = 1, left = "Lambda", bottom = "Soil permeability (SPI)")

# efectos: alt (con perm = {min, 25%, median, 75%, max})
temp.data <- data.frame(Alt = rep(seq(min(c(LaFad_Alt_TPI(0:1500, 0), Doble_Alt_TPI(0:1500, 0))), max(c(LaFad_Alt_TPI(0:1500, 0), Doble_Alt_TPI(0:1500, 0))), length.out = 50), 5), 
                        Perm = rep(fivenum(c(LaFad_Perm_Smooth(0:1500, 0), Doble_Perm_Smooth(0:1500, 0))), each = 50),
                        Lwidth = rep(1:5, each = 50)) %>%
                        bind_rows("PI" = ., "PM" = ., "PR" = ., .id = "Sp")

temp.matrix <- temp.data %>%
  mutate(V01 = TRUE, V02 = Sp == "PM", V03 = Sp == "PR", V04 = Perm, V05 = Alt, V06 = V02 * V04, V07 = V03 * V04, V08 = V02 * V05, V09 = V03 * V05, V10 = V04 * V05, V11 = V02 * V10, V12 = V03 * V10, .keep = "none")  
temp.data$Lambda <- exp(as.matrix(temp.matrix) %*% coef.data)[,1]

grid.arrange(ggplot(subset(temp.data, Sp == "PI"), aes(x = Alt, y = Lambda, group = as.factor(Perm), linewidth = Lwidth)) +
               geom_line(colour = "darkgreen", lineend = "round") +
               scale_linewidth(breaks = 1.5:5.5, range = c(0.3, 1.3)) +
               coord_cartesian(ylim = c(0, 0.007)) +
               geom_vline(xintercept = 0, colour = "blue", linewidth = 0.5, linetype = 2) +
               labs(subtitle = "PI", y  = NULL, x = NULL) +
               tema,
             ggplot(subset(temp.data, Sp == "PM"), aes(x = Alt, y = Lambda, group = as.factor(Perm), linewidth = Lwidth)) +
               geom_line(colour = "darkblue", lineend = "round") +
               scale_linewidth(breaks = 1.5:5.5, range = c(0.3, 1.3)) +
               coord_cartesian(ylim = c(0, 0.0015)) +
               geom_vline(xintercept = 0, colour = "blue", linewidth = 0.5, linetype = 2) +
               labs(subtitle = "PM", y  = NULL, x = NULL) +
               tema,
             ggplot(subset(temp.data, Sp == "PR"), aes(x = Alt, y = Lambda, group = as.factor(Perm), linewidth = Lwidth)) +
               geom_line(colour = "darkred", lineend = "round") +
               scale_linewidth(breaks = 1.5:5.5, range = c(0.3, 1.3)) +
               coord_cartesian(ylim = c(0, 0.0035)) +
               geom_vline(xintercept = 0, colour = "blue", linewidth = 0.5, linetype = 2) +
               labs(subtitle = "PR", y  = NULL, x = NULL) +
               tema,
             nrow = 1, left = "Lambda", bottom = "TPI [m]")

Figure 28: Same as Fig. 22 but for the active Pogonomyrmex colonies in 2004 along both 1.5 km segments of dirt roads.

# predicciones en la escala de los datos (simil effectfun)
# efectos: perm (con alt = {min, 25%, median, 75%, max})
temp.data <- data.frame(Alt = rep(fivenum(c(LaFad_Alt_TPI(0:1500, 0), Doble_Alt_TPI(0:1500, 0))), each = 50), 
                        Perm = rep(seq(min(c(LaFad_Perm_Smooth(0:1500, 0), Doble_Perm_Smooth(0:1500, 0))), max(c(LaFad_Perm_Smooth(0:1500, 0), Doble_Perm_Smooth(0:1500, 0))), length.out = 50), 5),
                        Lwidth = rep(1:5, each = 50)) %>%
                        bind_rows("PI" = ., "PM" = ., "PR" = ., .id = "Sp")

coef.data <- coef(Poisson2Ambas2004)
temp.matrix <- temp.data %>%
  mutate(V01 = TRUE, V02 = Sp == "PM", V03 = Sp == "PR", V04 = Perm, V05 = Perm ^ 2, V06 = Alt, V07 = Alt ^ 2, V08 = V02 * V04, V09 = V03 * V04, V10 = V02 * V05, V11 = V03 * V05, V12 = V02 * V06, V13 = V03 * V06, V14 = V02 * V07, V15 = V03 * V07, .keep = "none")
temp.data$Lambda <- exp(as.matrix(temp.matrix) %*% coef.data)[,1]

grid.arrange(ggplot(subset(temp.data, Sp == "PI"), aes(x = Perm, y = Lambda, group = as.factor(Alt), linewidth = Lwidth)) +
               geom_line(colour = "darkgreen", lineend = "round") +
               scale_linewidth(breaks = 1.5:5.5, range = c(0.3, 1.3)) +
               geom_vline(xintercept = c(0.7, -0.86), colour = "blue", linewidth = 0.5, linetype = 2) +
               coord_cartesian(ylim = c(0, 0.0045)) +
               labs(subtitle = "PI", y  = NULL, x = NULL) +
               tema,
             ggplot(subset(temp.data, Sp == "PM"), aes(x = Perm, y = Lambda, group = as.factor(Alt), linewidth = Lwidth)) +
               geom_line(colour = "darkblue", lineend = "round") +
               scale_linewidth(breaks = 1.5:5.5, range = c(0.3, 1.3)) +
               coord_cartesian(ylim = c(0, 0.0012)) +
               geom_vline(xintercept = c(0.7, -0.86), colour = "blue", linewidth = 0.5, linetype = 2) +
               labs(subtitle = "PM", y  = NULL, x = NULL) +
               tema,
             ggplot(subset(temp.data, Sp == "PR"), aes(x = Perm, y = Lambda, group = as.factor(Alt), linewidth = Lwidth)) +
               geom_line(colour = "darkred", lineend = "round") +
               scale_linewidth(breaks = 1.5:5.5, range = c(0.3, 1.3)) +
               coord_cartesian(ylim = c(0, 0.0012)) +
               geom_vline(xintercept = c(0.7, -0.86), colour = "blue", linewidth = 0.5, linetype = 2) +
               labs(subtitle = "PR", y  = NULL, x = NULL) +
               tema,
             nrow = 1, left = "Lambda", bottom = "Soil permeability (SPI)")

# efectos: alt (con perm = {min, 25%, median, 75%, max})
temp.data <- data.frame(Alt = rep(seq(min(c(LaFad_Alt_TPI(0:1500, 0), Doble_Alt_TPI(0:1500, 0))), max(c(LaFad_Alt_TPI(0:1500, 0), Doble_Alt_TPI(0:1500, 0))), length.out = 50), 5), 
                        Perm = rep(fivenum(c(LaFad_Perm_Smooth(0:1500, 0), Doble_Perm_Smooth(0:1500, 0))), each = 50),
                        Lwidth = rep(1:5, each = 50)) %>%
                        bind_rows("PI" = ., "PM" = ., "PR" = ., .id = "Sp")

temp.matrix <- temp.data %>%
  mutate(V01 = TRUE, V02 = Sp == "PM", V03 = Sp == "PR", V04 = Perm, V05 = Perm ^ 2, V06 = Alt, V07 = Alt ^ 2, V08 = V02 * V04, V09 = V03 * V04, V10 = V02 * V05, V11 = V03 * V05, V12 = V02 * V06, V13 = V03 * V06, V14 = V02 * V07, V15 = V03 * V07, .keep = "none")
temp.data$Lambda <- exp(as.matrix(temp.matrix) %*% coef.data)[,1]

grid.arrange(ggplot(subset(temp.data, Sp == "PI"), aes(x = Alt, y = Lambda, group = as.factor(Perm), linewidth = Lwidth)) +
               geom_line(colour = "darkgreen", lineend = "round") +
               scale_linewidth(breaks = 1.5:5.5, range = c(0.3, 1.3)) +
               coord_cartesian(ylim = c(0, 0.0045)) +
               geom_vline(xintercept = 0, colour = "blue", linewidth = 0.5, linetype = 2) +
               labs(subtitle = "PI", y  = NULL, x = NULL) +
               tema,
             ggplot(subset(temp.data, Sp == "PM"), aes(x = Alt, y = Lambda, group = as.factor(Perm), linewidth = Lwidth)) +
               geom_line(colour = "darkblue", lineend = "round") +
               scale_linewidth(breaks = 1.5:5.5, range = c(0.3, 1.3)) +
               coord_cartesian(ylim = c(0, 0.0012)) +
               geom_vline(xintercept = 0, colour = "blue", linewidth = 0.5, linetype = 2) +
               labs(subtitle = "PM", y  = NULL, x = NULL) +
               tema,
             ggplot(subset(temp.data, Sp == "PR"), aes(x = Alt, y = Lambda, group = as.factor(Perm), linewidth = Lwidth)) +
               geom_line(colour = "darkred", lineend = "round") +
               scale_linewidth(breaks = 1.5:5.5, range = c(0.3, 1.3)) +
               coord_cartesian(ylim = c(0, 0.0012)) +
               geom_vline(xintercept = 0, colour = "blue", linewidth = 0.5, linetype = 2) +
               labs(subtitle = "PR", y  = NULL, x = NULL) +
               tema,
             nrow = 1, left = "Lambda", bottom = "TPI [m]")

Figure 29: Same as Fig. 22 but for the active Pogonomyrmex colonies in 2004 along both 1.5 km segments of dirt roads.

¿Do ant patterns and edaphic associations persist?: resamples 2007 and 2010

In spite of a big decline in Pogonomyrmex spp active nest with time, notoriously lower in 2010 than 2004 and 2007 (Fig. 2), their spatial distribution remained relatively similar on both dirt roads (Fig. 30).

Within those general similarities, the spatial pattern in 2010 looked more aggregated than 2004 and 2007 (even when the power to detect non-random patterns falls with low colony numbers, i.e., the envelopes under the null model, or critical pointwise areas, become wider), probably as a consequence of the important decline of the more homogeneously distributed species (PR) and the increase in proportional contribution of the two species with aggregated spatial patterns (PI and PM).

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

#loop of transects with vector of values for vertical text position of tests results within plots
loop.texts <- data.frame(Main = c("Intensities (LaFad 2004–2010)", "", "", "Intensities (Doble 2004–2010)", "", ""), Lambda = c(0.09, 0.085, 0.08, 0.09, 0.085, 0.08), K = c(-25, -30, -35, -25, -30, -35), PCF = c(0.7, 0.64, 0.58, 0.7, 0.64, 0.58))
loop.ylims <- list(Lambda = list(c(0, 0.09),  c(0, 0.09)), K = list(c(-40, 50),  c(-40, 50)), PCF = list(c(0.5, 1.5),  c(0.5, 1.5)))
loop.colours <- data.frame(Cols = c("#6600FF", "#009933", "#FF6600"))
loop.count <- 0
big.loop <- 0

for(i in PicadasAll) {
  big.loop <- big.loop + 1
  for(j in 1:length(i)) {
    plot(rhohat(unmark(i[[j]]), "x"), legend = FALSE, main = loop.texts$Main[loop.count + j], cex.axis = .8, col = loop.colours$Cols[j], shadecol = alpha(loop.colours$Cols[j], 0.25), ylim = loop.ylims$Lambda[[big.loop]], add = j > 1)
    text(1525, y = loop.texts$Lambda[loop.count + j], paste("p = ", round(cdf.test(unmark(i[[j]]), "x", test = "ks", jitter = F)$p.value, 3)), adj = c(1, 1), col = loop.colours$Cols[j])
  }
  for(j in 1:length(i)) {
    #temp guarda los valores de la función para poder usarlos dos veces (plot + mad.test). fix.n = TRUE ver Baddeley et al 2014 Ecol.Monog.
    temp <- envelope(unmark(i[[j]]), linearK, nsim = 199, nrank = 5, r = seq(0, 350, 10), verbose = FALSE, fix.n = TRUE, savefuns = TRUE, warnWbw = FALSE)
    plot(temp, cbind(obs - theo, lo = lo - theo, hi = hi - theo) ~ r, main = "K-r ~ r", legend = FALSE, cex.axis = 0.8, col = loop.colours$Cols[j], shadecol = alpha(loop.colours$Cols[j], 0.25), ylim = loop.ylims$K[[big.loop]], add = j > 1)
    text(350, y = loop.texts$K[loop.count + j], paste("p = ", round(mad.test(temp, use.theory = TRUE)$p.value, 3)), adj = c(1, 1), col = loop.colours$Cols[j])
    abline(h = 0, lty = 2, col = "red", lwd = 0.5)
  }
  for(j in 1:length(i)) {
    temp <- envelope(unmark(i[[j]]), linearpcf, nsim = 199, nrank = 5, r = seq(0, 350, 10), verbose = FALSE, fix.n = TRUE, savefuns = TRUE, warnWbw = FALSE)
    plot(temp, main = "pcf ~ r", legend = FALSE, cex.axis = 0.8, col = loop.colours$Cols[j], shadecol = alpha(loop.colours$Cols[j], 0.25), ylim = loop.ylims$PCF[[big.loop]], add = j > 1)
    text(350, y = loop.texts$PCF[loop.count + j], paste("p = ", round(mad.test(temp, use.theory = TRUE)$p.value, 3)), adj = c(1, 1), col = loop.colours$Cols[j])
  }
  loop.count <- loop.count + length(i)
}

## Ties were encountered

Figure 30: Same as Fig. 3 but for the three periods in which active colonies were censused in the same two 1.5 km segments of dirt roads: 2004, 2007, and 2010.

In fact, in most cases their variations in linear density along the roads showed important differences to expectations if they were randomly and independently distributed (=CSRI; Fig. 31). The patterns of aggregation at scales of tenths of meters are similar in 2004 and 2007, but in some cases they are even more marked in 2010: in La Fad both PM and PR showed intraspecific patterns clearly aggregated at distances of tenths of meters, as well as PM in Doble.

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

loop.graphs <- list(Main = c("LaFad 2004–2010: PI", "LaFad 2004–2010: PM", "LaFad 2004–2010: PR", "Doble 2004–2010: PI", "Doble 2004–2010: PM", "Doble 2004–2010: PR"), Lambda = list(c(0, 0.045), c(0, 0.028), c(0, 0.07), c(0, 0.085), c(0, 0.025), c(0, 0.06)), K = list(c(-70, 150), c(-100, 250), c(-80, 230), c(-90, 220), c(-120, 200), c(-180, 180)), PCF = list(c(0, 2.5), c(0, 3.5), c(0.3, 2.7), c(0, 2.7), c(0, 2.8), c(0, 3.5)), Lambda.texts = list(c(0.045, 0.04, 0.035), c(0.026, 0.023, 0.02), c(0.066, 0.058, 0.05), c(0.08, 0.07, 0.06), c(0.024, 0.021, 0.018), c(0.056, 0.05, 0.044)), K.texts = list(c(150, 125, 100), c(-15, -50, -80), c(200, 170, 140), c(-10, -40, -70), c(-20, -55, -90), c(-70, -110, -150)), PCF.texts = list(c(2.5, 2.25, 2), c(3.5, 3.1, 2.7), c(2.6, 2.3, 2), c(2.7, 2.4, 2.1), c(2.8, 2.5, 2.2), c(3.4, 3, 2.6)))
loop.spp <- data.frame(Sp = c("PI", "PM", "PR"), Cols = c("#6600FF", "#009933", "#FF6600"))
graph.count <- 0

#intensity graphs
for(i in 1:length(PicadasAll)) { #picadas (panel)
  for(j in 1:length(loop.spp$Sp)) { #especies (plot)
    for(k in 1:length(PicadasAll[[i]])) { #años (curvas)
      graph.count <- (i - 1) * length(loop.spp$Sp) + j
      plot(rhohat(unmark(subset(PicadasAll[[i]][[k]], marks == loop.spp$Sp[j], drop = TRUE)), "x"), legend = FALSE, main = loop.graphs$Main[graph.count], cex.axis = .8, col = loop.spp$Cols[k], shadecol = alpha(loop.spp$Cols[k], 0.25), ylim = loop.graphs$Lambda[[graph.count]], add = k > 1)
      text(1525, y = loop.graphs$Lambda.texts[[graph.count]][k], paste("p = ", round(cdf.test(unmark(subset(PicadasAll[[i]][[k]], marks == loop.spp$Sp[j], drop = TRUE)), "x", test = "ks", jitter = F)$p.value, 3)), adj = c(1, 1), col = loop.spp$Cols[k])
    }
  }
#K graphs
  for(j in 1:length(loop.spp$Sp)) { #especies (plot)
    for(k in 1:length(PicadasAll[[i]])) { #años (curvas)
      graph.count <- (i - 1) * length(loop.spp$Sp) + j
      temp <- envelope(subset(PicadasAll[[i]][[k]], marks == loop.spp$Sp[j]), linearK, nsim = 199, nrank = 5, r = seq(0, 350, 10), verbose = FALSE, fix.n = TRUE, savefuns = TRUE, warnWbw = FALSE)
      plot(temp, cbind(obs - theo, lo = lo - theo, hi = hi - theo) ~ r, main = "K-r ~ r", legend = FALSE, cex.axis = 0.8, col = loop.spp$Cols[k], shadecol = alpha(loop.spp$Cols[k], 0.25), ylim = loop.graphs$K[[graph.count]], add = k > 1)
      text(350, y = loop.graphs$K.texts[[graph.count]][k], paste("p = ", round(mad.test(temp, use.theory = TRUE)$p.value, 3)), adj = c(1, 1), col = loop.spp$Cols[k])
      abline(h = 0, lty = 2, col = "red", lwd = 0.5)
    }
  }
#PCF graphs
  for(j in 1:length(loop.spp$Sp)) { #especies (plot)
    for(k in 1:length(PicadasAll[[i]])) { #años (curvas)
      graph.count <- (i - 1) * length(loop.spp$Sp) + j
      temp <- envelope(subset(PicadasAll[[i]][[k]], marks == loop.spp$Sp[j]), linearpcf, nsim = 199, nrank = 5, r = seq(0, 350, 10), verbose = FALSE, fix.n = TRUE, savefuns = TRUE, warnWbw = FALSE)
      plot(temp, main = "pcf ~ r", legend = FALSE, cex.axis = 0.8, col = loop.spp$Cols[k], shadecol = alpha(loop.spp$Cols[k], 0.25), ylim = loop.graphs$PCF[[graph.count]], add = k > 1)
      text(350, y = loop.graphs$PCF.texts[[graph.count]][k], paste("p = ", round(mad.test(temp, use.theory = TRUE)$p.value, 3)), adj = c(1, 1), col = loop.spp$Cols[k])
    }
  }
}

Figure 31: Same as Fig. 3 but for the three periods in which active colonies were censused in the same two 1.5 km segments of dirt roads: 2004, 2007, and 2010.

The spatial patterns of active colonies observed in 2007 and 2010 are also compatible with Poisson point processes in which the intensity changes along the road segments, with no evidence of further distance-dependent positive or negative interactions among neighbouring colonies (Fig. 32); only in a single case (PI in La Fadina in 2007) active colonies happened to appear more aggregated at very short distances than a Poisson model would expect. The sectors of roads with high and low densities coincide among sampled years for PI and PM but are more variable for PR (particularly after its big decline in active colony density in 2010).

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

loop.graphs <- list(Main = c("LaFad 2004–2010: PI", "LaFad 2004–2010: PM", "LaFad 2004–2010: PR", "Doble 2004–2010: PI", "Doble 2004–2010: PM", "Doble 2004–2010: PR"), K = list(c(-300, 300), c(-400, 600), c(-350, 400), c(-300, 500), c(-300, 550), c(-300, 500)), PCF = list(c(-0.3, 3), c(-0.4, 5), c(-0.3, 3), c(-0.2, 3.5), c(-0.2, 3.5), c(0, 5)), K.texts = list(c(300, 250, 200), c(600, 500, 400), c(400, 325, 250), c(500, 425, 350), c(550, 450, 350), c(500, 425, 350)), PCF.texts = list(c(3, 2.75, 2.5), c(5, 4.5, 4), c(3, 2.75, 2.5), c(3.5, 3.2, 2.9), c(3.5, 3.2, 2.9), c(5, 4.5, 4)))
loop.spp <- data.frame(Sp = c("PI", "PM", "PR"), Cols = c("#6600FF", "#009933", "#FF6600"))
graph.count <- 0

for(i in 1:length(PicadasAll)) { #picadas (panel)

  #intensity graphs
  for(j in 1:length(loop.spp$Sp)) { #especies (plot)
    for(k in 1:length(PicadasAll[[i]])) { #años (curvas)
      graph.count <- (i - 1) * length(loop.spp$Sp) + j
      
      #estimate density with gaussian kernel and fixed sigma and define model for sp j in year k
      est_density <- densityfun(unmark(subset(PicadasAll[[i]][[k]], marks == loop.spp$Sp[j], drop = TRUE)), sigma = 100, dx = 5)
      templppm <- lppm(unmark(subset(PicadasAll[[i]][[k]], marks == loop.spp$Sp[j], drop = TRUE)) ~ log(est_density))
      #plot the model with K-S testing fit of observed against model
      plot(templppm, legend = FALSE, main = loop.graphs$Main[graph.count], cex.axis = .8, style = "w", adjust = 30, col = alpha(loop.spp$Cols[k], 0.3), add = k > 1)
      text(1525, y = 240 - 40 * k, paste("p = ", round(cdf.test(templppm, "x", test = "ks", jitter = F)$p.value, 3)), adj = c(1, 1), col = loop.spp$Cols[k])
    }
  }
  #K graphs
  for(j in 1:length(loop.spp$Sp)) { #especies (plot)
    for(k in 1:length(PicadasAll[[i]])) { #años (curvas)
      graph.count <- (i - 1) * length(loop.spp$Sp) + j
      
       #estimate density with gaussian kernel and fixed sigma and define model for sp j in year k
      est_density <- densityfun(unmark(subset(PicadasAll[[i]][[k]], marks == loop.spp$Sp[j], drop = TRUE)), sigma = 100, dx = 5)
      templppm <- lppm(unmark(subset(PicadasAll[[i]][[k]], marks == loop.spp$Sp[j], drop = TRUE)) ~ log(est_density))     
      #pointwise envelope of inhomogeneous linear K function (λ estimated for each simulated pattern)
      temp <- envelope(templppm, linearKinhom, sigma = 100, r = seq(0, 350, 10), normpower = 2, nsim = 199, nrank = 5, verbose = FALSE, fix.n = TRUE, savefuns = TRUE, warnWbw = FALSE)
      #plot pointwise envelope with (single-staged) MAD test according to linearKinhom [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 = loop.spp$Cols[k], shadecol = alpha(loop.spp$Cols[k], 0.3), ylim = loop.graphs$K[[graph.count]], add = k > 1)
      text(350, y = loop.graphs$K.texts[[graph.count]][k], paste("p = ", round(mad.test(temp, use.theory = FALSE)$p.value, 3)), adj = c(1, 1), col = loop.spp$Cols[k])
      #(two-staged) MAD global test based on linearKinhom [Baddeley et al 2017]
      text(0, y = loop.graphs$K.texts[[graph.count]][k], paste("p = ", round(bits.test(templppm, fun = linearKinhom, sigma = 100, exponent = Inf, rinterval = c(0, 350), nsim = 19, leaveoneout = FALSE, use.theory = FALSE, verbose = FALSE, warnWbw = FALSE)$p.value, 3)), adj = c(0, 1), col = loop.spp$Cols[k], font = 2)
      abline(h = 0, lty = 2, col = "red", lwd = 0.5)
    }
  }
  #PCF graphs
  for(j in 1:length(loop.spp$Sp)) { #especies (plot)
    for(k in 1:length(PicadasAll[[i]])) { #años (curvas)
      graph.count <- (i - 1) * length(loop.spp$Sp) + j
      
       #estimate density with gaussian kernel and fixed sigma and define model for sp j in year k
      est_density <- densityfun(unmark(subset(PicadasAll[[i]][[k]], marks == loop.spp$Sp[j], drop = TRUE)), sigma = 100, dx = 5)
      templppm <- lppm(unmark(subset(PicadasAll[[i]][[k]], marks == loop.spp$Sp[j], drop = TRUE)) ~ log(est_density))     
      #pointwise envelope of inhomogeneous linear K function (λ estimated for each simulated pattern)
      temp <- envelope(templppm, linearpcfinhom, sigma = 100, r = seq(0, 350, 10), normpower = 2, nsim = 199, nrank = 5, verbose = FALSE, fix.n = TRUE, savefuns = TRUE, warnWbw = FALSE)
      #plot pointwise envelope with (single-staged) MAD test according to linearKinhom [conservative according to Baddeley et al 2017]
      plot(temp, cbind(obs - mmean + 1, lo = lo - mmean + 1, hi = hi - mmean + 1) ~ r, main = "pcfinhom-r ~ r", legend = FALSE, cex.axis = 0.8, col = loop.spp$Cols[k], shadecol = alpha(loop.spp$Cols[k], 0.3), ylim = loop.graphs$PCF[[graph.count]], add = k > 1)
      text(350, y = loop.graphs$PCF.texts[[graph.count]][k], paste("p = ", round(mad.test(temp, use.theory = FALSE)$p.value, 3)), adj = c(1, 1), col = loop.spp$Cols[k])
      #(two-staged) MAD global test based on linearKinhom [Baddeley et al 2017]
      text(0, y = loop.graphs$PCF.texts[[graph.count]][k], paste("p = ", round(bits.test(templppm, fun = linearpcfinhom, sigma = 100, exponent = Inf, rinterval = c(0, 350), nsim = 19, leaveoneout = FALSE, use.theory = FALSE, verbose = FALSE, warnWbw = FALSE)$p.value, 3)), adj = c(0, 1), col = loop.spp$Cols[k], font = 2)
      abline(h = 1, lty = 2, col = "red", lwd = 0.5)
    }
  }
}

## Error in (function (x, bw = "nrd0", adjust = 1, kernel = c("gaussian",  : 
##   need at least 2 points to select a bandwidth automatically
## [retrying]

## Error in (function (x, bw = "nrd0", adjust = 1, kernel = c("gaussian",  : 
##   need at least 2 points to select a bandwidth automatically
## [retrying]

## Error in (function (x, bw = "nrd0", adjust = 1, kernel = c("gaussian",  : 
##   need at least 2 points to select a bandwidth automatically
## [retrying]

## Error in (function (x, bw = "nrd0", adjust = 1, kernel = c("gaussian",  : 
##   need at least 2 points to select a bandwidth automatically
## [retrying]

## Error in (function (x, bw = "nrd0", adjust = 1, kernel = c("gaussian",  : 
##   need at least 2 points to select a bandwidth automatically
## [retrying]

## Error in (function (x, bw = "nrd0", adjust = 1, kernel = c("gaussian",  : 
##   need at least 2 points to select a bandwidth automatically
## [retrying]

## Error in (function (x, bw = "nrd0", adjust = 1, kernel = c("gaussian",  : 
##   need at least 2 points to select a bandwidth automatically
## [retrying]

## Error in (function (x, bw = "nrd0", adjust = 1, kernel = c("gaussian",  : 
##   need at least 2 points to select a bandwidth automatically
## [retrying]

## Error in (function (x, bw = "nrd0", adjust = 1, kernel = c("gaussian",  : 
##   need at least 2 points to select a bandwidth automatically
## [retrying]

## Error in (function (x, bw = "nrd0", adjust = 1, kernel = c("gaussian",  : 
##   need at least 2 points to select a bandwidth automatically
## [retrying]

## Error in (function (x, bw = "nrd0", adjust = 1, kernel = c("gaussian",  : 
##   need at least 2 points to select a bandwidth automatically
## [retrying]

## Error in (function (x, bw = "nrd0", adjust = 1, kernel = c("gaussian",  : 
##   need at least 2 points to select a bandwidth automatically
## [retrying]

Figure 32: Same as Fig. 4 but for the three periods in which active colonies were censused: 2004, 2007, and 2010.

Road sectors associated with the prevalence of each of the species persist, particularly between 2004 and 2007, when linear densities of PM and PR were similar or slightly lower and PI slightly increased (Fig. 33). In 2010, with notoriously lower densities of all species, the sequence of species-specific sectors was even clearer, mostly because the more distributed species, PR, reduced its density. Notice also that some previous distinct sectors were “lost” as a consequence of no active nests remaining there.

par(mar = c(0, 2, 1, 0.5))
layout(1:8, heights = c(0.3, 1, 1, 1))

loop.spp <- data.frame(Sp = c("PI", "PM", "PR"), Cols = c("darkgreen", "darkblue", "darkred"))

for(i in 1:length(PicadasAll)) { #picadas (panel)
  plot.new()
  mtext(text = names(PicadasAll)[i], side = 1, cex = 1.25, line = -0.5)
  for(j in 1:length(PicadasAll[[i]])) { #años (plot)
    temp <- relrisk(PicadasAll[[i]][[j]], sigma = 25, finespacing = T)
    for(k in 1:length(temp)) { #sp (curvas)
      plot(temp[[k]], style = "width", scale = 50, legend = F, main = "", col = alpha(loop.spp$Cols[k], 0.5), add = k > 1)
      plot(split(PicadasAll[[i]][[j]])[[k]], add = T, col = loop.spp$Cols[k], pch = 19, size = 0.5, lty = 0)
      mtext(c(2004, 2007, 2010)[j], side = 2, las = 1, adj = 0.5)
    }
  }
}

Figure 33: Same as Fig. 6 but along time, adding the resampling of active colonies in the same road segments in 2007 and 2010.

The association between the intensity of the point process corresponding to each species and both edaphic factors remained consistent between 2004 and 2010 in spite of their changes in density, following the patterns already described (e.g., PI with low permeability vs PM with the most permeable soil surfaces; Fig. 34.

par(mar = c(2, 2, 1, .5))
layout(matrix(1:6, 3, 2))

loop.graphs <- list(Main = c("LaFad 2004–2010: PI", "LaFad 2004–2010: PM", "LaFad 2004–2010: PR", "Doble 2004–2010: PI", "Doble 2004–2010: PM", "Doble 2004–2010: PR"), LambdaPerm = list(c(0, 0.15), c(0, 0.04), c(0, 0.1), c(0, 0.1), c(0, 0.05), c(0, 0.05)), LambdaPerm.texts = list(c(0.15, 0.14, 0.13), c(0.04, 0.035, 0.03), c(0.1, 0.092, 0.08), c(0.1, 0.09, 0.08), c(0.05, 0.045, 0.04), c(0.05, 0.045, 0.04)), LambdaAlt = list(c(0, 0.15), c(0, 0.04), c(0, 0.1), c(0, 0.1), c(0, 0.05), c(0, 0.05)), LambdaAlt.texts = list(c(0.15, 0.14, 0.13), c(0.04, 0.035, 0.03), c(0.1, 0.092, 0.08), c(0.1, 0.09, 0.08), c(0.05, 0.045, 0.04), c(0.05, 0.045, 0.04)))
loop.spp <- data.frame(Sp = c("PI", "PM", "PR"), Cols = c("#6600FF", "#009933", "#FF6600"))
graph.count <- 0

#intensity graphs x Perm
for(i in 1:length(PicadasAll)) { #picadas (panel)
  for(j in 1:length(loop.spp$Sp)) { #especies (plot)
    for(k in 1:length(PicadasAll[[i]])) { #años (curvas)
      graph.count <- (i - 1) * length(loop.spp$Sp) + j
      plot(rhohat(unmark(subset(PicadasAll[[i]][[k]], marks == loop.spp$Sp[j], drop = TRUE)), Both_Perm_Smooth[[i]]), legend = FALSE, main = loop.graphs$Main[graph.count], cex.axis = .8, col = loop.spp$Cols[k], shadecol = alpha(loop.spp$Cols[k], 0.25), ylim = loop.graphs$LambdaPerm[[graph.count]], add = k > 1)
      abline(v = c(0.7, -0.86), lty = 2, col = "blue", lwd = 0.5)
      text(0.4, y = loop.graphs$LambdaPerm.texts[[graph.count]][k], paste("p = ", round(cdf.test(unmark(subset(PicadasAll[[i]][[k]], marks == loop.spp$Sp[j], drop = TRUE)), Both_Perm_Smooth[[i]], test = "ks", jitter = F)$p.value, 3)), adj = c(0, 1), col = loop.spp$Cols[k])
    }
  }
  #intensity graphs x AltTPI
  for(j in 1:length(loop.spp$Sp)) { #especies (plot)
    for(k in 1:length(PicadasAll[[i]])) { #años (curvas)
      graph.count <- (i - 1) * length(loop.spp$Sp) + j
      plot(rhohat(unmark(subset(PicadasAll[[i]][[k]], marks == loop.spp$Sp[j], drop = TRUE)), Both_Alt_TPI[[i]]), legend = FALSE, main = loop.graphs$Main[graph.count], cex.axis = .8, col = loop.spp$Cols[k], shadecol = alpha(loop.spp$Cols[k], 0.25), ylim = loop.graphs$LambdaAlt[[graph.count]], add = k > 1)
      abline(v = 0, lty = 2, col = "blue", lwd = 0.5)
      text(0.1, y = loop.graphs$LambdaAlt.texts[[graph.count]][k], paste("p = ", round(cdf.test(unmark(subset(PicadasAll[[i]][[k]], marks == loop.spp$Sp[j], drop = TRUE)), Both_Alt_TPI[[i]], test = "ks", jitter = F)$p.value, 3)), adj = c(0, 1), col = loop.spp$Cols[k])
    }
  }
}

Figure 34: Resource selection functions, or intensities of heterogeneous Poisson point process models for the active colonies of each Pogonomyrmex spp. in 2004, 2007, and 2010 (solid line: fixed-bandwith kernel density estimator; shaded areas: pointwise confidence intervals) as functions of soil permeability (SPI; left) and relative altitude (Topography Index [m]; right) along the two segments of dirt roads (La Fadina and Doble). Discontinuous vertical lines show (left) proposed infiltration time thresholds between three typical habitats (see Fig. 11) and (right) neutral altitude (TPI = 0). P values result from Kolmogorov-Smirnov goodness-of-fit tests against (homogeneous = no relationship) CSR models with same mean intensities.

Poisson point pattern models

In spite of density changes and their spatial consequences, the structure of the most parsimonious Poisson models for 2007 and 2010, with the two segments considered as spatial replicates of the same underlying spatial processes, is the same for both years and very similar to that already described for 2004: both predictive edaphic factors (soil permeability always a more important influence than topography), in these cases interacting, with distinct parameters for each of the three Pogonomyrmex species:

Ambas2007 <- hyperframe(Nidos = solist(quadscheme(LaFad2007_Pogos), quadscheme(Doble2007_Pogos)),
                        Perm = solist(LaFad_Perm_Smooth, Doble_Perm_Smooth),
                        Alt = solist(LaFad_Alt_TPI, Doble_Alt_TPI), 
                        Picada = as.factor(c("LF", "D")))

Ambas2010 <- hyperframe(Nidos = solist(quadscheme(LaFad2010_Pogos), quadscheme(Doble2010_Pogos)),
                        Perm = solist(LaFad_Perm_Smooth, Doble_Perm_Smooth),
                        Alt = solist(LaFad_Alt_TPI, Doble_Alt_TPI), 
                        Picada = as.factor(c("LF", "D")))
                        
cat("2007\n")
#(1) agregando a la permeabilidad o a la altura como predictores lineales de la variación espacial de la intensidad  (i.e., afectando a las 3 spp por igual) o con interacción entre edáficas y especies [= test de interacción], de modo lineal o cuadrático
cat("\n[1] Model with heterogeneous intensities varying in the same way linearly with each edaphic factor (Models 1) or with species-specific linear (Models 2) or quadratic (Models 3) relationships\n")
anova(mppm(data = Ambas2007, formula = Nidos ~ marks + Perm),
      mppm(data = Ambas2007, formula = Nidos ~ marks * Perm),
      mppm(data = Ambas2007, formula = Nidos ~ marks * polynom(Perm, 2)),
      test = "LR")
anova(mppm(data = Ambas2007, formula = Nidos ~ marks + Alt),
      mppm(data = Ambas2007, formula = Nidos ~ marks * Alt),
      mppm(data = Ambas2007, formula = Nidos ~ marks * polynom(Alt, 2)),
      test = "LR")

#(2) agregando a la topografía sobre el modelo que ya tiene a permeabilidad * especie
cat("\n[2] Model with species-specific quadratic effects of soil permeability (Model 1), with an additional independent linear effect of topography (Model 2), with both species-specific effects interacting (Model 3), and adding a quadratic term for topography (Model 4)\n")
anova(mppm(data = Ambas2007, formula = Nidos ~ marks * polynom(Perm, 2)),
      mppm(data = Ambas2007, formula = Nidos ~ marks * polynom(Perm, 2) + marks * Alt),
      mppm(data = Ambas2007, formula = Nidos ~ marks * polynom(Perm, 2) * Alt),
      mppm(data = Ambas2007, formula = Nidos ~ marks * polynom(Perm, 2) * polynom(Alt, 2)),
      test = "LR")

cat("\n2010\n")
#(1) agregando a la permeabilidad o a la altura como predictores lineales de la variación espacial de la intensidad  (i.e., afectando a las 3 spp por igual) o con interacción entre edáficas y especies [= test de interacción], de modo lineal o cuadrático
cat("\n[1] Model with heterogeneous intensities varying in the same way linearly with each edaphic factor (Models 1) or with species-specific linear (Models 2) or quadratic (Models 3) relationships\n")
anova(mppm(data = Ambas2010, formula = Nidos ~ marks + Perm),
      mppm(data = Ambas2010, formula = Nidos ~ marks * Perm),
      mppm(data = Ambas2010, formula = Nidos ~ marks * polynom(Perm, 2)),
      test = "LR")
anova(mppm(data = Ambas2010, formula = Nidos ~ marks + Alt),
      mppm(data = Ambas2010, formula = Nidos ~ marks * Alt),
      mppm(data = Ambas2010, formula = Nidos ~ marks * polynom(Alt, 2)),
      test = "LR")

#(2) agregando a la topografía sobre el modelo que ya tiene a permeabilidad * especie
cat("\n[2] Model with species-specific quadratic effects of soil permeability (Model 1), with an additional independent linear effect of topography (Model 2), with both species-specific effects interacting (Model 3), and adding a quadratic term for topography (Model 4)\n")
anova(mppm(data = Ambas2010, formula = Nidos ~ marks * polynom(Perm, 2)),
      mppm(data = Ambas2010, formula = Nidos ~ marks * polynom(Perm, 2) + marks * Alt),
      mppm(data = Ambas2010, formula = Nidos ~ marks * polynom(Perm, 2) * Alt),
      mppm(data = Ambas2010, formula = Nidos ~ marks * polynom(Perm, 2) * polynom(Alt, 2)),
      test = "LR")

PoissonAmbas2007 <- mppm(data = Ambas2007, formula = Nidos ~ marks * polynom(Perm, 2) * Alt)
PoissonAmbas2010 <- mppm(data = Ambas2010, formula = Nidos ~ marks * polynom(Perm, 2) * Alt)

## 2007
## 
## [1] Model with heterogeneous intensities varying in the same way linearly with each edaphic factor (Models 1) or with species-specific linear (Models 2) or quadratic (Models 3) relationships
## Analysis of Deviance Table
## 
## Model 1: Nidos ~ marks + Perm     Poisson
## Model 2: Nidos ~ marks * Perm     Poisson
## Model 3: Nidos ~ marks * polynom(Perm, 2)     Poisson
##   Npar Df Deviance  Pr(>Chi)    
## 1    4                          
## 2    6  2   48.077 3.632e-11 ***
## 3    9  3   10.876   0.01242 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Analysis of Deviance Table
## 
## Model 1: Nidos ~ marks + Alt      Poisson
## Model 2: Nidos ~ marks * Alt      Poisson
## Model 3: Nidos ~ marks * polynom(Alt, 2)      Poisson
##   Npar Df Deviance     Pr(>Chi)    
## 1    4                             
## 2    6  2   27.872 0.0000008862 ***
## 3    9  3    6.428      0.09254 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## [2] Model with species-specific quadratic effects of soil permeability (Model 1), with an additional independent linear effect of topography (Model 2), with both species-specific effects interacting (Model 3), and adding a quadratic term for topography (Model 4)
## Analysis of Deviance Table
## 
## Model 1: Nidos ~ marks * polynom(Perm, 2)     Poisson
## Model 2: Nidos ~ marks * polynom(Perm, 2) + marks * Alt   Poisson
## Model 3: Nidos ~ marks * polynom(Perm, 2) * Alt   Poisson
## Model 4: Nidos ~ marks * polynom(Perm, 2) * polynom(Alt, 2)   Poisson
##   Npar Df Deviance Pr(>Chi)   
## 1    9                        
## 2   12  3   11.471 0.009434 **
## 3   18  6   22.011 0.001205 **
## 4   27  9   13.805 0.129443   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 2010
## 
## [1] Model with heterogeneous intensities varying in the same way linearly with each edaphic factor (Models 1) or with species-specific linear (Models 2) or quadratic (Models 3) relationships
## Analysis of Deviance Table
## 
## Model 1: Nidos ~ marks + Perm     Poisson
## Model 2: Nidos ~ marks * Perm     Poisson
## Model 3: Nidos ~ marks * polynom(Perm, 2)     Poisson
##   Npar Df Deviance       Pr(>Chi)    
## 1    4                               
## 2    6  2   40.746 0.000000001419 ***
## 3    9  3   14.693       0.002098 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Analysis of Deviance Table
## 
## Model 1: Nidos ~ marks + Alt      Poisson
## Model 2: Nidos ~ marks * Alt      Poisson
## Model 3: Nidos ~ marks * polynom(Alt, 2)      Poisson
##   Npar Df Deviance    Pr(>Chi)    
## 1    4                            
## 2    6  2  24.5713 0.000004618 ***
## 3    9  3   2.6262      0.4529    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## [2] Model with species-specific quadratic effects of soil permeability (Model 1), with an additional independent linear effect of topography (Model 2), with both species-specific effects interacting (Model 3), and adding a quadratic term for topography (Model 4)
## Analysis of Deviance Table
## 
## Model 1: Nidos ~ marks * polynom(Perm, 2)     Poisson
## Model 2: Nidos ~ marks * polynom(Perm, 2) + marks * Alt   Poisson
## Model 3: Nidos ~ marks * polynom(Perm, 2) * Alt   Poisson
## Model 4: Nidos ~ marks * polynom(Perm, 2) * polynom(Alt, 2)   Poisson
##   Npar Df Deviance Pr(>Chi)  
## 1    9                       
## 2   12  3   9.2806  0.02578 *
## 3   18  6  13.1759  0.04033 *
## 4   27  9   6.5827  0.68047  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Predictions from those Poisson point process models on how their intensity change along the dirt roads according to the measured edaphic variables (soil permeability and topography) are strikingly similar for the three compared years in spite of some structural differences in the selected model for different years (Fig. 35). For two of the three species, linear and quadratic effects resulted in rather “parallel” predictions. Succinctly,

• PI increased strongly in soils with low permeability (topography has marginal importance, increasing slightly in road depressions when permeability has higher values),
• PM increased strongly with soil permeability and with higher relative elevations, and
• PR always showed a pronounced effect of the quadratic term for permeability, peaking at typical values for the algarrobal and decreasing towards both lower and higher values, and a tendency to avoid high relative topography in 2004 and 2007.

# predict y effectfun no funcionan con mppm sobre lpp: hay que hacerlo como producto matricial con la matriz de diseño
# predicciones en la escala de los datos (simil effectfun)

# efectos: perm (con alt = quantile(25%, median, 75%))
temp.data <- data.frame(Alt = rep(rep(quantile(c(LaFad_Alt_TPI(0:1500, 0), Doble_Alt_TPI(0:1500, 0)), probs = c(0.25, 0.5, 0.75)), each = 50), 3), 
                        Perm = rep(rep(seq(min(c(LaFad_Perm_Smooth(0:1500, 0), Doble_Perm_Smooth(0:1500, 0))), max(c(LaFad_Perm_Smooth(0:1500, 0), Doble_Perm_Smooth(0:1500, 0))), length.out = 50), 3), 3),
                        Year = rep(c("2004", "2007", "2010"), each = 150), 
                        Lwidth = rep(rep(1:3, each = 50), 3)) %>%
             bind_rows("PI" = ., "PM" = ., "PR" = ., .id = "Sp")

#lambdas de cada año (potencialmente diferentes matrices de diseño de sus modelos)
coef.data <- coef(Poisson2Ambas2004)
temp.matrix <- temp.data %>%
  filter(Year == "2004") %>%
  mutate(V01 = TRUE, V02 = Sp == "PM", V03 = Sp == "PR", V04 = Perm, V05 = Perm ^ 2, V06 = Alt, V07 = Alt ^ 2, V08 = V02 * V04, V09 = V03 * V04, V10 = V02 * V05, V11 = V03 * V05, V12 = V02 * V06, V13 = V03 * V06, V14 = V02 * V07, V15 = V03 * V07, .keep = "none")
temp.data$Lambda[temp.data$Year == "2004"] <- exp(as.matrix(temp.matrix) %*% coef.data)[,1]

coef.data <- coef(PoissonAmbas2007)
temp.matrix <- temp.data %>%
  filter(Year == "2007") %>%
  mutate(V01 = TRUE, V02 = Sp == "PM", V03 = Sp == "PR", V04 = Perm, V05 = Perm ^ 2, V06 = Alt, V07 = V02 * V04, V08 = V03 * V04, V09 = V02 * V05, V10 = V03 * V05, V11 = V02 * V06, V12 = V03 * V06, V13 = V04 * V06, V14 = V05 * V06, V15 = V02 * V13, V16 = V03 * V13, V17 = V02 * V14, V18 = V03 * V14,.keep = "none")
temp.data$Lambda[temp.data$Year == "2007"] <- exp(as.matrix(temp.matrix) %*% coef.data)[,1]

coef.data <- coef(PoissonAmbas2010)
temp.matrix <- temp.data %>%
  filter(Year == "2010") %>%
  mutate(V01 = TRUE, V02 = Sp == "PM", V03 = Sp == "PR", V04 = Perm, V05 = Perm ^ 2, V06 = Alt, V07 = V02 * V04, V08 = V03 * V04, V09 = V02 * V05, V10 = V03 * V05, V11 = V02 * V06, V12 = V03 * V06, V13 = V04 * V06, V14 = V05 * V06, V15 = V02 * V13, V16 = V03 * V13, V17 = V02 * V14, V18 = V03 * V14,.keep = "none")
temp.data$Lambda[temp.data$Year == "2010"] <- exp(as.matrix(temp.matrix) %*% coef.data)[,1]

grid.arrange(ggplot(subset(temp.data, Sp == "PI"), aes(x = Perm, y = Lambda, group = interaction(Year, Alt), linewidth = Lwidth, colour = as.factor(Year))) + geom_line(lineend = "round") + scale_linewidth(breaks = 1.5:5.5, range = c(0.3, 1.3)) + scale_colour_manual(values = c("#6600FF75", "#00993375", "#FF660075")) + coord_cartesian(ylim = c(0, 0.004)) + geom_vline(xintercept = c(0.7, -0.86), colour = "blue", linewidth = 0.5, linetype = 2) + labs(title = "PI", y  = "Lambda", x = "Perm") + scale_y_continuous(labels = label_number()) + tema,
             ggplot(subset(temp.data, Sp == "PM"), aes(x = Perm, y = Lambda, group = interaction(Year, Alt), linewidth = Lwidth, colour = as.factor(Year))) + geom_line(lineend = "round") + scale_linewidth(breaks = 1.5:5.5, range = c(0.3, 1.3)) + scale_colour_manual(values = c("#6600FF75", "#00993375", "#FF660075")) + coord_cartesian(ylim = c(0, 0.0008)) + geom_vline(xintercept = c(0.7, -0.86), colour = "blue", linewidth = 0.5, linetype = 2) + labs(title = "PM", y  = NULL, x = "Perm") + scale_y_continuous(labels = label_number()) + tema,
             ggplot(subset(temp.data, Sp == "PR"), aes(x = Perm, y = Lambda, group = interaction(Year, Alt), linewidth = Lwidth, colour = as.factor(Year))) + geom_line(lineend = "round") + scale_linewidth(breaks = 1.5:5.5, range = c(0.3, 1.3)) + scale_colour_manual(values = c("#6600FF75", "#00993375", "#FF660075")) + coord_cartesian(ylim = c(0, 0.0015)) + geom_vline(xintercept = c(0.7, -0.86), colour = "blue", linewidth = 0.5, linetype = 2) + labs(title = "PR", y  = NULL, x = "Perm") + scale_y_continuous(labels = label_number()) + tema,
             nrow = 1)

# efectos: alt (con perm = {25%, median, 75%})
temp.data <- data.frame(Perm = rep(rep(quantile(c(LaFad_Perm_Smooth(0:1500, 0), Doble_Perm_Smooth(0:1500, 0)), probs = c(0.25, 0.5, 0.75)), each = 50), 3), 
                        Alt = rep(rep(seq(min(c(LaFad_Alt_TPI(0:1500, 0), Doble_Alt_TPI(0:1500, 0))), max(c(LaFad_Alt_TPI(0:1500, 0), Doble_Alt_TPI(0:1500, 0))), length.out = 50), 3), 3),
                        Year = rep(c("2004", "2007", "2010"), each = 150), 
                        Lwidth = rep(rep(1:3, each = 50), 3)) %>%
             bind_rows("PI" = ., "PM" = ., "PR" = ., .id = "Sp")

#lambdas de cada año (potencialmente diferentes matrices de diseño de sus modelos)
coef.data <- coef(Poisson2Ambas2004)
temp.matrix <- temp.data %>%
  filter(Year == "2004") %>%
  mutate(V01 = TRUE, V02 = Sp == "PM", V03 = Sp == "PR", V04 = Perm, V05 = Perm ^ 2, V06 = Alt, V07 = Alt ^ 2, V08 = V02 * V04, V09 = V03 * V04, V10 = V02 * V05, V11 = V03 * V05, V12 = V02 * V06, V13 = V03 * V06, V14 = V02 * V07, V15 = V03 * V07, .keep = "none")
temp.data$Lambda[temp.data$Year == "2004"] <- exp(as.matrix(temp.matrix) %*% coef.data)[,1]

coef.data <- coef(PoissonAmbas2007)
temp.matrix <- temp.data %>%
  filter(Year == "2007") %>%
  mutate(V01 = TRUE, V02 = Sp == "PM", V03 = Sp == "PR", V04 = Perm, V05 = Perm ^ 2, V06 = Alt, V07 = V02 * V04, V08 = V03 * V04, V09 = V02 * V05, V10 = V03 * V05, V11 = V02 * V06, V12 = V03 * V06, V13 = V04 * V06, V14 = V05 * V06, V15 = V02 * V13, V16 = V03 * V13, V17 = V02 * V14, V18 = V03 * V14,.keep = "none")
temp.data$Lambda[temp.data$Year == "2007"] <- exp(as.matrix(temp.matrix) %*% coef.data)[,1]

coef.data <- coef(PoissonAmbas2010)
temp.matrix <- temp.data %>%
  filter(Year == "2010") %>%
  mutate(V01 = TRUE, V02 = Sp == "PM", V03 = Sp == "PR", V04 = Perm, V05 = Perm ^ 2, V06 = Alt, V07 = V02 * V04, V08 = V03 * V04, V09 = V02 * V05, V10 = V03 * V05, V11 = V02 * V06, V12 = V03 * V06, V13 = V04 * V06, V14 = V05 * V06, V15 = V02 * V13, V16 = V03 * V13, V17 = V02 * V14, V18 = V03 * V14,.keep = "none")
temp.data$Lambda[temp.data$Year == "2010"] <- exp(as.matrix(temp.matrix) %*% coef.data)[,1]

grid.arrange(ggplot(subset(temp.data, Sp == "PI"), aes(x = Alt, y = Lambda, group = interaction(Year, Perm), linewidth = Lwidth, colour = as.factor(Year))) + geom_line(lineend = "round") + scale_linewidth(breaks = 1.5:5.5, range = c(0.3, 1.3)) + scale_colour_manual(values = c("#6600FF75", "#00993375", "#FF660075")) + coord_cartesian(ylim = c(0, 0.004)) + geom_vline(xintercept = 0, colour = "blue", linewidth = 0.5, linetype = 2) + labs(title = "PI", y  = "Intensity", x = "TPI") + scale_y_continuous(labels = label_number()) + tema,
             ggplot(subset(temp.data, Sp == "PM"), aes(x = Alt, y = Lambda, group = interaction(Year, Perm), linewidth = Lwidth, colour = as.factor(Year))) + geom_line(lineend = "round") + scale_linewidth(breaks = 1.5:5.5, range = c(0.3, 1.3)) + scale_colour_manual(values = c("#6600FF75", "#00993375", "#FF660075")) + coord_cartesian(ylim = c(0, 0.0008)) + geom_vline(xintercept = 0, colour = "blue", linewidth = 0.5, linetype = 2) + labs(title = "PM", y  = NULL, x = "TPI") + scale_y_continuous(labels = label_number()) + tema,
             ggplot(subset(temp.data, Sp == "PR"), aes(x = Alt, y = Lambda, group = interaction(Year, Perm), linewidth = Lwidth, colour = as.factor(Year))) + geom_line(lineend = "round") + scale_linewidth(breaks = 1.5:5.5, range = c(0.3, 1.3)) + scale_colour_manual(values = c("#6600FF75", "#00993375", "#FF660075")) + coord_cartesian(ylim = c(0, 0.0015)) + geom_vline(xintercept = 0, colour = "blue", linewidth = 0.5, linetype = 2) + labs(title = "PR", y  = NULL, x = "TPI") + scale_y_continuous(labels = label_number()) + tema,
             nrow = 1)

Figure 35: Partial effects of soil permeability (above) and topographic position (below), according to the selected Poisson models, for three different values of the other variable (25%, 50% and 75% quantiles, in order of increasing line width) on the intensity of the point process determining the location of active nests of each Pogonomyrmex ant species in 2004, 2007, and 2010 along a 1.5 km segment of La Fadina dirt road (Ñacuñán Reserve, central Monte desert).

— A single multispecies Poisson model across time and space

Building on conclusions on previous partial models, we can fit a single Poisson model for the three Pogonomyrmex species, assuming roads as replicates and explicitly including temporal variation to account for every sampled point pattern.

Along the fitting process below we can see that [1] the density of active Pogonomyrmex colonies changed significantly with time, with some (minor) variation among the three species [2] but in a similar way for both roads, and that soil characteristics [3] were both strong predictors of the spatial variations of each species (permeability being the strongest), with [4] interacting effects that did not change significantly with time.

#Modelos Poisson con transectas como réplicas y agregando el tiempo como factor cruzado con Spp 
# Picada and Year should be random, but n = 2 and 3 respectively
# Years on same Road should/could be correlated
# Problem: deviance tests (anova.mppm) do not work with random factors…
# (BTW, if mixed check correct full model in lme ¿sería medidas repetidas con random = ~ 1 | Picada, ~ Year | Picada, o ~ efectos | Picada/Year ?)

Todo <- hyperframe(Nidos = solist(quadscheme(LaFad2004_Pogos), quadscheme(LaFad2007_Pogos), quadscheme(LaFad2010_Pogos), quadscheme(Doble2004_Pogos), quadscheme(Doble2007_Pogos), quadscheme(Doble2010_Pogos)),
                  Perm = solist(LaFad_Perm_Smooth, LaFad_Perm_Smooth, LaFad_Perm_Smooth, Doble_Perm_Smooth, Doble_Perm_Smooth, Doble_Perm_Smooth),
                  Alt = solist(LaFad_Alt_TPI, LaFad_Alt_TPI, LaFad_Alt_TPI, Doble_Alt_TPI, Doble_Alt_TPI, Doble_Alt_TPI), 
                  Picada = as.factor(c("LF", "LF", "LF", "D", "D", "D")),
                  Year = as.factor(c("2004", "2007", "2010", "2004", "2007", "2010")))

#comparación de modelos Poisson anidados, stepwise forward
#(1) comparación modelo Poisson homogéneo por especie (= intensidad fija) vs. intensidad variable en el tiempo por spp y/o por picada
cat("[1] Homogeneous models with constant per species intensity (Model 1), with density varying with time for all species in the same way (Model 2) or with species-specific temporal changes (Model 3) \n")
PoissonTodosNull <- mppm(data = Todo, formula = Nidos ~ marks)
PoissonTodosYear <- mppm(data = Todo, formula = Nidos ~ marks * Year)

anova(PoissonTodosNull, 
      mppm(data = Todo, formula = Nidos ~ marks + Year), 
      PoissonTodosYear,
      test = "LRT")

#(2) comparación modelo homogeneo general o picada-dependiente
cat("\n[2] Species-specific intensities changing with time (Model 1), with parallel changes in each road (Model 2: intercepts model), or with different temporal changes for each road (Model 3: slopes model) \n")
PoissonTodosYearRoad <- mppm(data = Todo, formula = Nidos ~ Picada * (marks * Year))

anova(PoissonTodosYear,
      mppm(data = Todo, formula = Nidos ~ Picada + (marks * Year)),
      PoissonTodosYearRoad,
      test = "LRT")

#(3) comparación modelo homogeneo con heterogeneos dependiendo del suelo
cat("\n[3] Homogeneous species-specific intensities changing with time (Model 1) or with intensity variations depending on a single, time-invariant linear relationship with the edaphic variable (Models 2), or in species-specific ways (Models 3), or with a quadratic relationship (Models 4) that changes with time (Models 5) \n")
PoissonTodosPerm <- mppm(data = Todo, formula = Nidos ~ (marks * Year) + (marks * polynom(Perm, 2)))
PoissonTodosAlt <- mppm(data = Todo, formula = Nidos ~ (marks * Year) + (marks * polynom(Alt, 2)))

anova(PoissonTodosYear, 
      mppm(data = Todo, formula = Nidos ~ (marks * Year) + Perm),
      mppm(data = Todo, formula = Nidos ~ (marks * Year) + (marks * Perm)),
      PoissonTodosPerm,
      mppm(data = Todo, formula = Nidos ~ marks * Year * polynom(Perm, 2)),
      test = "LRT")
      
anova(mppm(data = Todo, formula = Nidos ~ marks * Year), 
      mppm(data = Todo, formula = Nidos ~ (marks * Year) + Alt),
      mppm(data = Todo, formula = Nidos ~ (marks * Year) + (marks * Alt)),
      PoissonTodosAlt,
      mppm(data = Todo, formula = Nidos ~ marks * Year * polynom(Alt, 2)),
      test = "LRT")

#(4) agregando a la topografía sobre el modelo que ya tiene a permeabilidad y rechequeando constancia temporal
cat("\n[4] Model with species-specific quadratic effects of soil permeability (Model 1), with an additional independent quadratic effect of topography (Model 2), with both species-specific effects interacting (Model 3), and effects that change with time (Model 4)\n")
PoissonTodosBothIndep <- mppm(data = Todo, formula = Nidos ~ (marks * Year) + (marks * polynom(Perm, 2)) + (marks * polynom(Alt, 2)))
PoissonTodosBothInter <- mppm(data = Todo, formula = Nidos ~ (marks * Year) + (marks * polynom(Perm, 2) * polynom(Alt, 2)))
PoissonTodosFull <- mppm(data = Todo, formula = Nidos ~ marks * Year * polynom(Perm, 2) * polynom(Alt, 2))

anova(PoissonTodosPerm,
      PoissonTodosBothIndep,
      PoissonTodosBothInter,
      PoissonTodosFull,
      test = "LR")

#(5) best model backwards?
# tomando a los factores fijos cruzados [no del todo correcto porque la permeabilidad y altura no están realmente replicados entre años, solo sus efectos], "se van" todas las interacciones que involucran Year con Alt o Perm => ¿efectos de edáficas constantes en el tiempo o solo artefacto por sus valores constantes?
#cat("\n Stepwise model selection [backward+forward] based on AIC \n")
#step(mppm(data = Todo, formula = Nidos ~ marks * Year * polynom(Alt,2) * polynom(Perm,2)), direction = "both")

## [1] Homogeneous models with constant per species intensity (Model 1), with density varying with time for all species in the same way (Model 2) or with species-specific temporal changes (Model 3) 
## Analysis of Deviance Table
## 
## Model 1: Nidos ~ marks    Poisson
## Model 2: Nidos ~ marks + Year     Poisson
## Model 3: Nidos ~ marks * Year     Poisson
##   Npar Df Deviance    Pr(>Chi)    
## 1    3                            
## 2    5  2   25.668 0.000002669 ***
## 3    9  4   10.301     0.03565 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## [2] Species-specific intensities changing with time (Model 1), with parallel changes in each road (Model 2: intercepts model), or with different temporal changes for each road (Model 3: slopes model) 
## Analysis of Deviance Table
## 
## Model 1: Nidos ~ marks * Year     Poisson
## Model 2: Nidos ~ Picada + (marks * Year)      Poisson
## Model 3: Nidos ~ Picada * (marks * Year)      Poisson
##   Npar Df Deviance Pr(>Chi)
## 1    9                     
## 2   10  1   0.7401   0.3896
## 3   18  8   7.2037   0.5148
## 
## [3] Homogeneous species-specific intensities changing with time (Model 1) or with intensity variations depending on a single, time-invariant linear relationship with the edaphic variable (Models 2), or in species-specific ways (Models 3), or with a quadratic relationship (Models 4) that changes with time (Models 5) 
## Analysis of Deviance Table
## 
## Model 1: Nidos ~ marks * Year     Poisson
## Model 2: Nidos ~ (marks * Year) + Perm    Poisson
## Model 3: Nidos ~ (marks * Year) + (marks * Perm)      Poisson
## Model 4: Nidos ~ (marks * Year) + (marks * polynom(Perm, 2))      Poisson
## Model 5: Nidos ~ marks * Year * polynom(Perm, 2)      Poisson
##   Npar Df Deviance       Pr(>Chi)    
## 1    9                               
## 2   10  1   35.517 0.000000002529 ***
## 3   12  2  149.538      < 2.2e-16 ***
## 4   15  3   35.251 0.000000107799 ***
## 5   27 12    8.218         0.7679    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Analysis of Deviance Table
## 
## Model 1: Nidos ~ marks * Year     Poisson
## Model 2: Nidos ~ (marks * Year) + Alt     Poisson
## Model 3: Nidos ~ (marks * Year) + (marks * Alt)   Poisson
## Model 4: Nidos ~ (marks * Year) + (marks * polynom(Alt, 2))   Poisson
## Model 5: Nidos ~ marks * Year * polynom(Alt, 2)   Poisson
##   Npar Df Deviance   Pr(>Chi)    
## 1    9                           
## 2   10  1    9.126    0.00252 ** 
## 3   12  2   85.693  < 2.2e-16 ***
## 4   15  3   24.767 0.00001727 ***
## 5   27 12   11.611    0.47740    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## [4] Model with species-specific quadratic effects of soil permeability (Model 1), with an additional independent quadratic effect of topography (Model 2), with both species-specific effects interacting (Model 3), and effects that change with time (Model 4)
## Analysis of Deviance Table
## 
## Model 1: Nidos ~ (marks * Year) + (marks * polynom(Perm, 2))      Poisson
## Model 2: Nidos ~ (marks * Year) + (marks * polynom(Perm, 2)) + (marks * polynom(Alt, 2))      Poisson
## Model 3: Nidos ~ (marks * Year) + (marks * polynom(Perm, 2) * polynom(Alt, 2))    Poisson
## Model 4: Nidos ~ marks * Year * polynom(Perm, 2) * polynom(Alt, 2)    Poisson
##   Npar Df Deviance      Pr(>Chi)    
## 1   15                              
## 2   21  6   43.158 0.00000010856 ***
## 3   33 12   57.428 0.00000006616 ***
## 4   81 48   34.171        0.9339    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The influence of the edaphic variables in the best model that assumes roads as replicates of the same species-specific Poisson spatial processes (Fig. 36) are very similar to those obtained by fitting every year as a different model (Fig. 35).

# predict y effectfun no funcionan con mppm sobre lpp: hay que hacerlo como combinación lineal de los coeficientes o con producto matricial armando las matrices de diseño y de coeficientes)
# predicciones en la escala de los datos (simil effectfun)
# efectos: perm (con alt = quantile(25%, median, 75%)); Year categórica
# triplicado con bind_rows para hacer una versión por Especie
temp.data <- data.frame(Alt = rep(rep(quantile(c(LaFad_Alt_TPI(0:1500, 0), Doble_Alt_TPI(0:1500, 0)), probs = c(0.25, 0.5, 0.75)), each = 50), 3), 
                        Perm = rep(rep(seq(min(c(LaFad_Perm_Smooth(0:1500, 0), Doble_Perm_Smooth(0:1500, 0))), max(c(LaFad_Perm_Smooth(0:1500, 0), Doble_Perm_Smooth(0:1500, 0))), length.out = 50), 3), 3),
                        Year = rep(c("2004", "2007", "2010"), each = 150), 
                        Lwidth = rep(rep(1:3, each = 50), 3)) %>%
  bind_rows("PI" = ., "PM" = ., "PR" = ., .id = "Sp")

coef.data <- coef(PoissonTodosBothInter)

temp.matrix <- temp.data %>%
  mutate(V01 = TRUE, V02 = Sp == "PM", V03 = Sp == "PR", V04 = Year == "2007", V05 = Year == "2010", V06 = Perm, V07 = Perm^2, V08 = Alt, V09 = Alt^2,
         V10 = V02 * V04, V11 = V03 * V04, V12 = V02 * V05, V13 = V03 * V05, V14 = V02 * V06, V15 = V03 * V06, V16 = V02 * V07, V17 = V03 * V07, V18 = V02 * V08, V19 = V03 * V08, V20 = V02 * V09, V21 = V03 * V09,
         V22 = V06 * V08, V23 = V07 * V08, V24 = V06 * V09, V25 = V07 * V09,
         V26 = V02 * V22, V27 = V03 * V22, V28 = V02 * V23, V29 = V03 * V23, V30 = V02 * V24, V31 = V03 * V24, V32 = V02 * V25, V33 = V03 * V25, .keep = "none")

temp.data$Lambda <- exp(as.matrix(temp.matrix) %*% coef.data)[,1]

grid.arrange(ggplot(subset(temp.data, Sp == "PI"), aes(x = Perm, y = Lambda, group = interaction(Year, Alt), linewidth = Lwidth, colour = as.factor(Year))) + geom_line(lineend = "round") + scale_linewidth(breaks = 1.5:5.5, range = c(0.3, 1.3)) + scale_colour_manual(values = c("#6600FF75", "#00993375", "#FF660075")) + coord_cartesian(ylim = c(0, 0.004)) + geom_vline(xintercept = c(0.7, -0.86), colour = "blue", linewidth = 0.5, linetype = 2) + labs(title = "PI", y  = NULL, x = NULL) + scale_y_continuous(labels = label_number()) + tema,
             ggplot(subset(temp.data, Sp == "PM"), aes(x = Perm, y = Lambda, group = interaction(Year, Alt), linewidth = Lwidth, colour = as.factor(Year))) + geom_line(lineend = "round") + scale_linewidth(breaks = 1.5:5.5, range = c(0.3, 1.3)) + scale_colour_manual(values = c("#6600FF75", "#00993375", "#FF660075")) + coord_cartesian(ylim = c(0, 0.0008)) + geom_vline(xintercept = c(0.7, -0.86), colour = "blue", linewidth = 0.5, linetype = 2) + labs(title = "PM", y  = NULL, x = NULL) + scale_y_continuous(labels = label_number()) + tema,
             ggplot(subset(temp.data, Sp == "PR"), aes(x = Perm, y = Lambda, group = interaction(Year, Alt), linewidth = Lwidth, colour = as.factor(Year))) + geom_line(lineend = "round") + scale_linewidth(breaks = 1.5:5.5, range = c(0.3, 1.3)) + scale_colour_manual(values = c("#6600FF75", "#00993375", "#FF660075")) + coord_cartesian(ylim = c(0, 0.0015)) + geom_vline(xintercept = c(0.7, -0.86), colour = "blue", linewidth = 0.5, linetype = 2) + labs(title = "PR", y  = NULL, x = NULL) + scale_y_continuous(labels = label_number()) + tema,
             nrow = 1, left = "Intensity (λ)", bottom = "Soil permeability (SPI)")

# efectos: alt (con perm = {25%, median, 75%}); Year categórica
temp.data <- data.frame(Perm = rep(rep(quantile(c(LaFad_Perm_Smooth(0:1500, 0), Doble_Perm_Smooth(0:1500, 0)), probs = c(0.25, 0.5, 0.75)), each = 50), 3), 
                        Alt = rep(rep(seq(min(c(LaFad_Alt_TPI(0:1500, 0), Doble_Alt_TPI(0:1500, 0))), max(c(LaFad_Alt_TPI(0:1500, 0), Doble_Alt_TPI(0:1500, 0))), length.out = 50), 3), 3),
                        Year = rep(c("2004", "2007", "2010"), each = 150), 
                        Lwidth = rep(rep(1:3, each = 50), 3))  %>%
  bind_rows("PI" = ., "PM" = ., "PR" = ., .id = "Sp")

temp.matrix <- temp.data %>%
  mutate(V01 = TRUE, V02 = Sp == "PM", V03 = Sp == "PR", V04 = Year == "2007", V05 = Year == "2010", V06 = Perm, V07 = Perm^2, V08 = Alt, V09 = Alt^2,
         V10 = V02 * V04, V11 = V03 * V04, V12 = V02 * V05, V13 = V03 * V05, V14 = V02 * V06, V15 = V03 * V06, V16 = V02 * V07, V17 = V03 * V07, V18 = V02 * V08, V19 = V03 * V08, V20 = V02 * V09, V21 = V03 * V09,
         V22 = V06 * V08, V23 = V07 * V08, V24 = V06 * V09, V25 = V07 * V09,
         V26 = V02 * V22, V27 = V03 * V22, V28 = V02 * V23, V29 = V03 * V23, V30 = V02 * V24, V31 = V03 * V24, V32 = V02 * V25, V33 = V03 * V25, .keep = "none")

temp.data$Lambda <- exp(as.matrix(temp.matrix) %*% coef.data)[,1]

grid.arrange(ggplot(subset(temp.data, Sp == "PI"), aes(x = Alt, y = Lambda, group = interaction(Year, Perm), linewidth = Lwidth, colour = as.factor(Year))) + geom_line(lineend = "round") + scale_linewidth(breaks = 1.5:5.5, range = c(0.3, 1.3)) + scale_colour_manual(values = c("#6600FF75", "#00993375", "#FF660075")) + coord_cartesian(ylim = c(0, 0.004)) + geom_vline(xintercept = 0, colour = "blue", linewidth = 0.5, linetype = 2) + labs(title = "PI", y  = NULL, x = NULL) + scale_y_continuous(labels = label_number()) + tema,
             ggplot(subset(temp.data, Sp == "PM"), aes(x = Alt, y = Lambda, group = interaction(Year, Perm), linewidth = Lwidth, colour = as.factor(Year))) + geom_line(lineend = "round") + scale_linewidth(breaks = 1.5:5.5, range = c(0.3, 1.3)) + scale_colour_manual(values = c("#6600FF75", "#00993375", "#FF660075")) + coord_cartesian(ylim = c(0, 0.0008)) + geom_vline(xintercept = 0, colour = "blue", linewidth = 0.5, linetype = 2) + labs(title = "PM", y  = NULL, x = NULL) + scale_y_continuous(labels = label_number()) + tema,
             ggplot(subset(temp.data, Sp == "PR"), aes(x = Alt, y = Lambda, group = interaction(Year, Perm), linewidth = Lwidth, colour = as.factor(Year))) + geom_line(lineend = "round") + scale_linewidth(breaks = 1.5:5.5, range = c(0.3, 1.3)) + scale_colour_manual(values = c("#6600FF75", "#00993375", "#FF660075")) + coord_cartesian(ylim = c(0, 0.0015)) + geom_vline(xintercept = 0, colour = "blue", linewidth = 0.5, linetype = 2) + labs(title = "PR", y  = NULL, x = NULL) + scale_y_continuous(labels = label_number()) + tema,
             nrow = 1, left = "Intensity (λ)", bottom = "TPI [m]")

Figure 36: Same as Fig. 35 for the best single model that assumes both road segments as replicates and includes species-specific coefficientes for Year and both edaphic characteristics as predictors.

Summary of the main models:

Model	Structure	Intensities	AIC
NULL	spp	βPI: 0.0004 βPM: 0.00019 βPR: 0.00037	6234.73
YEAR	spp×year	βPI2004: 0.0004 βPI2007: 0.00049 βPI2010: 0.00032 βPM2004: 0.00023 βPM2007: 0.00021 βPM2010: 0.00013 βPR2004: 0.00055 βPR2007: 0.00039 βPR2010: 0.00016	6210.76
ROAD	spp×year×road	LaFad: βPI2004: 0.00037 βPI2007: 0.00053 βPI2010: 0.00027 βPM2004: 0.00023 βPM2007: 0.00022 βPM2010: 0.00012 βPR2004: 0.00062 βPR2007: 0.00043 βPR2010: 0.00023 Doble: βPI2004: 0.0004 βPI2007: 0.00045 βPI2010: 0.00037 βPM2004: 0.00023 βPM2007: 0.0002 βPM2010: 0.00015 βPR2004: 0.00048 βPR2007: 0.00035 βPR2010: 0.00008	6220.82
PERM2	sppxyear+spp×SPI2	βPI2004: 0.00023+0.09732SPI+0.53692SPI² βPI2007: 0.00028+0.09732SPI+0.53692SPI² βPI2010: 0.00018+0.09732SPI+0.53692SPI² βPM2004: 0.00019+10.25538SPI+0.33766SPI² βPM2007: 0.00017+10.25538SPI+0.33766SPI² βPM2010: 0.00011+10.25538SPI+0.33766SPI² βPR2004: 0.00074+0.39268SPI+0.19794SPI² βPR2007: 0.00053+0.39268SPI+0.19794SPI² βPR2010: 0.00021+0.39268SPI+0.19794SPI²	6002.46
TOPO2	sppxyear+spp×TPI2	βPI2004: 0.00037+0.04066TPI+0.11649TPI² βPI2007: 0.00045+0.04066TPI+0.11649TPI² βPI2010: 0.00029+0.04066TPI+0.11649TPI² βPM2004: 0.00018+5531.81387TPI+0TPI² βPM2007: 0.00016+5531.81387TPI+0TPI² βPM2010: 0.0001+5531.81387TPI+0TPI² βPR2004: 0.00065+0.13026TPI+0.00049TPI² βPR2007: 0.00046+0.13026TPI+0.00049TPI² βPR2010: 0.00019+0.13026TPI+0.00049TPI²	6103.18
INDEP	spp×year+spp×(SPI2+TPI2)	βPI2004: 0.00025+0.1167SPI+0.62096SPI²+0.72595TPI+0.57932TPI² βPI2007: 0.00031+0.1167SPI+0.62096SPI²+0.72595TPI+0.57932TPI² βPI2010: 0.0002+0.1167SPI+0.62096SPI²+0.72595TPI+0.57932TPI² βPM2004: 0.00016+4.33337SPI+0.60064SPI²+567.2474TPI+0.00003TPI² βPM2007: 0.00014+4.33337SPI+0.60064SPI²+567.2474TPI+0.00003TPI² βPM2010: 0.00009+4.33337SPI+0.60064SPI²+567.2474TPI+0.00003TPI² βPR2004: 0.00079+0.54089SPI+0.22232SPI²+0.14353TPI+0.02482TPI² βPR2007: 0.00056+0.54089SPI+0.22232SPI²+0.14353TPI+0.02482TPI² βPR2010: 0.00023+0.54089SPI+0.22232SPI²+0.14353TPI+0.02482TPI²	5971.3
INTER	spp×year+spp×SPI2×TPI2	see Fig. 36	5937.87
FULL	spp×year×SPI2×TPI2	(too long)	5999.7

Soil texture in monospecific areas

Most soil samples collected near Pogonomyrmex colonies in the dirt roads were dominated by sand particles (classifiable as Sandy, Loamy Sand or Sandy Loam; Fig. 37). Samples with smaller average particle size (more silt and/or clay) appear always associated with clusters (road sectors) of PI colonies, while samples near PM colonies seem consistently the sandiest sampled.

grid.arrange(
  # A. As stacked bars
  ggplot(subset(Textura_long, Especie != "Random"), aes(x = Muestra, y = Percentage, fill = factor(Tamaño, levels = c("Arena", "Limo", "Arcilla")))) +
  geom_col() + 
  facet_wrap(~ Especie, scales = "free_x", strip.position = "bottom") + 
  labs(x = "Soil samples", y = "Soil composition (%)") +
  scale_fill_manual(name = "Particle size", labels = c("Sand", "Silt", "Clay"), values = c("lightgoldenrod1", "seashell3", "sienna2")) +
  tema + theme(legend.position = "right", strip.background = element_blank(), strip.placement = "outside", strip.text = element_text(size = rel(1.1))),
  # B. Ternary plot
  ggplot(subset(Textura, Especie != "Random"), aes(x = Arena, z = Limo, y = Arcilla, colour = Especie)) +
  coord_tern() +
  geom_point(size = 0.5) +
  geom_point(data = subset(TexturaXCluster, Especie != "Random"), size = 2, shape = 21, stroke = 0.75) +
  scale_colour_manual(values = c("darkgreen", "darkblue")) +
  labs(x = "Sand", y = "Clay", z = "Silt", xarrow = "Sand (%)", yarrow = "Clay (%)", zarrow = "Silt (%)") + 
  tema +
  theme_showarrows() +
  theme_clockwise(),
nrow = 1, widths = c(1.5, 1))

Figure 37: Soil composition by particle size (sand, silt and clay) of superficial and subsuperficial (0–30 cm) samples associated to Pogonomyrmex inermis (green) and P. mendozanus (blue) active colonies in dirt roads. Left: Percentual composition of individual soil samples (one per colony). Right: Ternary plot, with individual samples as dots and open circles for average percentual values for each spatial cluster of (3) colonies sampled.

Mixed models with Beta distribution to account for the response variable (proportions of each particle size) and nested sampling design (3 colonies sampled in 5 spatial clusters per species) showed both mean proportion and intersample dispersion (higher Beta dispersion parameter corresponds to lower variance) associated with the Species (fixed factor) of the colony sampled.

For the proportion of sand in the soil samples, deviance tests (LR or AIC) suggest the best model includes Species modifying the mean proportion (and also the dispersion parameter of the Beta distribution for a better fit: Fig. 38):

ArenaSpp <- glmmTMB(Arena/100 ~ Especie + (1 | Cluster), data = subset(Textura, Especie != "Random"), family = beta_family(link = "logit"))
ArenaSppHet <- glmmTMB(Arena/100 ~ Especie + (1 | Cluster), dispformula = ~ Especie, data = subset(Textura, Especie != "Random"), family = beta_family(link = "logit"))

#DHARMa objects: simulated residuals
ArenaSppResiduals <- simulateResiduals(ArenaSpp, n = 1000)
ArenaSppHetResiduals <- simulateResiduals(ArenaSppHet, n = 1000)

layout(matrix(1:4, 1, 4, byrow = TRUE))
plotQQunif(ArenaSppResiduals, main = "Sand ~ Spp (dispersion = ~1)") #not too bad, although not uniform
plotResiduals(ArenaSppResiduals, form = factor(c(rep("PM", 15), rep("PI", 15))))
plotQQunif(ArenaSppHetResiduals, main = "Sand ~ Spp (dispersion = ~Spp)") #better!
plotResiduals(ArenaSppHetResiduals, form = factor(c(rep("PM", 15), rep("PI", 15)))) #much better!

anova(ArenaSpp, ArenaSppHet) # heterogeneous is (also) much better fit

## Data: subset(Textura, Especie != "Random")
## Models:
## ArenaSpp: Arena/100 ~ Especie + (1 | Cluster), zi=~0, disp=~1
## ArenaSppHet: Arena/100 ~ Especie + (1 | Cluster), zi=~0, disp=~Especie
##             Df     AIC     BIC logLik deviance  Chisq Chi Df   Pr(>Chisq)    
## ArenaSpp     4 -17.223 -11.619 12.612  -25.223                               
## ArenaSppHet  5 -41.222 -34.216 25.611  -51.222 25.999      1 0.0000003416 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#also tested with different random parameters for each level of Especie [~ Especie + (0 + dummy(Especie, "PI") | Cluster) + (0 + dummy(Especie, "PM") | Cluster)] , but it was not better than (only) accounting for heterogeneous dispersions

drop1(ArenaSppHet, test = "Chisq") #test considering heterogeneous dispersions in the null, which make not much sense without Especie as fixed factor in the model…

## Single term deletions
## 
## Model:
## Arena/100 ~ Especie + (1 | Cluster)
##         Df     AIC    LRT  Pr(>Chi)    
## <none>     -41.222                     
## Especie  1 -31.601 11.622 0.0006518 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ArenaNull <- glmmTMB(Arena/100 ~ 1 + (1 | Cluster), data = subset(Textura, Especie != "Random"), family = beta_family(link = "logit"))

anova(ArenaNull, ArenaSppHet)

## Data: subset(Textura, Especie != "Random")
## Models:
## ArenaNull: Arena/100 ~ 1 + (1 | Cluster), zi=~0, disp=~1
## ArenaSppHet: Arena/100 ~ Especie + (1 | Cluster), zi=~0, disp=~Especie
##             Df     AIC     BIC logLik deviance  Chisq Chi Df   Pr(>Chisq)    
## ArenaNull    3 -14.769 -10.565 10.384  -20.769                               
## ArenaSppHet  5 -41.222 -34.216 25.611  -51.222 30.453      2 0.0000002439 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#best model. Higher dispersion parameter means lower variance: V = (mean - mean^2) / (disp + 1)
ArenaSppHet

## Formula:          Arena/100 ~ Especie + (1 | Cluster)
## Dispersion:                 ~Especie
## Data: subset(Textura, Especie != "Random")
##       AIC       BIC    logLik  df.resid 
## -41.22235 -34.21636  25.61117        25 
## Random-effects (co)variances:
## 
## Conditional model:
##  Groups  Name        Std.Dev.
##  Cluster (Intercept) 0.1584  
## 
## Number of obs: 30 / Conditional model: Cluster, 10
## 
## Fixed Effects:
## 
## Conditional model:
## (Intercept)    EspeciePM  
##      0.1261       1.1495  
## 
## Dispersion model:
## (Intercept)    EspeciePM  
##       1.141        3.462

ArenaSppHet_pred <- predict(ArenaSppHet, newdata = data.frame(Especie = c("PI", "PM")), re.form = NA, type = "response")

Figure 38: Tests of simulated residuals of GLMM models (logit link, assuming Beta distribution of error) fit to the proportion of sand in soil samples associated with Pogonomyrmex inermis and P. mendozanus active colonies in dirt roads, with Species as fixed factor and a random effect on the intercept of the spatial clusters of colonies. Left: Model with common single parameter for dispersion; right: Model with heterogeneous dispersion depending on Species.

Under the selected model, the mean proportion of sand in the soil samples was 47.1% higher near PM colonies (78.2%) than around PI colonies (53.1%).

Model fit was similar, with even more extreme results but in the opposite direction, for clay content (also in a model with heterogeneous dispersion: Fig 39):

ArcillaSpp <- glmmTMB(Arcilla/100 ~ Especie + (1 | Cluster), data = subset(Textura, Especie != "Random"), family = beta_family(link = "logit"))
ArcillaSppHet <- glmmTMB(Arcilla/100 ~ Especie + (1 | Cluster), dispformula = ~ Especie, data = subset(Textura, Especie != "Random"), family = beta_family(link = "logit"))

#DHARMa objects: simulated residuals
ArcillaSppResiduals <- simulateResiduals(ArcillaSpp, n = 1000)
ArcillaSppHetResiduals <- simulateResiduals(ArcillaSppHet, n = 1000)

layout(matrix(1:4, 1, 4, byrow = TRUE))
plotQQunif(ArcillaSppResiduals, main = "Clay ~ Spp (dispersion = ~1)") #not too bad, although not uniform
plotResiduals(ArcillaSppResiduals, form = factor(c(rep("PM", 15), rep("PI", 15))))
plotQQunif(ArcillaSppHetResiduals, main = "Clay ~ Spp (dispersion = ~Spp)") #better!
plotResiduals(ArcillaSppHetResiduals, form = factor(c(rep("PM", 15), rep("PI", 15)))) #much better!

anova(ArcillaSpp, ArcillaSppHet) # heterogeneous is (also) much better fit

## Data: subset(Textura, Especie != "Random")
## Models:
## ArcillaSpp: Arcilla/100 ~ Especie + (1 | Cluster), zi=~0, disp=~1
## ArcillaSppHet: Arcilla/100 ~ Especie + (1 | Cluster), zi=~0, disp=~Especie
##               Df      AIC     BIC logLik deviance  Chisq Chi Df    Pr(>Chisq)
## ArcillaSpp     4  -76.218 -70.614 42.109  -84.218                            
## ArcillaSppHet  5 -106.573 -99.567 58.286 -116.573 32.355      1 0.00000001285
##                  
## ArcillaSpp       
## ArcillaSppHet ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#also tested with different random parameters for each level of Especie [~ Especie + (0 + dummy(Especie, "PI") | Cluster) + (0 + dummy(Especie, "PM") | Cluster)] , but it was not better than (only) accounting for heterogeneous dispersions

drop1(ArcillaSppHet, test = "Chisq") #test considering heterogeneous dispersions in the null, which make not much sense without Especie as fixed factor in the model…

## Single term deletions
## 
## Model:
## Arcilla/100 ~ Especie + (1 | Cluster)
##         Df      AIC    LRT   Pr(>Chi)    
## <none>     -106.573                      
## Especie  1  -91.609 16.964 0.00003809 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ArcillaNull <- glmmTMB(Arcilla/100 ~ 1 + (1 | Cluster), data = subset(Textura, Especie != "Random"), family = beta_family(link = "logit"))

anova(ArcillaNull, ArcillaSppHet)

## Data: subset(Textura, Especie != "Random")
## Models:
## ArcillaNull: Arcilla/100 ~ 1 + (1 | Cluster), zi=~0, disp=~1
## ArcillaSppHet: Arcilla/100 ~ Especie + (1 | Cluster), zi=~0, disp=~Especie
##               Df      AIC     BIC logLik deviance  Chisq Chi Df    Pr(>Chisq)
## ArcillaNull    3  -73.038 -68.834 39.519  -79.038                            
## ArcillaSppHet  5 -106.573 -99.567 58.286 -116.573 37.535      2 0.00000000707
##                  
## ArcillaNull      
## ArcillaSppHet ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#best model. Higher dispersion parameter means lower variance: V = (mean - mean^2) / (disp + 1)
ArcillaSppHet

## Formula:          Arcilla/100 ~ Especie + (1 | Cluster)
## Dispersion:                   ~Especie
## Data: subset(Textura, Especie != "Random")
##        AIC        BIC     logLik   df.resid 
## -106.57296  -99.56698   58.28648         25 
## Random-effects (co)variances:
## 
## Conditional model:
##  Groups  Name        Std.Dev.  
##  Cluster (Intercept) 0.00003391
## 
## Number of obs: 30 / Conditional model: Cluster, 10
## 
## Fixed Effects:
## 
## Conditional model:
## (Intercept)    EspeciePM  
##     -1.5361      -0.8759  
## 
## Dispersion model:
## (Intercept)    EspeciePM  
##       2.515        3.622

ArcillaSppHet_pred <- predict(ArcillaSppHet, newdata = data.frame(Especie = c("PI", "PM")), re.form = NA, type = "response")

Figure 39: Same as Fig. 38 for proportion of clay in the soil samples.

Under the best fitted model, the mean proportion of clay in the soil samples was 115.3% higher around PI colonies (17.7%) than close to PM colonies (8.2%).

---

Conclusions

Colonies of the harvester ants in the Pogonomyrmex genus appeared rather homogeneously along the dirt roads crossing the algarrobal of Ñacuñán, but each of the three species showed intraspecific patterns of aggregation, and interespecific pattern of segregation between pairs of species, which resulted in 3–7 patches per species along each 1.5 km segment censused. These patched distributions probably result from species-specific responses to variable environmental conditions along the roads.

The two environmental variables that we estimated, soil permeability and relative altitude (which are partially correlated), both resulted relevant to explain the observed variations in active colony densities along each of the roads (strictly, in the intensity of the Poisson point pattern process that influence their location on both linear paths). This was strikingly clear for two of the three species (PI y PM) that showed contrasting associations and, as a consequence, resulted spatially segregated: PI selected more impermeable (clayish) soils that tend to accumulate in relative depressions, and PM more permeable (sandy) soils which predominate in relative highs. PR is more broadly distributed, asociated with the intermediate conditions typical of the algarrobal.
The habitat associations and the resulting spatial patterns of each Pogonomyrmex species along the dirt roads persisted in time in spite of strong declines in the density of active nests of the three species (particularly PR) between 2004 and 2010.