Resumen

Analizamos cómo varían la abundancia y arreglo espacial de colonias activas de Pheidole bergi con la heterogeneidad del hábitat a distintas escalas espaciales y con los cambios embientales de largo plazo en el desierto del Monte central (Reserva de la Biósfera de Ñacuñán y alrededores, Provincia de Mendoza, Argentina). La densidad de colonias en grillas permanentes dentro de un algarrobal protegido del pastoreo se redujeron notoriamente entre 2001 y 2019. La disposición espacial de colonias a escalas pequeñas no resultó aleatoria sino influida por la vegetación leñosa y la densidad de nidos (2001 > 2019). La mayor probabilidad de encontrar una colonia está en los bordes entre parches de vegetación leñosa y parches descubiertos, y baja notoriamente hacia los centros de los parches leñosos. Luego, detectamos interacciones entre nidos cercanos (hasta un radio de ~6 m) que, al empeorar las condiciones ambientales (sequía) y bajar la densidad de nidos (2019) resulta en menor agregamiento entre colonias cercanas, llegando incluso a patrones espaciales de inhibición o repulsión. Este raleo no aleatorio podría deberse a competencia directa entre colonias vecinas.

Introducción

Métodos

knitr::include_graphics("Fig_grillas.png")
Arriba izquierda: *Pheidole bergi* (obrera). Arriba derecha: entrada típica de colonia de *P.bergi*. Abajo izquierda: buscando y mapeando colonias en los parches arbustivos de parcelas permanentes en el algarrobal protegido de la Reserva de Ñacuñán (Mendoza, Argentina). Abajo derecha: heterogeneidad característica de los algarrobales del Monte central y muchos otros ambientes áridos y semiáridos: parches lindantes con y sin cobertura leñosa.

Figura 1: Arriba izquierda: Pheidole bergi (obrera). Arriba derecha: entrada típica de colonia de P.bergi. Abajo izquierda: buscando y mapeando colonias en los parches arbustivos de parcelas permanentes en el algarrobal protegido de la Reserva de Ñacuñán (Mendoza, Argentina). Abajo derecha: heterogeneidad característica de los algarrobales del Monte central y muchos otros ambientes áridos y semiáridos: parches lindantes con y sin cobertura leñosa.

Análisis y Resultados

Mapas de nidos y de vegetación leñosa

library(spatstat) # documented in "Spatial Point Patterns: Methodology and Applications with R" by Baddeley et al (2016)
library(png) # para importar las imágenes con readPNG
library(scales)

# NOTA: debido a los resultados obtenidos (que no difieren de manera relevante), en 2021 eliminamos toda la exploración de análisis hechos con la vegetación rasterizada a escala 100×100 y nos quedamos con la que tiene la misma escala que como fueron mapeados los nidos (= 50×50)

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

  # summary(J2001_ants_50) confirma a los objetos im creados como imágenes de 50×50 pixeles de 1×1 m dentro del rectángulo [0.5, 50.5]x[0.5, 50.5] m

  # vegetación leñosa a objetos spatstat::im (vegetación leñosa: pixel = 0 | sd: pixel = 1)
J2001_veg_50 <- im(transmat(readPNG("J2001_veg_50x50.png"), from = "European", to = "spatstat"), unitname = "m")
J2019_veg_50 <- im(transmat(readPNG("J2019_veg_50x50.png"), from = "European", to = "spatstat"), unitname = "m")
V2019_veg_50 <- im(transmat(readPNG("V2019_veg_50x50.png"), from = "European", to = "spatstat"), unitname = "m")

#escala 100×100 eliminada (2021)
# Imágenes raster (png) binarias en grillas 100x100 (pixel = 0.25 m^2^ => 10000 pixeles)
  # vegetación leñosa a objetos spatstat::im (vegetación: pixel = 0 | sd: pixel = 1)
#J2001_veg_100 <- im(transmat(readPNG("J2001_veg_100x100.png"), from = "European", to = "spatstat"), unitname = c("m", "m", 0.5))
#J2001_ants_100 <- im(transmat(readPNG("J2001_nidos_100x100.png"), from = "European", to = "spatstat"), unitname = c("m", "m", 0.5))
#J2019_veg_100 <- im(transmat(readPNG("J2019_veg_100x100.png"), from = "European", to = "spatstat"), unitname = c("m", "m", 0.5))
#J2019_ants_100 <- im(transmat(readPNG("J2019_nidos_100x100.png"), from = "European", to = "spatstat"), unitname = c("m", "m", 0.5))
#V2019_veg_100 <- im(transmat(readPNG("V2019_veg_100x100.png"), from = "European", to = "spatstat"), unitname = c("m", "m", 0.5))
#V2019_ants_100 <- im(transmat(readPNG("V2019_nidos_100x100.png"), from = "European", to = "spatstat"), unitname = c("m", "m", 0.5))
# no hay vegetación medida en V2001
#V2001_ants_100 <- im(transmat(readPNG("V2001_nidos_100x100.png"), from = "European", to = "spatstat"), unitname = c("m", "m", 0.5))

# Crea máscaras o ventanas de análisis de los píxeles CON nidos (spatstat::owin) (NA = sin nidos)
J2001_antmask_50 <- J2001_ants_50
is.na(J2001_antmask_50) <- eval.im(J2001_ants_50 == 1)
J2001_antmask_50 <- as.owin(J2001_antmask_50)

J2019_antmask_50 <- J2019_ants_50
is.na(J2019_antmask_50) <- eval.im(J2019_ants_50 == 1)
J2019_antmask_50 <- as.owin(J2019_antmask_50)

V2001_antmask_50 <- V2001_ants_50
is.na(V2001_antmask_50) <- eval.im(V2001_ants_50 == 1)
V2001_antmask_50 <- as.owin(V2001_antmask_50)

V2019_antmask_50 <- V2019_ants_50
is.na(V2019_antmask_50) <- eval.im(V2019_ants_50 == 1)
V2019_antmask_50 <- as.owin(V2019_antmask_50)

All_antmask_50 <- list("J 2001" = J2001_antmask_50, "J 2019" = J2019_antmask_50, "V 2001" = V2001_antmask_50, "V 2019" = V2019_antmask_50)

# y, lo mismo, pero como patrones de puntos (ubicando a los nidos en los centros de sus píxeles)
J2001_ppp <- pixelcentres(J2001_antmask_50)
J2019_ppp <- pixelcentres(J2019_antmask_50)
V2001_ppp <- pixelcentres(V2001_antmask_50)
V2019_ppp <- pixelcentres(V2019_antmask_50)

All_ppp <- list("J 2001" = J2001_ppp, "J 2019" = J2019_ppp, "V 2001" = V2001_ppp, "V 2019" = V2019_ppp)

# Crea máscaras o ventanas de análisis de los píxeles CON vegetación leñosa (spatstat::owin) (sd: NA)
J2001_vegmask_50 <- J2001_veg_50
is.na(J2001_vegmask_50) <- eval.im(J2001_veg_50 == 1)
J2001_vegmask_50 <- as.owin(J2001_vegmask_50)

J2019_vegmask_50 <- J2019_veg_50
is.na(J2019_vegmask_50) <- eval.im(J2019_veg_50 == 1)
J2019_vegmask_50 <- as.owin(J2019_vegmask_50)

V2019_vegmask_50 <- V2019_veg_50
is.na(V2019_vegmask_50) <- eval.im(V2019_veg_50 == 1)
V2019_vegmask_50 <- as.owin(V2019_vegmask_50)

#escala 100×100 eliminada (2021)
#J2001_vegmask_100 <- J2001_veg_100
#is.na(J2001_vegmask_100) <- eval.im(J2001_veg_100 == 1)
#J2001_vegmask_100 <- as.owin(J2001_vegmask_100)
#J2019_vegmask_100 <- J2019_veg_100
#is.na(J2019_vegmask_100) <- eval.im(J2019_veg_100 == 1)
#J2019_vegmask_100 <- as.owin(J2019_vegmask_100)
#V2019_vegmask_100 <- V2019_veg_100
#is.na(V2019_vegmask_100) <- eval.im(V2019_veg_100 == 1)
#V2019_vegmask_100 <- as.owin(V2019_vegmask_100)

# Si quisiera hacer máscaras SIN vegetación (spatstat::owin) puedo hacerlas como complemento
# J2001_sdmask_50 <- complement.owin(J2001_vegmask_50), pero no hace falta la conversión solo para estimar distancias a bordes de máscara: usar distmap(X, invert = TRUE) [de hecho invert usa internamente complement.owin])
# [old: ¡OJO PORQUE HAY DIFERENCIAS EN LAS DISTANCIAS EN LOS BORDES DE LAS GRILLAS (¿al invertir toma al borde como dato de la otra categoría?)! cf. plot(distmap(J2001_vegmask_50)) vs. plot(distmap(complement.owin(J2001_vegmask_50), invert = TRUE))); mejor no usar la opción invert = TRUE y hacer las máscaras de la otra categoría…] -> FIXED in spatstat.geom 2.3-2.011! 
J2001_sdmask_50 <- J2001_veg_50
is.na(J2001_sdmask_50) <- eval.im(J2001_veg_50 == 0)
J2001_sdmask_50 <- as.owin(J2001_sdmask_50)

J2019_sdmask_50 <- J2019_veg_50
is.na(J2019_sdmask_50) <- eval.im(J2019_veg_50 == 0)
J2019_sdmask_50 <- as.owin(J2019_sdmask_50)

V2019_sdmask_50 <- V2019_veg_50
is.na(V2019_sdmask_50) <- eval.im(V2019_veg_50 == 0)
V2019_sdmask_50 <- as.owin(V2019_sdmask_50)

#escala 100×100 eliminada (2021)
#genera vector con valores de índices seleccionables en cualquier vector con valores correspondientes a matriz 100×100 para poder usar en modelos nulos que elijan píxeles seleccionables en matriz 50×50 (grano de medición de nidos). Si en 50×50 son todos, en 100×100 son los impares de las filas impares
#Selectable_100x100 <- integer()
#for (y in seq.int(0, 98, 2)) {Selectable_100x100 = c(Selectable_100x100, seq.int(1, 99, 2) + y * 100)}

La abundancia de colonias se redujo consistentemente en las dos parcelas (V > J) entre 2001 y 2019; las ubicaciones de los píxeles (1 m2) ocupados por nidos de Pheidole bergi en cada grilla y fecha de muestreo se presenta en la Figura 2.

# el n en el título está calculado con el objeto, así que sirve para chequear
# también se podrían plotear las versiones de máscara binaria: J2001_antmask_50, etc.
par(mfrow = c(2, 2), mar = c(1.5, 1, 1.5, 1))
for(i in names(All_ants_50)) {
  plot(All_ants_50[[i]], col = c("black", "white"), ribbon = F, main = paste(i, " (n = ", sum(!All_ants_50[[i]]$v), ")", sep = ""))
}
Ubicación de los pixeles (1 m^2^, en negro) con nidos activos de *P. bergi* en las dos parcelas permanentes de 50x50 m^2^ (J y V) en el algarrobal protegido de la Reserva de la Biosfera de Ñacuñan (Mendoza, Argentina), en el desierto del Monte central, en febrero de 2001 y febrero de 2019.

Figura 2: Ubicación de los pixeles (1 m2, en negro) con nidos activos de P. bergi en las dos parcelas permanentes de 50x50 m2 (J y V) en el algarrobal protegido de la Reserva de la Biosfera de Ñacuñan (Mendoza, Argentina), en el desierto del Monte central, en febrero de 2001 y febrero de 2019. 

# las coordenadas de píxeles con nidos (por ejemplo para hacer patrones de puntos `ppp`) se pueden obtener con: which(!J2001_ants_50$v, arr.ind = TRUE)
# o directamente un objeto ppp desde máscara binaria con: pixelcentres(J2001_antmask_50)

Los patrones espaciales de los nidos activos de cada parcela en cada momento de muestreo, evaluados mediante las funciones de autocorrelación espacial (L de Ripley, Fig. 3: curvas observadas por encima de los límites de confianza de un proceso aleatorio indican agregamiento; por debajo, sobredispersión o uniformidad) y los diagramas de Fry (Fig. 4, que muestran que el supuesto de isotropismo es razonable), sugieren que su disposición no es completamente aleatoria: aparece más agregada que lo esperado por azar durante 2001 (en especial en J), con densidades de nidos más altas, y más aleatorias a sobredispersas en 2019 (en especial en V), con menos de la mitad de nidos activos.

par(mfcol = c(2, 4), mar = c(2, 2, 1.5, 1))
# problema: spatstat::envelope genera las simulaciones con puntos en lugares no posibles para nuestro muestreo (coordenadas decimales => puntos con "variación intrapixel")
# en cambio, el patrón de puntos en cada simulación (simulate =) se genera con permutaciones de los T/F de la máscara con nidos [i.e., equivale a reubicar los nidos en coordenadas posibles de una ventana igual a la nuestra; el n queda fijo => true MonteCarlo against CSR: Baddeley et al 2014,2017]

for(i in 1:length(All_ppp)) {
  temp <- envelope(All_ppp[[i]], Lest, nsim = 399, nrank = 10, simulate = expression(pixelcentres(owin(c(0.5, 50.5), c(0.5, 50.5), mask = matrix(sample(All_antmask_50[[i]]$m), 50, 50), unitname = "m"))), funargs = list(r = 0:12, rmax = 12), global = FALSE, verbose = FALSE, savefuns = TRUE, savepatterns = TRUE)
  plot(temp, . - mmean ~ r, legend = FALSE, main = paste(names(All_ppp[i]), ": L(r)"), cex.axis = 0.8)
  text(x = 12, y = 0.5, paste("p = ", round(mad.test(temp, use.theory = FALSE)$p.value, 3)), adj = c(1, 1))
  
  temp2 <- envelope(All_ppp[[i]], pcf, nsim = 399, nrank = 10, r = 0:12, rmax = 12, verbose = FALSE, simulate = temp, savefuns = TRUE)
  plot(temp2, . - mmean + 1 ~ r, legend = FALSE, main = paste(names(All_ppp[i]), ": g(r)"), cex.axis = 0.8)
  text(x = 12, y = 1.55, paste("p = ", round(mad.test(temp2, use.theory = FALSE)$p.value, 3)), adj = c(1, 1))
}
Autocorrelación espacial de las colonias en cada parcela y año según la función L de Ripley (= función K centrada y estandarizada) y la función de correlación de pares (g) en función de la distancia r entre colonias (en metros). Líneas: valores observados; sombreado: intervalo crítico del 95% asumiendo procesos Poisson homogéneos; el p-valor corresponde a una prueba global de Diferencia Máxima Absoluta poniendo a prueba que el proceso es CSR (i.e. Poisson homogéneo, sin autocorrelación espacial).

Figura 3: Autocorrelación espacial de las colonias en cada parcela y año según la función L de Ripley (= función K centrada y estandarizada) y la función de correlación de pares (g) en función de la distancia r entre colonias (en metros). Líneas: valores observados; sombreado: intervalo crítico del 95% asumiendo procesos Poisson homogéneos; el p-valor corresponde a una prueba global de Diferencia Máxima Absoluta poniendo a prueba que el proceso es CSR (i.e. Poisson homogéneo, sin autocorrelación espacial). 

par(mfcol = c(1, 4), mar = c(2, 2, 1.5, 1))
#Fry plots
for(i in names(All_ppp)) {
  fryplot(All_ppp[[i]], width = 50, axes = T, cex = .5, main = i)
}
Diagramas de Fry de los patrones de puntos dados por los píxeles con nidos de *P. bergi* en cada parcela y año. Estos diagramas son útiles para detectar anisotropía; un espacio vacío cerca del origen sugiere regularidad (inhibición entre puntos) y su forma asimétrica sugiere anisotropia. Una zona densa en el origen sugiere agregamiento y su escala espacial asociada. También sirven para detectar periodicidad o redondeo de las coordenadas espaciales (Baddeley et al. 2016).

Figura 4: Diagramas de Fry de los patrones de puntos dados por los píxeles con nidos de P. bergi en cada parcela y año. Estos diagramas son útiles para detectar anisotropía; un espacio vacío cerca del origen sugiere regularidad (inhibición entre puntos) y su forma asimétrica sugiere anisotropia. Una zona densa en el origen sugiere agregamiento y su escala espacial asociada. También sirven para detectar periodicidad o redondeo de las coordenadas espaciales (Baddeley et al. 2016).

Entre los potenciales causantes de esos patrones no aleatorios, uno de los más evidentes es la disposición de la vegetación leñosa (típico determinante de la heterogeneidad en parches en los ambientes áridos). Los mapas de vegetación leñosa alta (>1 m), rasterizados con el mismo grano que el mapa de nidos (1 m2) [a partir de dibujos vectoriales por JLdC tomados a campo con referencias físicas cada 5 m, y con resolución no determinada pero más detallada que 1 m], con los nidos superpuestos, se muestran en la Figura 5. Nótese la aparente ausencia de nidos asociados a píxeles en el interior de los parches más grandes de vegetación leñosa.

par(mfrow = c(2, 2), mar = c(1.5, 1, 1.5, 1))

plot(J2001_veg_50, col = c("#006400CC", "white"), ribbon = F, main = "J 2001")
plot(J2001_ants_50, col = c(rgb(0,0,0,.9), rgb(1,1,1,0)), ribbon = F, add = TRUE)

plot(J2019_veg_50, col = c("#006400CC", "white"), ribbon = F, main = "J 2019")
plot(J2019_ants_50, col = c(rgb(0,0,0,.9), rgb(1,1,1,0)), ribbon = F, add = TRUE)

plot(V2019_veg_50, col = c("#006400CC", "white"), ribbon = F, main = "nidos: V 2001, veg: V 2019")
plot(V2001_ants_50, col = c(rgb(0,0,0,.9), rgb(1,1,1,0)), ribbon = F, add = TRUE)

plot(V2019_veg_50, col = c("#006400CC", "white"), ribbon = F, main = "V 2019")
plot(V2019_ants_50, col = c(rgb(0,0,0,.9), rgb(1,1,1,0)), ribbon = F, add = TRUE)
Ubicación de los píxeles con nidos de *P. bergi* (negro) en cada parcela y año y de los parches de vegetación leñosa alta (> 1 m altura, en verde).

Figura 5: Ubicación de los píxeles con nidos de P. bergi (negro) en cada parcela y año y de los parches de vegetación leñosa alta (> 1 m altura, en verde). 

#escala 100×100 eliminada (2021)
#plot(J2001_veg_100, col = c("dark green", "white"), ribbon = F, main = "J 2001 (100×100)")
#plot(J2001_ants_100, col = c(rgb(0,0,0,.8), rgb(1,1,1,0)), ribbon = F, add = TRUE)
#plot(J2019_veg_100, col = c("dark green", "white"), ribbon = F, main = "J 2019 (100×100)")
#plot(J2019_ants_100, col = c(rgb(0,0,0,.8), rgb(1,1,1,0)), ribbon = F, add = TRUE)
#plot(V2019_veg_100, col = c("dark green", "white"), ribbon = F, main = "V 2019 (100×100)")
#plot(V2019_ants_100, col = c(rgb(0,0,0,.8), rgb(1,1,1,0)), ribbon = F, add = TRUE)

Distancias a bordes de parches (transición vegetación leñosa–sd) vs. modelos nulos

A partir de los mapas de vegetación leñosa se pueden estimar las distancias 1 desde (el centro de) cada píxel (con o sin nido) hasta los bordes de la vegetación 2. Si asignamos valores positivos de distancia a los pixeles sin vegetación y negativos a aquellos bajo arbustos/árboles, se obtiene un conjunto de valores discretos de distancia. Por ejemplo, para J2001 los valores posibles son:

# distancias a borde veg-sd de cada pixel
J2001_vegdist_50 <- distmap(J2001_vegmask_50) - distmap(J2001_sdmask_50)
J2019_vegdist_50 <- distmap(J2019_vegmask_50) - distmap(J2019_sdmask_50)
sort(unique(as.vector(J2001_vegdist_50$v)))
##  [1] -4.828427 -4.414214 -4.242641 -4.000000 -3.828427 -3.414214 -3.000000
##  [8] -2.828427 -2.414214 -2.000000 -1.414214 -1.000000  1.000000  1.414214
## [15]  2.000000  2.414214  2.828427  3.000000
# Los valores eran algo mayores a lo esperado hasta spatstat 2.3-3; corregido en spatstat.geom a partir de 2.3-2.004 (ver Tests-of-distmap-values.Rmd)

Con esas distancias se puede evaluar si los píxeles que tienen nidos activos de Pheidole bergi resultan una muestra sesgada de los píxeles disponibles en función de sus distancias hasta el borde de la vegetación leñosa al evaluarlo bajo distintos modelos nulos de procesos espaciales tipo Poisson, es decir, en los que se supone que los puntos son independientes entre sí (CSR o Poisson homogéneo, o Poisson heterogéneo con intensidad dependiente de covariables)

# J2001
# dataframe con los estadísticos derivados de las permutaciones (media y desvío) bajo cada modelo nulo, incluyendo al valor observado en la primera fila
# inicializando el dataframe con vectores vacíos
J2001_vegdist_50_null <- data.frame(CSR_mean = numeric(length = 2000), 
                                    CSR_sd = numeric(length = 2000), 
                                    CSR2_mean = numeric(length = 2000), 
                                    CSR2_sd = numeric(length = 2000), 
                                    HET_mean = numeric(length = 2000), 
                                    HET_sd = numeric(length = 2000), 
                                    HET2_mean = numeric(length = 2000), 
                                    HET2_sd = numeric(length = 2000))

#escala 100×100 eliminada (2021)
#J2001_vegdist_100_null <- data.frame(CSR_mean = numeric(length = 2000))

# código nuevo para obtener los dos parámetros evaluados (mean, sd) a partir de un mismo procedimiento de muestras aleatorias
# 50×50 modelo nulo CSR estricto: distribución de medias (y de sd) de c/u de 2000 muestras aleatorias independientes sin reemplazo de tamaño N (píxeles con nidos observados) del mapa de distancias
# equivalente a ppm(J2001_ppp ~ 1), con intensidad λ constante = nidos/píxeles [Baddeley et al: 9.3.1]
J2001_vegdist_50_null[c("CSR_mean", "CSR_sd")] <- t(replicate(2000, 
    do.call(function(x) c(mean(x), sd(x)), list(
    sample(J2001_vegdist_50$v, size = sum(!J2001_ants_50$v), replace = FALSE)))))

# NEW: 50×50 modelo nulo CSR amplio: idem CSR estricto pero CON reemplazo y de tamaño determinado por una distribución Poisson con media igual al número observado de píxeles con nidos
J2001_vegdist_50_null[c("CSR2_mean", "CSR2_sd")] <- t(replicate(2000, 
    do.call(function(x) c(mean(x), sd(x)), list(
    sample(J2001_vegdist_50$v, size = rpois(1, sum(!J2001_ants_50$v)), replace = TRUE)))))

# 50×50 modelo nulo Poisson heterogéneo estricto: idem pero con diferente intensidad según microhábitat: muestras de tamaño N1 (nidos en pixeles con vegetación [vegmask_50]) + N2 (nidos en pixeles SD [sdmask])
# subsetting de imagen con distancias con las máscaras de veg y sd (devuelve vector numérico con drop = TRUE [el default])
# equivalente a ppm(J2001_ppp ~ J2001_veg_50) o ppm(J2001_ppp ~ J2001_sdmask_50) con intensidad heterogénea en función de vegetación (λ = nidos/pixeles, con un valor para c/u de los dos niveles) [Baddeley et al: 9.3.3–4]
J2001_vegdist_50_null[c("HET_mean", "HET_sd")] <- t(replicate(2000, 
    do.call(function(x) c(mean(x), sd(x)), list(c(
    sample(J2001_vegdist_50[J2001_vegmask_50], size = sum(!J2001_ants_50[J2001_vegmask_50]), replace = FALSE), 
    sample(J2001_vegdist_50[J2001_sdmask_50], size = sum(!J2001_ants_50[J2001_sdmask_50]), replace = FALSE))))))

# NEW: 50×50 modelo nulo Poisson heterogéneo amplio: idem pero CON reemplazo y de tamaño determinado por distribuciones Poisson con media igual al número observado de píxeles en cada categoría (máscara) veg y sd
J2001_vegdist_50_null[c("HET2_mean", "HET2_sd")] <- t(replicate(2000, 
    do.call(function(x) c(mean(x), sd(x)), list(c(
    sample(J2001_vegdist_50[J2001_vegmask_50], size = rpois(1, sum(!J2001_ants_50[J2001_vegmask_50])), replace = TRUE), 
    sample(J2001_vegdist_50[J2001_sdmask_50], size = rpois(1, sum(!J2001_ants_50[J2001_sdmask_50])), replace = TRUE))))))

# incluye a los valores observados del estadístico en la fila 1 de cada vector
J2001_vegdist_50_null[1, c(1, 3, 5, 7)] <- mean(J2001_vegdist_50[!J2001_ants_50$v])
J2001_vegdist_50_null[1, c(2, 4, 6, 8)] <- sd(J2001_vegdist_50[!J2001_ants_50$v])

# escala 100×100 eliminada (2021) [oldcode; change to modern version as above if needeed]
## 100×100 modelo nulo CSR: distribución de medias (o sd) de c/u de 2000 muestras aleatorias de tamaño N (píxeles con nidos) tomadas sin reemplazo de grilla con grano 50×50 y extendidas para promediar los 4 valores correspondientes del mapa de distancias con grano 100×100
## función que junta los cuatro valores de distancia (en la matriz "distveg") en 100×100 que corresponden a la muestra de píxeles elegidos ("px") entre aquellos con índice elegible en 50×50
# Sample_100x100 <- function(distveg, px) c(distveg$v[px], distveg$v[px + 1], distveg$v[px + 100], distveg$v[px + 101])
## llama a la función sample100x100 para cada muestra aleatoria de píxeles con los valores de índice seleccionables en 100×100
# J2001_vegdist_100_null$CSR_mean <- replicate(2000, mean(Sample_100x100(J2001_vegdist_100, sample(Selectable_100x100, size = sum(!J2001_ants_50$v), replace = FALSE))))
# J2001_vegdist_100_null$CSR_mean[1] <- mean(J2001_vegdist_100[!J2001_ants_100$v])
## idem para sd
# J2001_vegdist_100_null$CSR_sd <- replicate(2000, sd(Sample_100x100(J2001_vegdist_100, sample(Selectable_100x100, size = sum(!J2001_ants_50$v), replace = FALSE))))
# J2001_vegdist_100_null$CSR_sd[1] <- sd(J2001_vegdist_100[!J2001_ants_100$v])
## 100×100 modelo nulo Poisson heterogéneo: idem pero con diferente intensidad según microhábitat: muestras de tamaño N1 (nidos en pixeles con vegetación [veg_50$v == 0]) + N2 (nidos en pixeles SD [veg_50$v == 1])
## llama a la función sample100x100 para cada muestra aleatoria de píxeles con los valores de índice seleccionables en 100×100 (incluyendo su condición de vegetación)
# J2001_vegdist_100_null$HET_mean <- replicate(2000, mean(c(Sample_100x100(J2001_vegdist_100, sample(Selectable_100x100[as.logical(J2001_veg_50$v)], size = sum(!J2001_ants_50$v & J2001_veg_50$v), replace = FALSE)), Sample_100x100(J2001_vegdist_100, sample(Selectable_100x100[!J2001_veg_50$v], size = sum(!J2001_ants_50$v & !J2001_veg_50$v), replace = FALSE)))))
# J2001_vegdist_100_null$HET_mean[1] <- mean(J2001_vegdist_100[!J2001_ants_100$v])
## idem para sd
# J2001_vegdist_100_null$HET_sd <- replicate(2000, sd(c(Sample_100x100(J2001_vegdist_100, sample(Selectable_100x100[as.logical(J2001_veg_50$v)], size = sum(!J2001_ants_50$v & J2001_veg_50$v), replace = FALSE)), Sample_100x100(J2001_vegdist_100, sample(Selectable_100x100[!J2001_veg_50$v], size = sum(!J2001_ants_50$v & !J2001_veg_50$v), replace = FALSE)))))
# J2001_vegdist_100_null$HET_sd[1] <- sd(J2001_vegdist_100[!J2001_ants_100$v])

# J2019
# dataframe con los estadísticos derivados de las permutaciones (media y desvío) bajo cada modelo nulo, incluyendo al valor observado en la primera fila
# inicializando el dataframe con vectores vacíos
J2019_vegdist_50_null <- data.frame(CSR_mean = numeric(length = 2000), 
                                    CSR_sd = numeric(length = 2000), 
                                    CSR2_mean = numeric(length = 2000), 
                                    CSR2_sd = numeric(length = 2000), 
                                    HET_mean = numeric(length = 2000), 
                                    HET_sd = numeric(length = 2000), 
                                    HET2_mean = numeric(length = 2000), 
                                    HET2_sd = numeric(length = 2000))

#escala 100×100 eliminada (2021)
#J2019_vegdist_100_null <- data.frame(CSR_mean = numeric(length = 2000))

# código nuevo para obtener los dos parámetros evaluados (mean, sd) a partir de un mismo procedimiento de muestras aleatorias
# 50×50 modelo nulo CSR estricto: distribución de medias (y de sd) de c/u de 2000 muestras aleatorias independientes sin reemplazo de tamaño N (píxeles con nidos observados) del mapa de distancias
J2019_vegdist_50_null[c("CSR_mean", "CSR_sd")] <- t(replicate(2000, 
    do.call(function(x) c(mean(x), sd(x)), list(
    sample(J2019_vegdist_50$v, size = sum(!J2019_ants_50$v), replace = FALSE)))))

# NEW: 50×50 modelo nulo CSR amplio: idem CSR estricto pero CON reemplazo y de tamaño determinado por una distribución Poisson con media igual al número observado de píxeles con nidos
J2019_vegdist_50_null[c("CSR2_mean", "CSR2_sd")] <- t(replicate(2000, 
    do.call(function(x) c(mean(x), sd(x)), list(
    sample(J2019_vegdist_50$v, size = rpois(1, sum(!J2019_ants_50$v)), replace = TRUE)))))

# 50×50 modelo nulo Poisson heterogéneo estricto: idem pero con diferente intensidad según microhábitat: muestras de tamaño N1 (nidos en pixeles con vegetación [vegmask_50]) + N2 (nidos en pixeles SD [sdmask])
# subsetting de imagen con distancias con las máscaras de veg y sd (devuelve vector numérico con drop = TRUE [el default])
J2019_vegdist_50_null[c("HET_mean", "HET_sd")] <- t(replicate(2000, 
    do.call(function(x) c(mean(x), sd(x)), list(c(
    sample(J2019_vegdist_50[J2019_vegmask_50], size = sum(!J2019_ants_50[J2019_vegmask_50]), replace = FALSE), 
    sample(J2019_vegdist_50[J2019_sdmask_50], size = sum(!J2019_ants_50[J2019_sdmask_50]), replace = FALSE))))))

# NEW: 50×50 modelo nulo Poisson heterogéneo amplio: idem pero CON reemplazo y de tamaño determinado por distribuciones Poisson con media igual al número observado de píxeles con nidos en cada categoría (máscara) veg y sd
J2019_vegdist_50_null[c("HET2_mean", "HET2_sd")] <- t(replicate(2000, 
    do.call(function(x) c(mean(x), sd(x)), list(c(
    sample(J2019_vegdist_50[J2019_vegmask_50], size = rpois(1, sum(!J2019_ants_50[J2019_vegmask_50])), replace = TRUE), 
    sample(J2019_vegdist_50[J2019_sdmask_50], size = rpois(1, sum(!J2019_ants_50[J2019_sdmask_50])), replace = TRUE))))))

# incluye a los valores observados del estadístico en la fila 1 de cada vector
J2019_vegdist_50_null[1, c(1, 3, 5, 7)] <- mean(J2019_vegdist_50[!J2019_ants_50$v])
J2019_vegdist_50_null[1, c(2, 4, 6, 8)] <- sd(J2019_vegdist_50[!J2019_ants_50$v])

En la Figura 6 se presentan las distribuciones de distancias al borde observadas en la grilla J en 2001 y 2019 junto a las distribuciones esperadas bajo cuatro modelos nulos de tipo Poisson que combinan dos condiciones: modelos espacialmente homogéneos o heterogéneos (procesos espaciales de intensidad constante o variable respecto a la presencia de vegetación leñosa >1 m de altura) y más estrictos o más amplios o generales (número fijo o variable de píxeles con nidos efectivamente observados dentro de la parcela estudiada):

  • CSR (homogéneo, estricto): 2000 iteraciones de muestras de igual tamaño a la observada (= número de píxeles con nidos) tomadas SIN repetición de la distribución de distancias disponibles en la parcela,
  • CSR+ (homogéneo, general): idem pero cada muestra CON REPETICIÓN y de un tamaño aleatorio proveniente de una distribución Poisson con media igual al número observado de píxeles con nidos,
  • HET (heterogéneo, estricto): idem CSR estricto pero la muestra de distancias en cada iteración está constituida por dos submuestras: una de distancias de píxeles con vegetación leñosa (valores negativos) del mismo tamaño que el número de nidos observados bajo vegetación leñosa, y otra (de valores positivos) de tamaño igual al número de píxeles con nidos observados en suelo desnudo (lo que es equivalente a un muestreo de distancias con probabilidad dependiente de la variable categórica “vegetación leñosa”),
  • HET+ (heterogéneo, general): idem pero con muestras con repetición y de un tamaño aleatorio de cada submuestra proveniente de sendas distribuciones Poisson con medias iguales al número de píxeles con nidos observados en cada una de las dos categorías dadas por la vegetación leñosa.
# matriz de 8 subfiguras (una columna para cada año)
par(mar = c(2, 4, 3, 2))
layout(matrix(c(1, 2, 3, 4, 5, 6, 7, 8), 4, 2, byrow = FALSE), heights = c(2, 1.5, 1.5, 1.5))

# columna J2001
# 1: mapa de distancias con nidos
plot(J2001_vegdist_50, col = interp.colourmap(colourmap(c("darkgreen", "yellow", "darkred"), inputs = c(-5, 0, 5)), n = 255), ribside = "left", main = "J2001: distancia a bordes")
plot(J2001_ants_50, col = c(rgb(0, 0, 0, .6), rgb(1, 1, 1, 0)), ribbon = F, add = TRUE)

# histogramas de frecuencias
#hist(J2001_vegdist_50, breaks = seq(-5.125, 5.125, 1.025), col = "darkgreen", ylab = "Número de pixeles", xlab = "Distancia a borde", main = "")
#hist(J2001_vegdist_50[as.logical(!J2001_ants_50$v)], breaks = seq(-5.125, 5.125, 1.025), col = "black", add = TRUE)

# 2: histogramas como proporciones (con p-valor de test K-S)
# Prueba de K-S
J2001_vegdist_50_KS <- ks.test(J2001_vegdist_50$v, J2001_vegdist_50[as.logical(!J2001_ants_50$v)], exact = TRUE)
# distribución total
hist(J2001_vegdist_50, breaks = seq(-5.005, 5.005, 1.001), col = c( interp.colours(c("darkgreen", "yellow"), 5), interp.colours(c("yellow", "darkred"), 5)), ylab = "Proporción de pixeles", xlab = "Distancia a borde", ylim = c(0, 0.5), probability = TRUE, main = "")
# distribución observada (nidos)
hist(J2001_vegdist_50[as.logical(!J2001_ants_50$v)], breaks = seq(-5.005, 5.005, 1.001), col = rgb(0, 0, 0, .7), border = "black", density = 30, add = TRUE, probability = TRUE)
text(5, 0.5, paste("KS test:\nD =", round(J2001_vegdist_50_KS$stat, digits = 2), "\nP =", round(J2001_vegdist_50_KS$p, digits = 3)), adj = c(1, 1))

# 3: boxplot de la distribución de medias CSR, CSR2, HET y HET2
boxplot(c(J2001_vegdist_50_null$CSR_mean,
          J2001_vegdist_50_null$CSR2_mean,
          J2001_vegdist_50_null$HET_mean,
          J2001_vegdist_50_null$HET2_mean)
        ~ c(rep("CSR", 2000), rep("CSR+", 2000), rep("HET", 2000), rep("HET+", 2000)),
        horizontal = TRUE, ylim = c(-2, 2), xlim = c(4.75, 0.25), xlab = "", las = 1, main = "Media", ann = FALSE)
# valor observado (primer fila en cualquiera de las columnas)
abline(v = J2001_vegdist_50_null$CSR_mean[1], col = "red", lty = 3)
# pseudo-p values from permutations (double the number of more extreme values on the smallest [<0.5] tail)
text(2, 1, paste("P =", 2 * min(mean(J2001_vegdist_50_null$CSR_mean <= J2001_vegdist_50_null$CSR_mean[1]), mean(J2001_vegdist_50_null$CSR_mean >= J2001_vegdist_50_null$CSR_mean[1]))), adj = c(1, 0.5))

text(2, 2, paste("P =", 2 * min(mean(J2001_vegdist_50_null$CSR2_mean <= J2001_vegdist_50_null$CSR2_mean[1]), mean(J2001_vegdist_50_null$CSR2_mean >= J2001_vegdist_50_null$CSR2_mean[1]))), adj = c(1, 0.5))

text(2, 3, paste("P =", 2 * min(mean(J2001_vegdist_50_null$HET_mean <= J2001_vegdist_50_null$HET_mean[1]), mean(J2001_vegdist_50_null$HET_mean >= J2001_vegdist_50_null$HET_mean[1]))), adj = c(1, 0.5))

text(2, 4, paste("P =", 2 * min(mean(J2001_vegdist_50_null$HET2_mean <= J2001_vegdist_50_null$HET2_mean[1]), mean(J2001_vegdist_50_null$HET2_mean >= J2001_vegdist_50_null$HET2_mean[1]))), adj = c(1, 0.5))

# 4: boxplot de la distribución de sd CSR, CSR2, HET y HET2
boxplot(c(J2001_vegdist_50_null$CSR_sd,
          J2001_vegdist_50_null$CSR2_sd,
          J2001_vegdist_50_null$HET_sd,
          J2001_vegdist_50_null$HET2_sd)
        ~ c(rep("CSR", 2000), rep("CSR+", 2000), rep("HET", 2000), rep("HET+", 2000)),
        horizontal = TRUE, ylim = c(0, 3.5), xlim = c(4.75, 0.25), xlab = "", las = 1, main = "Desvío estándar", ann = FALSE)
# valor observado (primer fila en cualquiera de las columnas)
abline(v = J2001_vegdist_50_null$CSR_sd[1], col = "red", lty = 3)
# pseudo-p value from permutations (double the number of more extreme values on the smallest [<0.5] tail)
text(3.5, 1, paste("P =", 2 * min(mean(J2001_vegdist_50_null$CSR_sd <= J2001_vegdist_50_null$CSR_sd[1]), mean(J2001_vegdist_50_null$CSR_sd >= J2001_vegdist_50_null$CSR_sd[1]))), adj = c(1, 0.5))

text(3.5, 2, paste("P =", 2 * min(mean(J2001_vegdist_50_null$CSR2_sd <= J2001_vegdist_50_null$CSR2_sd[1]), mean(J2001_vegdist_50_null$CSR2_sd >= J2001_vegdist_50_null$CSR2_sd[1]))), adj = c(1, 0.5))

text(3.5, 3, paste("P =", 2 * min(mean(J2001_vegdist_50_null$HET_sd <= J2001_vegdist_50_null$HET_sd[1]), mean(J2001_vegdist_50_null$HET_sd >= J2001_vegdist_50_null$HET_sd[1]))), adj = c(1, 0.5))

text(3.5, 4, paste("P =", 2 * min(mean(J2001_vegdist_50_null$HET2_sd <= J2001_vegdist_50_null$HET2_sd[1]), mean(J2001_vegdist_50_null$HET2_sd >= J2001_vegdist_50_null$HET2_sd[1]))), adj = c(1, 0.5))

# escala 100×100 eliminada (2021)
# todo idem con 100×100
# plot(J2001_vegdist_100, col = interp.colourmap(colourmap(c("darkgreen", "yellow", "darkred"), inputs = c(-10, 0, 10)), n = 255), main = "J2001: distancia a bordes (100×100)")
# plot(J2001_ants_100, col = c(rgb(0, 0, 0, .6), rgb(1, 1, 1, 0)), ribbon = F, add = TRUE)
#hist(J2001_vegdist_100, breaks = seq(-9.135, 9.135, 1.015), col = "darkgreen", main = "", ylab = "Número de pixeles", xlab = "Distancia a borde")
#hist(J2001_vegdist_100[as.logical(!J2001_ants_100$v)], breaks = seq(-9.135, 9.135, 1.015), col = "black", add = TRUE)
# J2001_vegdist_100_KS <- ks.test(J2001_vegdist_100$v, J2001_vegdist_100[as.logical(!J2001_ants_100$v)])
# hist(J2001_vegdist_100, breaks = seq(-9.135, 9.135, 1.015), col = "darkgreen", main = "", ylab = "Proporción de pixeles", xlab = "Distancia a borde", ylim = c(0, 0.35), probability = TRUE)
# hist(J2001_vegdist_100[as.logical(!J2001_ants_100$v)], breaks = seq(-9.135, 9.135, 1.015), col = rgb(0, 0, 0, .7), border = "black", density = 30, add = TRUE, probability = TRUE)
# text(9, 0.35, paste("KS test:\nD =", round(J2001_vegdist_50_KS$stat, digits = 2), "\nP =", round(J2001_vegdist_50_KS$p, digits = 3)), adj = c(1, 1))
# boxplot(c(J2001_vegdist_100_null$CSR_mean, J2001_vegdist_100_null$HET_mean) ~ c(rep("CSR", 2000), rep("HET", 2000)), horizontal = TRUE, ylim = c(-9.135, 9.135), main = "Media")
# abline(v = J2001_vegdist_100_null$CSR_mean[1], col = "red", lty = 3)
# pseudo-p value from permutations (double the number of more extreme values on the smallest [<0.5] tail)
# text(9, 1, paste("P =", 2 * min(mean(J2001_vegdist_100_null$CSR_mean <= J2001_vegdist_100_null$CSR_mean[1]), mean(J2001_vegdist_100_null$CSR_mean >= J2001_vegdist_100_null$CSR_mean[1]))), adj = c(1, 0.5))
# text(9, 2, paste("P =", 2 * min(mean(J2001_vegdist_100_null$HET_mean <= J2001_vegdist_100_null$HET_mean[1]), mean(J2001_vegdist_100_null$HET_mean >= J2001_vegdist_100_null$HET_mean[1]))), adj = c(1, 0.5))
# boxplot de la distribución de sd CSR y HTR marcando el observado
# boxplot(c(J2001_vegdist_100_null$CSR_sd, J2001_vegdist_100_null$HET_sd) ~ c(rep("CSR", 2000), rep("HET", 2000)), horizontal = TRUE, ylim = c(0, 4), main = "Desvío estándar")
# abline(v = J2001_vegdist_100_null$CSR_sd[1], col = "red", lty = 3)
# pseudo-p value from permutations (double the number of more extreme values on the smallest [<0.5] tail)
# text(4, 1, paste("P =", 2 * min(mean(J2001_vegdist_100_null$CSR_sd <= J2001_vegdist_100_null$CSR_sd[1]), mean(J2001_vegdist_100_null$CSR_sd >= J2001_vegdist_100_null$CSR_sd[1]))), adj = c(1, 0.5))
# text(4, 2, paste("P =", 2 * min(mean(J2001_vegdist_100_null$HET_sd <= J2001_vegdist_100_null$HET_sd[1]), mean(J2001_vegdist_100_null$HET_sd >= J2001_vegdist_100_null$HET_sd[1]))), adj = c(1, 0.5))

# columna J2019
# 5: mapa de distancias con nidos
plot(J2019_vegdist_50, col = interp.colourmap(colourmap(c("darkgreen", "yellow", "darkred"), inputs = c(-5, 0, 5)), n = 255), ribside = "left", main = "J2019: distancia a bordes")
plot(J2019_ants_50, col = c(rgb(0, 0, 0, .6), rgb(1, 1, 1, 0)), ribbon = F, add = TRUE)

# histogramas de frecuencias
#hist(J2019_vegdist_50, breaks = seq(-5.125, 5.125, 1.025), col = "darkgreen", ylab = "Número de pixeles", xlab = "Distancia a borde", main = "")
#hist(J2019_vegdist_50[as.logical(!J2019_ants_50$v)], breaks = seq(-5.125, 5.125, 1.025), col = "black", add = TRUE)

# 6: histogramas como proporciones (con p-valor de test K-S)
# Prueba de K-S
J2019_vegdist_50_KS <- ks.test(J2019_vegdist_50$v, J2019_vegdist_50[as.logical(!J2019_ants_50$v)], exact = TRUE)
# distribución total
hist(J2019_vegdist_50, breaks = seq(-5.005, 5.005, 1.001), col = c( interp.colours(c("darkgreen", "yellow"), 5), interp.colours(c("yellow", "darkred"), 5)), ylab = "Proporción de pixeles", xlab = "Distancia a borde", ylim = c(0, 0.5), probability = TRUE, main = "")
# distribución observada (nidos)
hist(J2019_vegdist_50[as.logical(!J2019_ants_50$v)], breaks = seq(-5.005, 5.005, 1.001), col = rgb(0, 0, 0, .7), border = "black", density = 30, add = TRUE, probability = TRUE)
text(5, 0.5, paste("KS test:\nD =", round(J2019_vegdist_50_KS$stat, digits = 2), "\nP =", round(J2019_vegdist_50_KS$p, digits = 3)), adj = c(1, 1))

# 7: boxplot de la distribución de medias CSR, CSR2, HET y HET2
boxplot(c(J2019_vegdist_50_null$CSR_mean,
          J2019_vegdist_50_null$CSR2_mean,
          J2019_vegdist_50_null$HET_mean,
          J2019_vegdist_50_null$HET2_mean)
        ~ c(rep("CSR", 2000), rep("CSR+", 2000), rep("HET", 2000), rep("HET+", 2000)),
        horizontal = TRUE, ylim = c(-2, 2), xlim = c(4.75, 0.25), xlab = "", las = 1, main = "Media", ann = FALSE)
# valor observado (primer fila en cualquiera de las columnas)
abline(v = J2019_vegdist_50_null$CSR_mean[1], col = "red", lty = 3)
# pseudo-p values from permutations (double the number of more extreme values on the smallest [<0.5] tail)
text(2, 1, paste("P =", 2 * min(mean(J2019_vegdist_50_null$CSR_mean <= J2019_vegdist_50_null$CSR_mean[1]), mean(J2019_vegdist_50_null$CSR_mean >= J2019_vegdist_50_null$CSR_mean[1]))), adj = c(1, 0.5))

text(2, 2, paste("P =", 2 * min(mean(J2019_vegdist_50_null$CSR2_mean <= J2019_vegdist_50_null$CSR2_mean[1]), mean(J2019_vegdist_50_null$CSR2_mean >= J2019_vegdist_50_null$CSR2_mean[1]))), adj = c(1, 0.5))

text(2, 3, paste("P =", 2 * min(mean(J2019_vegdist_50_null$HET_mean <= J2019_vegdist_50_null$HET_mean[1]), mean(J2019_vegdist_50_null$HET_mean >= J2019_vegdist_50_null$HET_mean[1]))), adj = c(1, 0.5))

text(2, 4, paste("P =", 2 * min(mean(J2019_vegdist_50_null$HET2_mean <= J2019_vegdist_50_null$HET2_mean[1]), mean(J2019_vegdist_50_null$HET2_mean >= J2019_vegdist_50_null$HET2_mean[1]))), adj = c(1, 0.5))

# 8: boxplot de la distribución de sd CSR, CSR2, HET y HET2
boxplot(c(J2019_vegdist_50_null$CSR_sd,
          J2019_vegdist_50_null$CSR2_sd,
          J2019_vegdist_50_null$HET_sd,
          J2019_vegdist_50_null$HET2_sd)
        ~ c(rep("CSR", 2000), rep("CSR+", 2000), rep("HET", 2000), rep("HET+", 2000)),
        horizontal = TRUE, ylim = c(0, 3.5), xlim = c(4.75, 0.25), xlab = "", las = 1, main = "Desvío estándar", ann = FALSE)
# valor observado (primer fila en cualquiera de las columnas)
abline(v = J2019_vegdist_50_null$CSR_sd[1], col = "red", lty = 3)
# pseudo-p value from permutations (double the number of more extreme values on the smallest [<0.5] tail)
text(3.5, 1, paste("P =", 2 * min(mean(J2019_vegdist_50_null$CSR_sd <= J2019_vegdist_50_null$CSR_sd[1]), mean(J2019_vegdist_50_null$CSR_sd >= J2019_vegdist_50_null$CSR_sd[1]))), adj = c(1, 0.5))

text(3.5, 2, paste("P =", 2 * min(mean(J2019_vegdist_50_null$CSR2_sd <= J2019_vegdist_50_null$CSR2_sd[1]), mean(J2019_vegdist_50_null$CSR2_sd >= J2019_vegdist_50_null$CSR2_sd[1]))), adj = c(1, 0.5))

text(3.5, 3, paste("P =", 2 * min(mean(J2019_vegdist_50_null$HET_sd <= J2019_vegdist_50_null$HET_sd[1]), mean(J2019_vegdist_50_null$HET_sd >= J2019_vegdist_50_null$HET_sd[1]))), adj = c(1, 0.5))

text(3.5, 4, paste("P =", 2 * min(mean(J2019_vegdist_50_null$HET2_sd <= J2019_vegdist_50_null$HET2_sd[1]), mean(J2019_vegdist_50_null$HET2_sd >= J2019_vegdist_50_null$HET2_sd[1]))), adj = c(1, 0.5))
Arriba: Disposición espacial de los valores de distancia al borde calculados para la grilla J en 2001 (izquierda) y 2019 (derecha) según el mapa de vegetación leñosa, con los píxeles ocupados por nidos en negro. Luego, histogramas con la distribución de frecuencias (como proporciones) del conjunto de distancias disponibles (n = 2500 píxeles, en color) y de los correspondientes a sitios con nidos (rayado); los estadísticos corresponden a una prueba de Kolmogorov-Smirnov que pone a prueba la hipótesis nula de que ambos grupos provienen de una distribución común. Más abajo, estadísticos de posición (media) y dispersión (DE) de las distancias observadas (= píxeles con nidos; línea punteada roja) respecto a las distribuciones obtenidas bajo cuatro modelos nulos [ver texto] y la probabilidad (pseudo-p, estimada empíricamente a dos colas) de obtener un valor igual o más extremo que el observado bajo ese modelo nulo.

Figura 6: Arriba: Disposición espacial de los valores de distancia al borde calculados para la grilla J en 2001 (izquierda) y 2019 (derecha) según el mapa de vegetación leñosa, con los píxeles ocupados por nidos en negro. Luego, histogramas con la distribución de frecuencias (como proporciones) del conjunto de distancias disponibles (n = 2500 píxeles, en color) y de los correspondientes a sitios con nidos (rayado); los estadísticos corresponden a una prueba de Kolmogorov-Smirnov que pone a prueba la hipótesis nula de que ambos grupos provienen de una distribución común. Más abajo, estadísticos de posición (media) y dispersión (DE) de las distancias observadas (= píxeles con nidos; línea punteada roja) respecto a las distribuciones obtenidas bajo cuatro modelos nulos [ver texto] y la probabilidad (pseudo-p, estimada empíricamente a dos colas) de obtener un valor igual o más extremo que el observado bajo ese modelo nulo.

En la parcela J la distribución de distancias observadas entre los nidos y los bordes de la vegetación leñosa difieren notoriamente de la de las distancias disponibles (Fig. 6). En comparación con lo disponible o lo esperado por un modelo de independencia espacial completa (CSR), los píxeles con nidos son un subgrupo menos variable y con valores de distancia promedio más positivos (= hacia las áreas de suelo desnudo) respecto a lo esperado por azar (aún con una distribución de distancias menos acotada y un número de nidos similar pero aleatorio: CSR+). Eso parece deberse a una escasez relativa de nidos asociados los valores de distancia más negativos en contraste con una abundancia mayor a la esperada en las primeras categorías de distancias positivas. Esto es, un número relativamente alto de nidos en el suelo desnudo cercano a la vegetación leñosa a expensas de un número menor a lo esperado hacia adentro de los parches de vegetación.

Al incorporar intensidades diferentes del proceso espacial en los modelos aleatorios en función de la categoría de vegetación (modelo Poisson heterogéneo), como si hubiese una selección de sitios de nidificación a ese nivel, el sesgo de distancias observado disminuye, en especial bajo el modelo que aleatoriza el número efectivo de píxeles con nidos en cada tipo de vegetación (HET+) y, por lo tanto, amplía la variación entre los estadísticos calculados por cada una de las iteraciones. Los patrones son algo más notorios en 2001 (cuando los nidos eran más abundantes) que en 2019, donde la dispersión de las distancias de los nidos a los bordes de la vegetación parecen algo menos sesgadas.

En el caso de la otra parcela, solo está disponible el mapa de la vegetación leñosa presente en 2019. Para los nidos presentes en 2019 se pueden calcular las distancias al borde de la vegetación leñosa al igual que para la parcela anterior (con precisión dada por el grano del muestreo: 1 m2). Para los nidos presentes en 2001, en cambio, solo se pueden estimar distancias aproximadas suponiendo que ese tipo de vegetación perenne no se modificó demasiado. Es importante resaltar que la correlación entre las distancias desde cada punto muestreado al borde de vegetación en la grilla J en 2001 y en 2019 es r = 0.67, por lo que las resultados sobre V2001 que involucran a la vegetación deben ser tomadas con más reparos. Las distancias posibles al borde la vegetación en la grilla V y las distribuciones observadas y estimadas bajo cada uno de los modelos nulos planteados son:

# distancias a borde veg-sd de cada pixel
V2019_vegdist_50 <- distmap(V2019_vegmask_50) - distmap(V2019_sdmask_50)
sort(unique(as.vector(V2019_vegdist_50$v)))
##  [1] -6.414214 -6.242641 -5.828427 -5.656854 -5.414214 -5.242641 -5.000000
##  [8] -4.828427 -4.414214 -4.242641 -4.000000 -3.828427 -3.414214 -3.000000
## [15] -2.828427 -2.414214 -2.000000 -1.414214 -1.000000  1.000000  1.414214
## [22]  2.000000  2.414214  2.828427  3.000000  3.414214
## Modelos nulos
# como solo hay un set de distancias (V2019) ambos periodos (2001 y 2019) generarán distribuciones de distancias muestreando el mismo conjunto de distancias con tamaños de muestra diferentes
# para V2001
V2001_vegdist_50_null <- data.frame(CSR_mean = numeric(length = 2000), 
                                    CSR_sd = numeric(length = 2000), 
                                    CSR2_mean = numeric(length = 2000), 
                                    CSR2_sd = numeric(length = 2000), 
                                    HET_mean = numeric(length = 2000), 
                                    HET_sd = numeric(length = 2000), 
                                    HET2_mean = numeric(length = 2000), 
                                    HET2_sd = numeric(length = 2000))

# 50×50 modelo nulo CSR estricto
V2001_vegdist_50_null[c("CSR_mean", "CSR_sd")] <- t(replicate(2000, 
    do.call(function(x) c(mean(x), sd(x)), list(
    sample(V2019_vegdist_50$v, size = sum(!V2001_ants_50$v), replace = FALSE)))))

# 50×50 modelo nulo CSR amplio
V2001_vegdist_50_null[c("CSR2_mean", "CSR2_sd")] <- t(replicate(2000, 
    do.call(function(x) c(mean(x), sd(x)), list(
    sample(V2019_vegdist_50$v, size = rpois(1, sum(!V2001_ants_50$v)), replace = TRUE)))))

# 50×50 modelo nulo Poisson heterogéneo estricto
V2001_vegdist_50_null[c("HET_mean", "HET_sd")] <- t(replicate(2000, 
    do.call(function(x) c(mean(x), sd(x)), list(c(
    sample(V2019_vegdist_50[V2019_vegmask_50], size = sum(!V2001_ants_50[V2019_vegmask_50]), replace = FALSE), 
    sample(V2019_vegdist_50[V2019_sdmask_50], size = sum(!V2001_ants_50[V2019_sdmask_50]), replace = FALSE))))))

# 50×50 modelo nulo Poisson heterogéneo amplio
V2001_vegdist_50_null[c("HET2_mean", "HET2_sd")] <- t(replicate(2000, 
    do.call(function(x) c(mean(x), sd(x)), list(c(
    sample(V2019_vegdist_50[V2019_vegmask_50], size = rpois(1, sum(!V2001_ants_50[V2019_vegmask_50])), replace = TRUE), 
    sample(V2019_vegdist_50[V2019_sdmask_50], size = rpois(1, sum(!V2001_ants_50[V2019_sdmask_50])), replace = TRUE))))))

# incluye a los valores observados del estadístico en la fila 1 de cada vector
V2001_vegdist_50_null[1, c(1, 3, 5, 7)] <- mean(V2019_vegdist_50[!V2001_ants_50$v])
V2001_vegdist_50_null[1, c(2, 4, 6, 8)] <- sd(V2019_vegdist_50[!V2001_ants_50$v])

# idem para V2019
V2019_vegdist_50_null <- data.frame(CSR_mean = numeric(length = 2000), 
                                    CSR_sd = numeric(length = 2000), 
                                    CSR2_mean = numeric(length = 2000), 
                                    CSR2_sd = numeric(length = 2000), 
                                    HET_mean = numeric(length = 2000), 
                                    HET_sd = numeric(length = 2000), 
                                    HET2_mean = numeric(length = 2000), 
                                    HET2_sd = numeric(length = 2000))

# 50×50 modelo nulo CSR estricto
V2019_vegdist_50_null[c("CSR_mean", "CSR_sd")] <- t(replicate(2000, 
    do.call(function(x) c(mean(x), sd(x)), list(
    sample(V2019_vegdist_50$v, size = sum(!V2019_ants_50$v), replace = FALSE)))))

# 50×50 modelo nulo CSR amplio
V2019_vegdist_50_null[c("CSR2_mean", "CSR2_sd")] <- t(replicate(2000, 
    do.call(function(x) c(mean(x), sd(x)), list(
    sample(V2019_vegdist_50$v, size = rpois(1, sum(!V2019_ants_50$v)), replace = TRUE)))))

# 50×50 modelo nulo Poisson heterogéneo estricto
V2019_vegdist_50_null[c("HET_mean", "HET_sd")] <- t(replicate(2000, 
    do.call(function(x) c(mean(x), sd(x)), list(c(
    sample(V2019_vegdist_50[V2019_vegmask_50], size = sum(!V2019_ants_50[V2019_vegmask_50]), replace = FALSE), 
    sample(V2019_vegdist_50[V2019_sdmask_50], size = sum(!V2019_ants_50[V2019_sdmask_50]), replace = FALSE))))))

# 50×50 modelo nulo Poisson heterogéneo amplio
V2019_vegdist_50_null[c("HET2_mean", "HET2_sd")] <- t(replicate(2000, 
    do.call(function(x) c(mean(x), sd(x)), list(c(
    sample(V2019_vegdist_50[V2019_vegmask_50], size = rpois(1, sum(!V2019_ants_50[V2019_vegmask_50])), replace = TRUE), 
    sample(V2019_vegdist_50[V2019_sdmask_50], size = rpois(1, sum(!V2019_ants_50[V2019_sdmask_50])), replace = TRUE))))))

# incluye a los valores observados del estadístico en la fila 1 de cada vector
V2019_vegdist_50_null[1, c(1, 3, 5, 7)] <- mean(V2019_vegdist_50[!V2019_ants_50$v])
V2019_vegdist_50_null[1, c(2, 4, 6, 8)] <- sd(V2019_vegdist_50[!V2019_ants_50$v])

## FIGURAS (una columna para cada año)
par(mar = c(2, 4, 3, 2))
layout(matrix(c(1, 2, 3, 4, 5, 6, 7, 8), 4, 2, byrow = FALSE), heights = c(2, 1.5, 1.5, 1.5))

# columna V2001
# 1: mapa de distancias con nidos
plot(V2019_vegdist_50, col = interp.colourmap(colourmap(c("darkgreen", "yellow", "darkred"), inputs = c(-7, 0, 7)), n = 255), ribside = "left", main = "V2001: distancia a bordes en 2019")
plot(V2001_ants_50, col = c(rgb(0, 0, 0, .6), rgb(1, 1, 1, 0)), ribbon = F, add = TRUE)

# 2: histogramas como proporciones (con p-valor de test K-S)
# Prueba de K-S
V2001_vegdist_50_KS <- ks.test(V2019_vegdist_50$v, V2019_vegdist_50[as.logical(!V2001_ants_50$v)], exact = TRUE)
# distribución total
hist(V2019_vegdist_50, breaks = seq(-7.007, 7.007, 1.001), col = c( interp.colours(c("darkgreen", "yellow"), 7), interp.colours(c("yellow", "darkred"), 4)), ylab = "Proporción de pixeles", xlab = "Distancia a borde", ylim = c(0, 0.5), probability = TRUE, main = "")
# distribución observada (nidos)
hist(V2019_vegdist_50[as.logical(!V2001_ants_50$v)], breaks = seq(-7.007, 7.007, 1.001), col = rgb(0, 0, 0, .7), border = "black", density = 30, add = TRUE, probability = TRUE)
text(7, 0.5, paste("KS test:\nD =", round(V2001_vegdist_50_KS$stat, digits = 2), "\nP =", round(V2001_vegdist_50_KS$p, digits = 3)), adj = c(1, 1))

# 3: boxplot de la distribución de medias CSR, CSR2, HET y HET2
boxplot(c(V2001_vegdist_50_null$CSR_mean,
          V2001_vegdist_50_null$CSR2_mean,
          V2001_vegdist_50_null$HET_mean,
          V2001_vegdist_50_null$HET2_mean)
        ~ c(rep("CSR", 2000), rep("CSR+", 2000), rep("HET", 2000), rep("HET+", 2000)),
        horizontal = TRUE, ylim = c(-2, 2), xlim = c(4.75, 0.25), las = 1, main = "Media", ann = FALSE)
# valor observado (primer fila en cualquiera de las columnas)
abline(v = V2001_vegdist_50_null$CSR_mean[1], col = "red", lty = 3)
# pseudo-p values from permutations (double the number of more extreme values on the smallest [<0.5] tail)
text(2, 1, paste("P =", 2 * min(mean(V2001_vegdist_50_null$CSR_mean <= V2001_vegdist_50_null$CSR_mean[1]), mean(V2001_vegdist_50_null$CSR_mean >= V2001_vegdist_50_null$CSR_mean[1]))), adj = c(1, 0.5))

text(2, 2, paste("P =", 2 * min(mean(V2001_vegdist_50_null$CSR2_mean <= V2001_vegdist_50_null$CSR2_mean[1]), mean(V2001_vegdist_50_null$CSR2_mean >= V2001_vegdist_50_null$CSR2_mean[1]))), adj = c(1, 0.5))

text(2, 3, paste("P =", 2 * min(mean(V2001_vegdist_50_null$HET_mean <= V2001_vegdist_50_null$HET_mean[1]), mean(V2001_vegdist_50_null$HET_mean >= V2001_vegdist_50_null$HET_mean[1]))), adj = c(1, 0.5))

text(2, 4, paste("P =", 2 * min(mean(V2001_vegdist_50_null$HET2_mean <= V2001_vegdist_50_null$HET2_mean[1]), mean(V2001_vegdist_50_null$HET2_mean >= V2001_vegdist_50_null$HET2_mean[1]))), adj = c(1, 0.5))

# 4: boxplot de la distribución de sd CSR, CSR2, HET y HET2
boxplot(c(V2001_vegdist_50_null$CSR_sd,
          V2001_vegdist_50_null$CSR2_sd,
          V2001_vegdist_50_null$HET_sd,
          V2001_vegdist_50_null$HET2_sd)
        ~ c(rep("CSR", 2000), rep("CSR+", 2000), rep("HET", 2000), rep("HET+", 2000)),
        horizontal = TRUE, ylim = c(0, 3.5), xlim = c(4.75, 0.25), las = 1, main = "Desvío estándar", ann = FALSE)
# valor observado (primer fila en cualquiera de las columnas)
abline(v = V2001_vegdist_50_null$CSR_sd[1], col = "red", lty = 3)
# pseudo-p value from permutations (double the number of more extreme values on the smallest [<0.5] tail)
text(3.5, 1, paste("P =", 2 * min(mean(V2001_vegdist_50_null$CSR_sd <= V2001_vegdist_50_null$CSR_sd[1]), mean(V2001_vegdist_50_null$CSR_sd >= V2001_vegdist_50_null$CSR_sd[1]))), adj = c(1, 0.5))

text(3.5, 2, paste("P =", 2 * min(mean(V2001_vegdist_50_null$CSR2_sd <= V2001_vegdist_50_null$CSR2_sd[1]), mean(V2001_vegdist_50_null$CSR2_sd >= V2001_vegdist_50_null$CSR2_sd[1]))), adj = c(1, 0.5))

text(3.5, 3, paste("P =", 2 * min(mean(V2001_vegdist_50_null$HET_sd <= V2001_vegdist_50_null$HET_sd[1]), mean(V2001_vegdist_50_null$HET_sd >= V2001_vegdist_50_null$HET_sd[1]))), adj = c(1, 0.5))

text(3.5, 4, paste("P =", 2 * min(mean(V2001_vegdist_50_null$HET2_sd <= V2001_vegdist_50_null$HET2_sd[1]), mean(V2001_vegdist_50_null$HET2_sd >= V2001_vegdist_50_null$HET2_sd[1]))), adj = c(1, 0.5))

# columna V2019
# 5: mapa de distancias con nidos
plot(V2019_vegdist_50, col = interp.colourmap(colourmap(c("darkgreen", "yellow", "darkred"), inputs = c(-7, 0, 7)), n = 255), ribside = "left", main = "V2019: distancia a bordes")
plot(V2019_ants_50, col = c(rgb(0, 0, 0, .6), rgb(1, 1, 1, 0)), ribbon = F, add = TRUE)

# 6: histogramas como proporciones (con p-valor de test K-S)
# Prueba de K-S
V2019_vegdist_50_KS <- ks.test(V2019_vegdist_50$v, V2019_vegdist_50[as.logical(!V2019_ants_50$v)], exact = TRUE)
# distribución total
hist(V2019_vegdist_50, breaks = seq(-7.007, 7.007, 1.001), col = c( interp.colours(c("darkgreen", "yellow"), 7), interp.colours(c("yellow", "darkred"), 4)), ylab = "Proporción de pixeles", xlab = "Distancia a borde", ylim = c(0, 0.5), probability = TRUE, main = "")
# distribución observada (nidos)
hist(V2019_vegdist_50[as.logical(!V2019_ants_50$v)], breaks = seq(-7.007, 7.007, 1.001), col = rgb(0, 0, 0, .7), border = "black", density = 30, add = TRUE, probability = TRUE)
text(7, 0.5, paste("KS test:\nD =", round(V2019_vegdist_50_KS$stat, digits = 2), "\nP =", round(V2019_vegdist_50_KS$p, digits = 3)), adj = c(1, 1))

# 7: boxplot de la distribución de medias CSR, CSR2, HET y HET2
boxplot(c(V2019_vegdist_50_null$CSR_mean,
          V2019_vegdist_50_null$CSR2_mean,
          V2019_vegdist_50_null$HET_mean,
          V2019_vegdist_50_null$HET2_mean)
        ~ c(rep("CSR", 2000), rep("CSR+", 2000), rep("HET", 2000), rep("HET+", 2000)),
        horizontal = TRUE, ylim = c(-2, 2), xlim = c(4.75, 0.25), las = 1, main = "Media", ann = FALSE)
# valor observado (primer fila en cualquiera de las columnas)
abline(v = V2019_vegdist_50_null$CSR_mean[1], col = "red", lty = 3)
# pseudo-p values from permutations (double the number of more extreme values on the smallest [<0.5] tail)
text(2, 1, paste("P =", 2 * min(mean(V2019_vegdist_50_null$CSR_mean <= V2019_vegdist_50_null$CSR_mean[1]), mean(V2019_vegdist_50_null$CSR_mean >= V2019_vegdist_50_null$CSR_mean[1]))), adj = c(1, 0.5))

text(2, 2, paste("P =", 2 * min(mean(V2019_vegdist_50_null$CSR2_mean <= V2019_vegdist_50_null$CSR2_mean[1]), mean(V2019_vegdist_50_null$CSR2_mean >= V2019_vegdist_50_null$CSR2_mean[1]))), adj = c(1, 0.5))

text(2, 3, paste("P =", 2 * min(mean(V2019_vegdist_50_null$HET_mean <= V2019_vegdist_50_null$HET_mean[1]), mean(V2019_vegdist_50_null$HET_mean >= V2019_vegdist_50_null$HET_mean[1]))), adj = c(1, 0.5))

text(2, 4, paste("P =", 2 * min(mean(V2019_vegdist_50_null$HET2_mean <= V2019_vegdist_50_null$HET2_mean[1]), mean(V2019_vegdist_50_null$HET2_mean >= V2019_vegdist_50_null$HET2_mean[1]))), adj = c(1, 0.5))

# 8: boxplot de la distribución de sd CSR, CSR2, HET y HET2
boxplot(c(V2019_vegdist_50_null$CSR_sd,
          V2019_vegdist_50_null$CSR2_sd,
          V2019_vegdist_50_null$HET_sd,
          V2019_vegdist_50_null$HET2_sd)
        ~ c(rep("CSR", 2000), rep("CSR+", 2000), rep("HET", 2000), rep("HET+", 2000)),
        horizontal = TRUE, ylim = c(0, 3.5), xlim = c(4.75, 0.25), las = 1, main = "Desvío estándar", ann = FALSE)
# valor observado (primer fila en cualquiera de las columnas)
abline(v = V2019_vegdist_50_null$CSR_sd[1], col = "red", lty = 3)
# pseudo-p value from permutations (double the number of more extreme values on the smallest [<0.5] tail)
text(3.5, 1, paste("P =", 2 * min(mean(V2019_vegdist_50_null$CSR_sd <= V2019_vegdist_50_null$CSR_sd[1]), mean(V2019_vegdist_50_null$CSR_sd >= V2019_vegdist_50_null$CSR_sd[1]))), adj = c(1, 0.5))

text(3.5, 2, paste("P =", 2 * min(mean(V2019_vegdist_50_null$CSR2_sd <= V2019_vegdist_50_null$CSR2_sd[1]), mean(V2019_vegdist_50_null$CSR2_sd >= V2019_vegdist_50_null$CSR2_sd[1]))), adj = c(1, 0.5))

text(3.5, 3, paste("P =", 2 * min(mean(V2019_vegdist_50_null$HET_sd <= V2019_vegdist_50_null$HET_sd[1]), mean(V2019_vegdist_50_null$HET_sd >= V2019_vegdist_50_null$HET_sd[1]))), adj = c(1, 0.5))

text(3.5, 4, paste("P =", 2 * min(mean(V2019_vegdist_50_null$HET2_sd <= V2019_vegdist_50_null$HET2_sd[1]), mean(V2019_vegdist_50_null$HET2_sd >= V2019_vegdist_50_null$HET2_sd[1]))), adj = c(1, 0.5))
Idem Fig. \@ref(fig:FiguraDistanciasJ) para la parcela V.

Figura 7: Idem Fig. 6 para la parcela V.

En este caso, (y dándole más énfasis a los patrones 2019 dado que la vegetación fue medida en ese momento) la distribución de distancias observadas no difiere de las disponibles pero se mantiene la tendencia hacia mayor uso de suelo desnudo y de los sitios más cercanos al borde de la vegetación (7). Las distancias de los píxeles con nidos parecen ser un subgrupo con valores algo más positivos que lo esperado por modelos espacialmente homogéneos (CSR, CSR+), un sesgo que disminuye al considerar lo predicho por modelos heterogéneos (HET, HET+), y siempre menos variable que lo esperado por azar bajo todos los modelos aleatorios (un sesgo que disminuye en magnitud en 2001, cuando había más nidos).

#Listas para bucle de subfiguras
All_vegdist_50 <- list("J 2001" = J2001_vegdist_50, "J 2019" = J2019_vegdist_50, "V 2001" = V2019_vegdist_50, "V 2019" = V2019_vegdist_50)
All_vegdist_50_KS <- list("J 2001" = NULL, "J 2019" = NULL, "V 2001" = NULL, "V 2019" = NULL)
All_vegdist_50_null <- list("J 2001" = J2001_vegdist_50_null, "J 2019" = J2019_vegdist_50_null, "V 2001" = V2001_vegdist_50_null, "V 2019" = V2019_vegdist_50_null)

# matriz de 8 subfiguras (una columna para cada año)
layout(matrix(c(1:16), 4, 4, byrow = FALSE), heights = c(2.5, 2, 1, 1))

for(i in 1:length(All_ants_50)) {

  # 1: mapa de distancias con nidos
  par(mar = c(0, 1.3, 1, 0))
  plot(All_vegdist_50[[i]], col = interp.colourmap(colourmap(c("darkgreen", "yellow", "darkred"), inputs = c(-6.5, 0, 3.5)), n = 255), ribside = "left", ribsep = 0.05, main = names(All_ants_50)[i])
  plot(All_ants_50[[i]], col = c(rgb(0, 0, 0, .8), rgb(1, 1, 1, 0)), ribbon = F, add = TRUE)
  
  # 2: histogramas como proporciones (con p-valor de test K-S)
  par(mar = c(4, 3.2, 0, 0.5), mgp = c(2, 0.75, 0))
  # distribución total
  hist(All_vegdist_50[[i]], breaks = seq(-7.007, 4.004, 1.001), col = c( interp.colours(c("darkgreen", "yellow"), 7), interp.colours(c("yellow", "darkred"), 4)), ylab = "Pixels (proportion)", xlab = "Distance to patch border (m)", cex.lab = 1.2, ylim = c(0, 0.5), probability = TRUE, main = "")
  # distribución observada (nidos)
  hist(All_vegdist_50[[i]][as.logical(!All_ants_50[[i]]$v)], breaks = seq(-7.007, 4.004, 1.001), col = rgb(0, 0, 0, .8), border = "black", density = 30, add = TRUE, probability = TRUE)
  # Prueba de K-S con p-valor exacto
  All_vegdist_50_KS[[i]] <- ks.test(All_vegdist_50[[i]]$v, All_vegdist_50[[i]][as.logical(!All_ants_50[[i]]$v)], exact = TRUE)
  text(-6.5, 0.5, paste("D =", round(All_vegdist_50_KS[[i]]$stat, digits = 2), "\nP =", round(All_vegdist_50_KS[[i]]$p, digits = 3)), adj = c(0, 1))
  
  # 3: boxplot de la distribución de medias CSR
  par(mar = c(3.5, 1.5, 0.5, 0.5), mgp = c(1.75, 0.75, 0))
  boxplot(All_vegdist_50_null[[i]]$CSR_mean ~ c(rep("CSR", 2000)),
          horizontal = TRUE, ylim = c(-2, 2), xlab = "Mean distance (m)", ylab = "", cex.lab = 1.2, las = 1, frame.plot = FALSE)
  # valor observado (primer fila en cualquiera de las columnas)
  abline(v = All_vegdist_50_null[[i]]$CSR_mean[1], col = "red", lty = 1, lwd = 1.5)
  # pseudo-p values from permutations (double the number of more extreme values on the smallest [<0.5] tail)
  text(2, 1, paste("P =", 2 * min(mean(All_vegdist_50_null[[i]]$CSR_mean <= All_vegdist_50_null[[i]]$CSR_mean[1]), mean(All_vegdist_50_null[[i]]$CSR_mean >= All_vegdist_50_null[[i]]$CSR_mean[1]))), adj = c(1, 0.5))
  
  # 4: boxplot de la distribución de sd CSR
  boxplot(All_vegdist_50_null[[i]]$CSR_sd ~ rep("CSR", 2000),
          horizontal = TRUE, ylim = c(0, 3.5), xlab = "SD distance (m)", ylab = "", cex.lab = 1.2, las = 1, frame.plot = FALSE)
  # valor observado (primer fila en cualquiera de las columnas)
  abline(v = All_vegdist_50_null[[i]]$CSR_sd[1], col = "red", lty = 1, lwd = 2)
  # pseudo-p value from permutations (double the number of more extreme values on the smallest [<0.5] tail)
  text(3.5, 1, paste("P =", 2 * min(mean(All_vegdist_50_null[[i]]$CSR_sd <= All_vegdist_50_null[[i]]$CSR_sd[1]), mean(All_vegdist_50_null[[i]]$CSR_sd >= All_vegdist_50_null[[i]]$CSR_sd[1]))), adj = c(1, 0.5))
}
Idem Figs. \@ref(fig:FiguraDistanciasJ) y \@ref(fig:DistanciasModelosFigsV) pero reordenadas, con escalas comunes y solo comparando contra CSR como modelo nulo.

Figura 8: Idem Figs. 6 y 7 pero reordenadas, con escalas comunes y solo comparando contra CSR como modelo nulo. 

#back to defaults
par(mar = c(5, 4, 4, 2) + 0.1, mgp = c(3, 1, 0))

Ajuste de modelos Poisson

Si ajustamos esos modelos por máxima verosimilitud (función spatstat::ppm), asumiendo que no hay interacción entre nidos, la/s intensidad/es estimadas para el modelo homogéneo y el heterogéneo (con dos categorías: áreas con vegetación leñosa o con suelo desnudo) para el caso de J2001 son:

# los parámetros para los puntos extras (dummy) necesarios para los ajuste se fijan en los posibles según este muestreo (2500 píxeles) mediante quadscheme(J2001_ppp, nd = 50), o pixelquad(J2001_ppp, J2001_ants_50), que ya establece la ventana de análisis y su resolución [ver Baddeley et al: 9.8.1.3–4]
# lo mismo se consigue en este caso indicando que divida el rectángulo en 50x50 directamente agregando la opción nd = 50 en el llamado a ppm: ppm(J2001_ppp ~ 1, nd = 50)

# Modelo CSR
ModJ2001_CSR <- ppm(quadscheme(J2001_ppp, nd = 50) ~ 1)
ModJ2001_CSR

# Modelo HET categórico: vegetación
ModJ2001_Veg <- ppm(quadscheme(J2001_ppp, nd = 50) ~ J2001_veg_50)
ModJ2001_Veg
anova(ModJ2001_CSR, ModJ2001_Veg, test = "LRT")
## Stationary Poisson process
## Fitted to point pattern dataset 'quadscheme(J2001_ppp, nd = 50)'
## Intensity: 0.0304
##              Estimate      S.E.   CI95.lo   CI95.hi Ztest      Zval
## log(lambda) -3.493313 0.1147079 -3.718136 -3.268489   *** -30.45399
## Nonstationary Poisson process
## Fitted to point pattern dataset 'quadscheme(J2001_ppp, nd = 50)'
## 
## Log intensity:  ~J2001_veg_50
## 
## Fitted trend coefficients:
##  (Intercept) J2001_veg_50 
##    -4.070341     1.097877 
## 
##               Estimate      S.E.   CI95.lo   CI95.hi Ztest       Zval
## (Intercept)  -4.070341 0.1961161 -4.454721 -3.685960   *** -20.754747
## J2001_veg_50  1.097877 0.2417882  0.623981  1.571773   ***   4.540656
## Analysis of Deviance Table
## 
## Model 1: ~1   Poisson
## Model 2: ~J2001_veg_50    Poisson
##   Npar Df Deviance  Pr(>Chi)    
## 1    1                          
## 2    2  1   22.079 2.616e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El test de cociente de verosimilitud entre los dos modelos (anidados) muestra que el modelo heterogéneo categórico (veg vs sd) tiene un ajuste mucho mejor que el homogéneo. La misma conclusión se obtiene comparando sus valores de AIC (CSR: 685; HET: 664.9). [de aquí en más compararemos el ajuste de los diferentes modelos solo mediante sus valores de AIC]

También podemos estimar modelos donde la heterogeneidad de la intensidad del proceso espacial esté determinada solo por la distancia al borde de la vegetación (BRD), o por la distancia al borde pero con distinta relación potencial dependiendo de si es hacia dentro del parche de vegetación leñosa o hacia el parche de suelo desnudo. Esto se puede modelar con relaciones distintas a cada lado (BRDHET) o mediante una función cuadrática sobre la variable continua de distancia (DIST²). La forma de las distintas relaciones estimadas entre la intensidad del proceso y esas variables explicativas se muestran en la Fig. 9.

# Modelo het continuo: distancia a borde (independiente del tipo de vegetación)
ModJ2001_AbsDist <- ppm(quadscheme(J2001_ppp, nd = 50) ~ abs(J2001_vegdist_50))

# figura de efectos predichos sobre la intensidad por c/modelo
par(mfrow = c(1, 4), mar = c(2, 2, 1.5, 1))
rangox = seq(min(J2001_vegdist_50), max(J2001_vegdist_50), 0.2)

barplot(c(exp(coef(ModJ2001_Veg)[1]), exp(sum(coef(ModJ2001_Veg)))), names.arg = c("veg", "sd"), ylab = "Intensidad (λ)", main = "J2001:HET", col = "white")
plot(effectfun(ModJ2001_AbsDist), main = "J2001:BRD")

# Modelos HET interacción (distancia y vegetación)
## interacción entre distancia absoluta y veg, con intercepto común (factor anidado en continua; ver Baddeley et al: 9.3.9.4)
ModJ2001_DistxVeg <- ppm(quadscheme(J2001_ppp, nd = 50) ~ abs(J2001_vegdist_50) / J2001_veg_50)
plot(x = rangox[rangox <= 0], y = exp(coef(ModJ2001_DistxVeg)[1] - coef(ModJ2001_DistxVeg)[2] * rangox[rangox <= 0]), type = "l", ylab = "Intensidad (λ)", main = "J2001:BRDHET", xlim = range(rangox), ylim = c(0, 0.07))
lines(x = rangox[rangox >= 0], y = exp(sum(coef(ModJ2001_DistxVeg)[1]) + sum(coef(ModJ2001_DistxVeg)[c(2, 3)]) * rangox[rangox >= 0]))

## cuadrática de la distancia con signo (=> máximo a cualquier valor intermedio, no necesariamente cero)
ModJ2001_Dist2 <- ppm(quadscheme(J2001_ppp, nd = 50) ~ polynom(J2001_vegdist_50, 2))
plot(effectfun(ModJ2001_Dist2), main = "J2001:DIST²")
Relaciones estimadas entre la intensidad del proceso espacial Poisson que determina la presencia de colonias de *P. bergi* en la parcela J durante 2001 y los parches de vegetación leñosa.

Figura 9: Relaciones estimadas entre la intensidad del proceso espacial Poisson que determina la presencia de colonias de P. bergi en la parcela J durante 2001 y los parches de vegetación leñosa. 

ModsInhomJ2001 <- list("J2001:HET" = ModJ2001_Veg,
                       "J2001:BRD" = ModJ2001_AbsDist,
                       "J2001:BRDHET" = ModJ2001_DistxVeg,
                       "J2001:DIST²" = ModJ2001_Dist2)

El resumen de los parámetros estimados por cada modelo es:

Modelos ajustados: J2001

Modelo Objeto Dependencia Intensidad AIC
CSR ModJ2001_CSR — 0.03 684.98
HET ModJ2001_Veg veg/sd λveg: 0.017 / λsd: 0.051 664.9
BRD ModJ2001_AbsDist abs(dist) 0.079 × 0.462m 678.61
BRDHET ModJ2001_DistxVeg abs(dist) * veg/sd λveg: 0.059 × 0.37m / λsd: 0.059 × 0.896m 659.74
DIST2 ModJ2001_Dist2 dist2 0.04 × 1.553m × 0.845m2 659.79

Cualquiera de los dos modelos que propone una relación asimétrica de los nidos con la distancia al borde dependiendo del tipo de cobertura (BRDHET, DIST2) ajusta mejor al patrón de puntos observado [en negrita] que los más simples (i.e., CSR: homogéneo, HET: solo tipo de cobertura, BRD: solo distancia).

Lo mismo puede hacerse para la distribución de nidos en el mismo lugar (J) en 2019 (Fig. 10).

par(mfrow = c(1, 4), mar = c(2, 2, 1.5, 1))
rangox = seq(min(J2019_vegdist_50), max(J2019_vegdist_50), 0.2)

# Modelo CSR
ModJ2019_CSR <- ppm(quadscheme(J2019_ppp, nd = 50) ~ 1)

# Modelo HET categórico: vegetación
ModJ2019_Veg <- ppm(quadscheme(J2019_ppp, nd = 50) ~ J2019_veg_50)
barplot(c(exp(coef(ModJ2019_Veg)[1]), exp(sum(coef(ModJ2019_Veg)))), names.arg = c("veg", "sd"), ylab = "Intensidad (λ)", main = "J2019:HET", col ="white")

# Modelo HET continuo: distancia a borde (independiente del tipo de vegetación)
ModJ2019_AbsDist <- ppm(quadscheme(J2019_ppp, nd = 50) ~ abs(J2019_vegdist_50))
plot(effectfun(ModJ2019_AbsDist), main = "J2019:BRD")

# Modelos HET interacción (distancia y vegetación)
## interacción distancia absoluta / veg
ModJ2019_DistxVeg <- ppm(quadscheme(J2019_ppp, nd = 50) ~ abs(J2019_vegdist_50) / J2019_veg_50)
plot(x = rangox[rangox <= 0], y = exp(coef(ModJ2019_DistxVeg)[1] - coef(ModJ2019_DistxVeg)[2] * rangox[rangox <= 0]), type = "l", ylab = "Intensidad (λ)", main = "J2019:BRDHET", xlim = range(rangox), ylim = c(0, 0.04))
lines(x = rangox[rangox >= 0], y = exp(sum(coef(ModJ2019_DistxVeg)[1]) + sum(coef(ModJ2019_DistxVeg)[c(2, 3)]) * rangox[rangox >= 0]))

## cuadrática de la distancia con signo (=> máximo a cualquier valor intermedio, no necesariamente cero)
ModJ2019_Dist2 <- ppm(quadscheme(J2019_ppp, nd = 50) ~ polynom(J2019_vegdist_50, 2))
plot(effectfun(ModJ2019_Dist2), main = "J2019:DIST²")
Idem Fig. \@ref(fig:MasModelosJ2001) para la parcela J durante 2019.

Figura 10: Idem Fig. 9 para la parcela J durante 2019. 

ModsInhomJ2019 <- list("J2019:HET" = ModJ2019_Veg,
                       "J2019:BRD" = ModJ2019_AbsDist,
                       "J2019:BRDHET" = ModJ2019_DistxVeg,
                       "J2019:DIST²" = ModJ2019_Dist2)

Modelos ajustados: J2019

Modelo Objeto Dependencia Intensidad AIC
CSR ModJ2019_CSR — 0.013 344.93
HET ModJ2019_Veg veg/sd λveg: 0.009 / λsd: 0.022 340.35
BRD ModJ2019_AbsDist abs(dist) 0.042 × 0.423m 339.3
BRDHET ModJ2019_DistxVeg abs(dist) * veg/sd λveg: 0.028 × 0.427m / λsd: 0.028 × 0.847m 335.65
DIST2 ModJ2019_Dist2 dist2 0.019 × 1.404m × 0.868m2 335.66

De nuevo, los modelos con relación asimétrica al borde de la vegetación son los que mejor ajustan [en negrita] al patrón observado (siempre suponiendo que no hay interacción relevante entre los nidos a esta escala).

Podemos, entonces, generar distribuciones de los estadísticos de posición y de dispersión bajo modelos nulos que incluyan esas relaciones entre la intensidad del proceso espacial tipo Poisson y la distancia/presencia de vegetación leñosa para poder compararlos con los valores observados:

  • BRD: la intensidad del proceso espacial (= la probabilidad de formar parte de la muestra elegida en cada iteración) es proporcional a la distancia al píxel más cercano con el otro tipo de vegetación (“distancia al borde”), independientemente de cuál sea (valor absoluto de distancia),
  • BRD+: idem BRD pero con muestras con repetición y un tamaño de muestra aleatorio proveniente de una distribución Poisson con media igual a la observada,
  • BRDHET / BRDHET+: idem BRD pero con (potencialmente) distinta relación con la distancia al borde dependiendo del tipo de vegetación,
  • DIST2 / DIST2+: la intensidad del proceso está relacionada por una función cuadrática con la distancia al borde, que tiene valores negativos hacia dentro de la vegetación leñosa y positivos hacia el parche de suelo desnudo.

La distribución de valores de distancia al borde de la vegetación de los píxeles con nidos, que es muy sesgada respecto a un modelo CSR, deja de presentar valores extremos de sus estadísticos de posición y dispersión al compararlos con lo esperado por modelos de tipo Poisson con intensidades variables en función de la combinación de la distancia al borde y del tipo de vegetación (Fig. 11. Esto apoya la inclusión de estas variables predictoras en los modelos para explicar los patrones espaciales observados de los nidos.

## Modelos nulos
# J2001: modelos nulos con Poisson ajustados por spatstat
J2001_vegdist_50_null <- cbind(J2001_vegdist_50_null, data.frame(BRD_mean = numeric(length = 2000), 
                                    BRD_sd = numeric(length = 2000), 
                                    BRD2_mean = numeric(length = 2000), 
                                    BRD2_sd = numeric(length = 2000), 
                                    BRDHET_mean = numeric(length = 2000), 
                                    BRDHET_sd = numeric(length = 2000), 
                                    BRDHET2_mean = numeric(length = 2000), 
                                    BRDHET2_sd = numeric(length = 2000), 
                                    DIST_mean = numeric(length = 2000), 
                                    DIST_sd = numeric(length = 2000), 
                                    DIST2_mean = numeric(length = 2000), 
                                    DIST2_sd = numeric(length = 2000)))

# la prob. de selección de c/píxel (valor del vector $v) está dada por la predicción del modelo ajustado en un grilla del mísmo número de celdas (50×50)
J2001_vegdist_50_null[c("BRD_mean", "BRD_sd")] <- t(replicate(2000, 
    do.call(function(x) c(mean(x), sd(x)), list(
    sample(J2001_vegdist_50$v, size = sum(!J2001_ants_50$v), replace = FALSE, prob = predict(ModJ2001_AbsDist, ngrid = 50)$v)))))
# el tamaño de muestra es Poisson con media = observados
J2001_vegdist_50_null[c("BRD2_mean", "BRD2_sd")] <- t(replicate(2000, 
    do.call(function(x) c(mean(x), sd(x)), list(
    sample(J2001_vegdist_50$v, size = rpois(1, sum(!J2001_ants_50$v)), replace = TRUE, prob = predict(ModJ2001_AbsDist, ngrid = 50)$v)))))

# idem con modelo ModJ2001_DistxVeg
J2001_vegdist_50_null[c("BRDHET_mean", "BRDHET_sd")] <- t(replicate(2000, 
    do.call(function(x) c(mean(x), sd(x)), list(
    sample(J2001_vegdist_50$v, size = sum(!J2001_ants_50$v), replace = FALSE, prob = predict(ModJ2001_DistxVeg, ngrid = 50)$v)))))

J2001_vegdist_50_null[c("BRDHET2_mean", "BRDHET2_sd")] <- t(replicate(2000, 
    do.call(function(x) c(mean(x), sd(x)), list(
    sample(J2001_vegdist_50$v, size = rpois(1, sum(!J2001_ants_50$v)), replace = TRUE, prob = predict(ModJ2001_DistxVeg, ngrid = 50)$v)))))

# idem con modelo ModJ2001_Dist2
J2001_vegdist_50_null[c("DIST_mean", "DIST_sd")] <- t(replicate(2000, 
    do.call(function(x) c(mean(x), sd(x)), list(
    sample(J2001_vegdist_50$v, size = sum(!J2001_ants_50$v), replace = FALSE, prob = predict(ModJ2001_Dist2, ngrid = 50)$v)))))

J2001_vegdist_50_null[c("DIST2_mean", "DIST2_sd")] <- t(replicate(2000, 
    do.call(function(x) c(mean(x), sd(x)), list(
    sample(J2001_vegdist_50$v, size = rpois(1, sum(!J2001_ants_50$v)), replace = TRUE, prob = predict(ModJ2001_Dist2, ngrid = 50)$v)))))

# incluye a los valores observados del estadístico en la fila 1 de cada vector
J2001_vegdist_50_null[1, c(9, 11, 13, 15, 17, 19)] <- mean(J2001_vegdist_50[!J2001_ants_50$v])
J2001_vegdist_50_null[1, c(10, 12, 14, 16, 18, 20)] <- sd(J2001_vegdist_50[!J2001_ants_50$v])

# J2019: modelos nulos con Poisson ajustados por spatstat
J2019_vegdist_50_null <- cbind(J2019_vegdist_50_null, data.frame(BRD_mean = numeric(length = 2000), 
                                    BRD_sd = numeric(length = 2000), 
                                    BRD2_mean = numeric(length = 2000), 
                                    BRD2_sd = numeric(length = 2000), 
                                    BRDHET_mean = numeric(length = 2000), 
                                    BRDHET_sd = numeric(length = 2000), 
                                    BRDHET2_mean = numeric(length = 2000), 
                                    BRDHET2_sd = numeric(length = 2000), 
                                    DIST_mean = numeric(length = 2000), 
                                    DIST_sd = numeric(length = 2000), 
                                    DIST2_mean = numeric(length = 2000), 
                                    DIST2_sd = numeric(length = 2000)))

# la prob. de selección de c/píxel (valor del vector $v) está dada por la predicción del modelo ajustado en un grilla del mísmo número de celdas (50×50)
J2019_vegdist_50_null[c("BRD_mean", "BRD_sd")] <- t(replicate(2000, 
    do.call(function(x) c(mean(x), sd(x)), list(
    sample(J2019_vegdist_50$v, size = sum(!J2019_ants_50$v), replace = FALSE, prob = predict(ModJ2019_AbsDist, ngrid = 50)$v)))))
# el tamaño de muestra es Poisson con media = observados
J2019_vegdist_50_null[c("BRD2_mean", "BRD2_sd")] <- t(replicate(2000, 
    do.call(function(x) c(mean(x), sd(x)), list(
    sample(J2019_vegdist_50$v, size = rpois(1, sum(!J2019_ants_50$v)), replace = TRUE, prob = predict(ModJ2019_AbsDist, ngrid = 50)$v)))))

# idem con modelo ModJ2019_DistxVeg
J2019_vegdist_50_null[c("BRDHET_mean", "BRDHET_sd")] <- t(replicate(2000, 
    do.call(function(x) c(mean(x), sd(x)), list(
    sample(J2019_vegdist_50$v, size = sum(!J2019_ants_50$v), replace = FALSE, prob = predict(ModJ2019_DistxVeg, ngrid = 50)$v)))))

J2019_vegdist_50_null[c("BRDHET2_mean", "BRDHET2_sd")] <- t(replicate(2000, 
    do.call(function(x) c(mean(x), sd(x)), list(
    sample(J2019_vegdist_50$v, size = rpois(1, sum(!J2019_ants_50$v)), replace = TRUE, prob = predict(ModJ2019_DistxVeg, ngrid = 50)$v)))))

# idem con modelo ModJ2019_Dist2
J2019_vegdist_50_null[c("DIST_mean", "DIST_sd")] <- t(replicate(2000, 
    do.call(function(x) c(mean(x), sd(x)), list(
    sample(J2019_vegdist_50$v, size = sum(!J2019_ants_50$v), replace = FALSE, prob = predict(ModJ2019_Dist2, ngrid = 50)$v)))))

J2019_vegdist_50_null[c("DIST2_mean", "DIST2_sd")] <- t(replicate(2000, 
    do.call(function(x) c(mean(x), sd(x)), list(
    sample(J2019_vegdist_50$v, size = rpois(1, sum(!J2019_ants_50$v)), replace = TRUE, prob = predict(ModJ2019_Dist2, ngrid = 50)$v)))))

# incluye a los valores observados del estadístico en la fila 1 de cada vector
J2019_vegdist_50_null[1, c(9, 11, 13, 15, 17, 19)] <- mean(J2019_vegdist_50[!J2019_ants_50$v])
J2019_vegdist_50_null[1, c(10, 12, 14, 16, 18, 20)] <- sd(J2019_vegdist_50[!J2019_ants_50$v])

## FIGURA
# matriz de 8 subfiguras (una columna para cada año)
par(mar = c(2, 4, 3, 2))
layout(matrix(c(1, 2, 3, 4, 5, 6, 7, 8), 4, 2, byrow = FALSE), heights = c(2, 1.5, 3, 3))

# columna J2001
# 1: mapa de distancias con nidos
plot(J2001_vegdist_50, col = interp.colourmap(colourmap(c("darkgreen", "yellow", "darkred"), inputs = c(-5, 0, 5)), n = 255), ribside = "left", main = "J2001: distancia a bordes")
plot(J2001_ants_50, col = c(rgb(0, 0, 0, .6), rgb(1, 1, 1, 0)), ribbon = F, add = TRUE)

# 2: histogramas como proporciones (con p-valor de test K-S)
# Prueba de K-S
J2001_vegdist_50_KS <- ks.test(J2001_vegdist_50$v, J2001_vegdist_50[as.logical(!J2001_ants_50$v)])
# distribución total
hist(J2001_vegdist_50, breaks = seq(-5.005, 5.005, 1.001), col = "darkgreen", ylab = "Proporción de pixeles", xlab = "Distancia a borde", ylim = c(0, 0.5), probability = TRUE, main = "")
# distribución observada (nidos)
hist(J2001_vegdist_50[as.logical(!J2001_ants_50$v)], breaks = seq(-5.005, 5.005, 1.001), col = rgb(0, 0, 0, .7), border = "black", density = 30, add = TRUE, probability = TRUE)
text(5, 0.5, paste("KS test:\nD =", round(J2001_vegdist_50_KS$stat, digits = 2), "\nP =", round(J2001_vegdist_50_KS$p, digits = 3)), adj = c(1, 1))

# 3: boxplot de la distribución de medias de todos los modelos
boxplot(c(J2001_vegdist_50_null$CSR_mean,
          J2001_vegdist_50_null$CSR2_mean,
          J2001_vegdist_50_null$HET_mean,
          J2001_vegdist_50_null$HET2_mean,
          J2001_vegdist_50_null$BRD_mean,
          J2001_vegdist_50_null$BRD2_mean,
          J2001_vegdist_50_null$BRDHET_mean,
          J2001_vegdist_50_null$BRDHET2_mean,
          J2001_vegdist_50_null$DIST_mean,
          J2001_vegdist_50_null$DIST2_mean)
        ~ factor(c(rep("CSR", 2000), rep("CSR+", 2000), rep("HET", 2000), rep("HET+", 2000), rep("BRD", 2000), rep("BRD+", 2000), rep("BRDHET", 2000), rep("BRDHET+", 2000), rep("DIST2", 2000), rep("DIST2+", 2000)), levels = c("CSR", "CSR+", "HET", "HET+", "BRD", "BRD+", "BRDHET", "BRDHET+", "DIST2", "DIST2+")),
        horizontal = TRUE, ylim = c(-2, 2), xlim = c(10.75, 0.25), las = 1, main = "Media", ann = FALSE)
# valor observado (primer fila en cualquiera de las columnas)
abline(v = J2001_vegdist_50_null$CSR_mean[1], col = "red", lty = 3)
# pseudo-p values from permutations (double the number of more extreme values on the smallest [<0.5] tail)
text(2, 1, paste("P =", 2 * min(mean(J2001_vegdist_50_null$CSR_mean <= J2001_vegdist_50_null$CSR_mean[1]), mean(J2001_vegdist_50_null$CSR_mean >= J2001_vegdist_50_null$CSR_mean[1]))), adj = c(1, 0.5))

text(2, 2, paste("P =", 2 * min(mean(J2001_vegdist_50_null$CSR2_mean <= J2001_vegdist_50_null$CSR2_mean[1]), mean(J2001_vegdist_50_null$CSR2_mean >= J2001_vegdist_50_null$CSR2_mean[1]))), adj = c(1, 0.5))

text(2, 3, paste("P =", 2 * min(mean(J2001_vegdist_50_null$HET_mean <= J2001_vegdist_50_null$HET_mean[1]), mean(J2001_vegdist_50_null$HET_mean >= J2001_vegdist_50_null$HET_mean[1]))), adj = c(1, 0.5))

text(2, 4, paste("P =", 2 * min(mean(J2001_vegdist_50_null$HET2_mean <= J2001_vegdist_50_null$HET2_mean[1]), mean(J2001_vegdist_50_null$HET2_mean >= J2001_vegdist_50_null$HET2_mean[1]))), adj = c(1, 0.5))

text(2, 5, paste("P =", 2 * min(mean(J2001_vegdist_50_null$BRD_mean <= J2001_vegdist_50_null$BRD_mean[1]), mean(J2001_vegdist_50_null$BRD_mean >= J2001_vegdist_50_null$BRD_mean[1]))), adj = c(1, 0.5))

text(2, 6, paste("P =", 2 * min(mean(J2001_vegdist_50_null$BRD2_mean <= J2001_vegdist_50_null$BRD2_mean[1]), mean(J2001_vegdist_50_null$BRD2_mean >= J2001_vegdist_50_null$BRD2_mean[1]))), adj = c(1, 0.5))

text(2, 7, paste("P =", 2 * min(mean(J2001_vegdist_50_null$BRDHET_mean <= J2001_vegdist_50_null$BRDHET_mean[1]), mean(J2001_vegdist_50_null$BRDHET_mean >= J2001_vegdist_50_null$BRDHET_mean[1]))), adj = c(1, 0.5))

text(2, 8, paste("P =", 2 * min(mean(J2001_vegdist_50_null$BRDHET2_mean <= J2001_vegdist_50_null$BRDHET2_mean[1]), mean(J2001_vegdist_50_null$BRDHET2_mean >= J2001_vegdist_50_null$BRDHET2_mean[1]))), adj = c(1, 0.5))

text(2, 9, paste("P =", 2 * min(mean(J2001_vegdist_50_null$DIST_mean <= J2001_vegdist_50_null$DIST_mean[1]), mean(J2001_vegdist_50_null$DIST_mean >= J2001_vegdist_50_null$DIST_mean[1]))), adj = c(1, 0.5))

text(2, 10, paste("P =", 2 * min(mean(J2001_vegdist_50_null$DIST2_mean <= J2001_vegdist_50_null$DIST2_mean[1]), mean(J2001_vegdist_50_null$DIST2_mean >= J2001_vegdist_50_null$DIST2_mean[1]))), adj = c(1, 0.5))

# 4: boxplot de la distribución de sd CSR, CSR2, HET y HET2
boxplot(c(J2001_vegdist_50_null$CSR_sd,
          J2001_vegdist_50_null$CSR2_sd,
          J2001_vegdist_50_null$HET_sd,
          J2001_vegdist_50_null$HET2_sd,
          J2001_vegdist_50_null$BRD_sd,
          J2001_vegdist_50_null$BRD2_sd,
          J2001_vegdist_50_null$BRDHET_sd,
          J2001_vegdist_50_null$BRDHET2_sd,
          J2001_vegdist_50_null$DIST_sd,
          J2001_vegdist_50_null$DIST2_sd)
        ~ factor(c(rep("CSR", 2000), rep("CSR+", 2000), rep("HET", 2000), rep("HET+", 2000), rep("BRD", 2000), rep("BRD+", 2000), rep("BRDHET", 2000), rep("BRDHET+", 2000), rep("DIST2", 2000), rep("DIST2+", 2000)), levels = c("CSR", "CSR+", "HET", "HET+", "BRD", "BRD+", "BRDHET", "BRDHET+", "DIST2", "DIST2+")),
        horizontal = TRUE, ylim = c(0, 3.5), xlim = c(10.75, 0.25), las = 1, main = "Desvío estándar", ann = FALSE)
# valor observado (primer fila en cualquiera de las columnas)
abline(v = J2001_vegdist_50_null$CSR_sd[1], col = "red", lty = 3)
# pseudo-p value from permutations (double the number of more extreme values on the smallest [<0.5] tail)
text(3.5, 1, paste("P =", 2 * min(mean(J2001_vegdist_50_null$CSR_sd <= J2001_vegdist_50_null$CSR_sd[1]), mean(J2001_vegdist_50_null$CSR_sd >= J2001_vegdist_50_null$CSR_sd[1]))), adj = c(1, 0.5))

text(3.5, 2, paste("P =", 2 * min(mean(J2001_vegdist_50_null$CSR2_sd <= J2001_vegdist_50_null$CSR2_sd[1]), mean(J2001_vegdist_50_null$CSR2_sd >= J2001_vegdist_50_null$CSR2_sd[1]))), adj = c(1, 0.5))

text(3.5, 3, paste("P =", 2 * min(mean(J2001_vegdist_50_null$HET_sd <= J2001_vegdist_50_null$HET_sd[1]), mean(J2001_vegdist_50_null$HET_sd >= J2001_vegdist_50_null$HET_sd[1]))), adj = c(1, 0.5))

text(3.5, 4, paste("P =", 2 * min(mean(J2001_vegdist_50_null$HET2_sd <= J2001_vegdist_50_null$HET2_sd[1]), mean(J2001_vegdist_50_null$HET2_sd >= J2001_vegdist_50_null$HET2_sd[1]))), adj = c(1, 0.5))

text(3.5, 5, paste("P =", 2 * min(mean(J2001_vegdist_50_null$BRD_sd <= J2001_vegdist_50_null$BRD_sd[1]), mean(J2001_vegdist_50_null$BRD_sd >= J2001_vegdist_50_null$BRD_sd[1]))), adj = c(1, 0.5))

text(3.5, 6, paste("P =", 2 * min(mean(J2001_vegdist_50_null$BRD2_sd <= J2001_vegdist_50_null$BRD2_sd[1]), mean(J2001_vegdist_50_null$BRD2_sd >= J2001_vegdist_50_null$BRD2_sd[1]))), adj = c(1, 0.5))

text(3.5, 7, paste("P =", 2 * min(mean(J2001_vegdist_50_null$BRDHET_sd <= J2001_vegdist_50_null$BRDHET_sd[1]), mean(J2001_vegdist_50_null$BRDHET_sd >= J2001_vegdist_50_null$BRDHET_sd[1]))), adj = c(1, 0.5))

text(3.5, 8, paste("P =", 2 * min(mean(J2001_vegdist_50_null$BRDHET2_sd <= J2001_vegdist_50_null$BRDHET2_sd[1]), mean(J2001_vegdist_50_null$BRDHET2_sd >= J2001_vegdist_50_null$BRDHET2_sd[1]))), adj = c(1, 0.5))

text(3.5, 9, paste("P =", 2 * min(mean(J2001_vegdist_50_null$DIST_sd <= J2001_vegdist_50_null$DIST_sd[1]), mean(J2001_vegdist_50_null$DIST_sd >= J2001_vegdist_50_null$DIST_sd[1]))), adj = c(1, 0.5))

text(3.5, 10, paste("P =", 2 * min(mean(J2001_vegdist_50_null$DIST2_sd <= J2001_vegdist_50_null$DIST2_sd[1]), mean(J2001_vegdist_50_null$DIST2_sd >= J2001_vegdist_50_null$DIST2_sd[1]))), adj = c(1, 0.5))

# columna J2019
# 5: mapa de distancias con nidos
plot(J2019_vegdist_50, col = interp.colourmap(colourmap(c("darkgreen", "yellow", "darkred"), inputs = c(-5, 0, 5)), n = 255), ribside = "left", main = "J2019: distancia a bordes")
plot(J2019_ants_50, col = c(rgb(0, 0, 0, .6), rgb(1, 1, 1, 0)), ribbon = F, add = TRUE)

# 6: histogramas como proporciones (con p-valor de test K-S)
# Prueba de K-S
J2019_vegdist_50_KS <- ks.test(J2019_vegdist_50$v, J2019_vegdist_50[as.logical(!J2019_ants_50$v)])
# distribución total
hist(J2019_vegdist_50, breaks = seq(-5.005, 5.005, 1.001), col = "darkgreen", ylab = "Proporción de pixeles", xlab = "Distancia a borde", ylim = c(0, 0.5), probability = TRUE, main = "")
# distribución observada (nidos)
hist(J2019_vegdist_50[as.logical(!J2019_ants_50$v)], breaks = seq(-5.005, 5.005, 1.001), col = rgb(0, 0, 0, .7), border = "black", density = 30, add = TRUE, probability = TRUE)
text(5, 0.5, paste("KS test:\nD =", round(J2019_vegdist_50_KS$stat, digits = 2), "\nP =", round(J2019_vegdist_50_KS$p, digits = 3)), adj = c(1, 1))

# 7: boxplot de la distribución de medias de todos los modelos
boxplot(c(J2019_vegdist_50_null$CSR_mean,
          J2019_vegdist_50_null$CSR2_mean,
          J2019_vegdist_50_null$HET_mean,
          J2019_vegdist_50_null$HET2_mean,
          J2019_vegdist_50_null$BRD_mean,
          J2019_vegdist_50_null$BRD2_mean,
          J2019_vegdist_50_null$BRDHET_mean,
          J2019_vegdist_50_null$BRDHET2_mean,
          J2019_vegdist_50_null$DIST_mean,
          J2019_vegdist_50_null$DIST2_mean)
        ~ factor(c(rep("CSR", 2000), rep("CSR+", 2000), rep("HET", 2000), rep("HET+", 2000), rep("BRD", 2000), rep("BRD+", 2000), rep("BRDHET", 2000), rep("BRDHET+", 2000), rep("DIST2", 2000), rep("DIST2+", 2000)), levels = c("CSR", "CSR+", "HET", "HET+", "BRD", "BRD+", "BRDHET", "BRDHET+", "DIST2", "DIST2+")),
        horizontal = TRUE, ylim = c(-2, 2), xlim = c(10.75, 0.25), las = 1, main = "Media", ann = FALSE)
# valor observado (primer fila en cualquiera de las columnas)
abline(v = J2019_vegdist_50_null$CSR_mean[1], col = "red", lty = 3)
# pseudo-p values from permutations (double the number of more extreme values on the smallest [<0.5] tail)
text(2, 1, paste("P =", 2 * min(mean(J2019_vegdist_50_null$CSR_mean <= J2019_vegdist_50_null$CSR_mean[1]), mean(J2019_vegdist_50_null$CSR_mean >= J2019_vegdist_50_null$CSR_mean[1]))), adj = c(1, 0.5))

text(2, 2, paste("P =", 2 * min(mean(J2019_vegdist_50_null$CSR2_mean <= J2019_vegdist_50_null$CSR2_mean[1]), mean(J2019_vegdist_50_null$CSR2_mean >= J2019_vegdist_50_null$CSR2_mean[1]))), adj = c(1, 0.5))

text(2, 3, paste("P =", 2 * min(mean(J2019_vegdist_50_null$HET_mean <= J2019_vegdist_50_null$HET_mean[1]), mean(J2019_vegdist_50_null$HET_mean >= J2019_vegdist_50_null$HET_mean[1]))), adj = c(1, 0.5))

text(2, 4, paste("P =", 2 * min(mean(J2019_vegdist_50_null$HET2_mean <= J2019_vegdist_50_null$HET2_mean[1]), mean(J2019_vegdist_50_null$HET2_mean >= J2019_vegdist_50_null$HET2_mean[1]))), adj = c(1, 0.5))

text(2, 5, paste("P =", 2 * min(mean(J2019_vegdist_50_null$BRD_mean <= J2019_vegdist_50_null$BRD_mean[1]), mean(J2019_vegdist_50_null$BRD_mean >= J2019_vegdist_50_null$BRD_mean[1]))), adj = c(1, 0.5))

text(2, 6, paste("P =", 2 * min(mean(J2019_vegdist_50_null$BRD2_mean <= J2019_vegdist_50_null$BRD2_mean[1]), mean(J2019_vegdist_50_null$BRD2_mean >= J2019_vegdist_50_null$BRD2_mean[1]))), adj = c(1, 0.5))

text(2, 7, paste("P =", 2 * min(mean(J2019_vegdist_50_null$BRDHET_mean <= J2019_vegdist_50_null$BRDHET_mean[1]), mean(J2019_vegdist_50_null$BRDHET_mean >= J2019_vegdist_50_null$BRDHET_mean[1]))), adj = c(1, 0.5))

text(2, 8, paste("P =", 2 * min(mean(J2019_vegdist_50_null$BRDHET2_mean <= J2019_vegdist_50_null$BRDHET2_mean[1]), mean(J2019_vegdist_50_null$BRDHET2_mean >= J2019_vegdist_50_null$BRDHET2_mean[1]))), adj = c(1, 0.5))

text(2, 9, paste("P =", 2 * min(mean(J2019_vegdist_50_null$DIST_mean <= J2019_vegdist_50_null$DIST_mean[1]), mean(J2019_vegdist_50_null$DIST_mean >= J2019_vegdist_50_null$DIST_mean[1]))), adj = c(1, 0.5))

text(2, 10, paste("P =", 2 * min(mean(J2019_vegdist_50_null$DIST2_mean <= J2019_vegdist_50_null$DIST2_mean[1]), mean(J2019_vegdist_50_null$DIST2_mean >= J2019_vegdist_50_null$DIST2_mean[1]))), adj = c(1, 0.5))

# 8: boxplot de la distribución de sd CSR, CSR2, HET y HET2
boxplot(c(J2019_vegdist_50_null$CSR_sd,
          J2019_vegdist_50_null$CSR2_sd,
          J2019_vegdist_50_null$HET_sd,
          J2019_vegdist_50_null$HET2_sd,
          J2019_vegdist_50_null$BRD_sd,
          J2019_vegdist_50_null$BRD2_sd,
          J2019_vegdist_50_null$BRDHET_sd,
          J2019_vegdist_50_null$BRDHET2_sd,
          J2019_vegdist_50_null$DIST_sd,
          J2019_vegdist_50_null$DIST2_sd)
        ~ factor(c(rep("CSR", 2000), rep("CSR+", 2000), rep("HET", 2000), rep("HET+", 2000), rep("BRD", 2000), rep("BRD+", 2000), rep("BRDHET", 2000), rep("BRDHET+", 2000), rep("DIST2", 2000), rep("DIST2+", 2000)), levels = c("CSR", "CSR+", "HET", "HET+", "BRD", "BRD+", "BRDHET", "BRDHET+", "DIST2", "DIST2+")),
        horizontal = TRUE, ylim = c(0, 3.5), xlim = c(10.75, 0.25), las = 1, main = "Desvío estándar", ann = FALSE)
# valor observado (primer fila en cualquiera de las columnas)
abline(v = J2019_vegdist_50_null$CSR_sd[1], col = "red", lty = 3)
# pseudo-p value from permutations (double the number of more extreme values on the smallest [<0.5] tail)
text(3.5, 1, paste("P =", 2 * min(mean(J2019_vegdist_50_null$CSR_sd <= J2019_vegdist_50_null$CSR_sd[1]), mean(J2019_vegdist_50_null$CSR_sd >= J2019_vegdist_50_null$CSR_sd[1]))), adj = c(1, 0.5))

text(3.5, 2, paste("P =", 2 * min(mean(J2019_vegdist_50_null$CSR2_sd <= J2019_vegdist_50_null$CSR2_sd[1]), mean(J2019_vegdist_50_null$CSR2_sd >= J2019_vegdist_50_null$CSR2_sd[1]))), adj = c(1, 0.5))

text(3.5, 3, paste("P =", 2 * min(mean(J2019_vegdist_50_null$HET_sd <= J2019_vegdist_50_null$HET_sd[1]), mean(J2019_vegdist_50_null$HET_sd >= J2019_vegdist_50_null$HET_sd[1]))), adj = c(1, 0.5))

text(3.5, 4, paste("P =", 2 * min(mean(J2019_vegdist_50_null$HET2_sd <= J2019_vegdist_50_null$HET2_sd[1]), mean(J2019_vegdist_50_null$HET2_sd >= J2019_vegdist_50_null$HET2_sd[1]))), adj = c(1, 0.5))

text(3.5, 5, paste("P =", 2 * min(mean(J2019_vegdist_50_null$BRD_sd <= J2019_vegdist_50_null$BRD_sd[1]), mean(J2019_vegdist_50_null$BRD_sd >= J2019_vegdist_50_null$BRD_sd[1]))), adj = c(1, 0.5))

text(3.5, 6, paste("P =", 2 * min(mean(J2019_vegdist_50_null$BRD2_sd <= J2019_vegdist_50_null$BRD2_sd[1]), mean(J2019_vegdist_50_null$BRD2_sd >= J2019_vegdist_50_null$BRD2_sd[1]))), adj = c(1, 0.5))

text(3.5, 7, paste("P =", 2 * min(mean(J2019_vegdist_50_null$BRDHET_sd <= J2019_vegdist_50_null$BRDHET_sd[1]), mean(J2019_vegdist_50_null$BRDHET_sd >= J2019_vegdist_50_null$BRDHET_sd[1]))), adj = c(1, 0.5))

text(3.5, 8, paste("P =", 2 * min(mean(J2019_vegdist_50_null$BRDHET2_sd <= J2019_vegdist_50_null$BRDHET2_sd[1]), mean(J2019_vegdist_50_null$BRDHET2_sd >= J2019_vegdist_50_null$BRDHET2_sd[1]))), adj = c(1, 0.5))

text(3.5, 9, paste("P =", 2 * min(mean(J2019_vegdist_50_null$DIST_sd <= J2019_vegdist_50_null$DIST_sd[1]), mean(J2019_vegdist_50_null$DIST_sd >= J2019_vegdist_50_null$DIST_sd[1]))), adj = c(1, 0.5))

text(3.5, 10, paste("P =", 2 * min(mean(J2019_vegdist_50_null$DIST2_sd <= J2019_vegdist_50_null$DIST2_sd[1]), mean(J2019_vegdist_50_null$DIST2_sd >= J2019_vegdist_50_null$DIST2_sd[1]))), adj = c(1, 0.5))
Idem Fig. \@ref(fig:FiguraDistanciasJ) pero incluyendo modelos nulos que incorporan a la distancia a los bordes entre parches de vegetación [ver texto].

Figura 11: Idem Fig. 6 pero incluyendo modelos nulos que incorporan a la distancia a los bordes entre parches de vegetación [ver texto].

Para el caso de la grilla V en 2019, los ajustes de esos mismos modelos Poisson resultan en intensidades estimadas como:

par(mfrow = c(1, 4), mar = c(2, 2, 1.5, 1))
rangox = seq(min(V2019_vegdist_50), max(V2019_vegdist_50), 0.2)

# Modelo CSR
ModV2019_CSR <- ppm(quadscheme(V2019_ppp, nd = 50) ~ 1)

# Modelo HET categórico: vegetación
ModV2019_Veg <- ppm(quadscheme(V2019_ppp, nd = 50) ~ V2019_veg_50)
barplot(c(exp(coef(ModV2019_Veg)[1]), exp(sum(coef(ModV2019_Veg)))), names.arg = c("veg", "sd"), ylab = "Intensidad (λ)", main = "V2019:HET", col = "white")

# Modelo het continuo: BRD: distancia a borde (independiente del tipo de vegetación)
ModV2019_AbsDist <- ppm(quadscheme(V2019_ppp, nd = 50) ~ abs(V2019_vegdist_50))
plot(effectfun(ModV2019_AbsDist), main = "V2019:BRD")

# Modelos het interacción (distancia y vegetación)
## BRDHET: interacción distancia absoluta / veg
ModV2019_DistxVeg <- ppm(quadscheme(V2019_ppp, nd = 50) ~ abs(V2019_vegdist_50) / V2019_veg_50)
plot(x = rangox[rangox <= 0], y = exp(coef(ModV2019_DistxVeg)[1] - coef(ModV2019_DistxVeg)[2] * rangox[rangox <= 0]), type = "l", ylab = "Intensidad (λ)", main = "V2019:BRDHET", xlim = range(rangox), ylim = c(0, 0.05))
lines(x = rangox[rangox >= 0], y = exp(sum(coef(ModV2019_DistxVeg)[1]) + sum(coef(ModV2019_DistxVeg)[c(2, 3)]) * rangox[rangox >= 0]))

## cuadrática de la distancia con signo (=> máximo a cualquier valor intermedio, no necesariamente cero)
ModV2019_Dist2 <- ppm(quadscheme(V2019_ppp, nd = 50) ~ polynom(V2019_vegdist_50, 2))
plot(effectfun(ModV2019_Dist2), main = "V2019:DIST²")
Idem Fig. \@ref(fig:MasModelosJ2001) para la parcela V durante 2019.

Figura 12: Idem Fig. 9 para la parcela V durante 2019. 

ModsInhomV2019 <- list("V2019:HET" = ModV2019_Veg,
                       "V2019:BRD" = ModV2019_AbsDist,
                       "V2019:BRDHET" = ModV2019_DistxVeg,
                       "V2019:DIST²" = ModV2019_Dist2)

Modelos ajustados: V2019

Modelo Objeto Dependencia Intensidad AIC
CSR ModV2019_CSR — 0.017 429.26
HET ModV2019_Veg veg/sd λveg: 0.015 / λsd: 0.019 430.64
BRD ModV2019_AbsDist abs(dist) 0.041 × 0.511m 425.95
BRDHET ModV2019_DistxVeg abs(dist) * veg/sd λveg: 0.04 × 0.503m / λsd: 0.04 × 0.547m 427.84
DIST2 ModV2019_Dist2 dist2 0.023 × 1.045m × 0.859m2 428.71

En esta grilla, que tiene algunos puntos más alejados de los bordes de la vegetación leñosa, el mejor ajuste corresponde al modelo que considera que la intensidad (densidad de nidos) es mayor en los bordes, con una relación descendiente al aumentar la distancia al borde (Fig. 12). El grado de ajuste es similar (ΔAIC < 2) con el modelo donde esa relación es potencialmente asimétrica hacia uno y otro lado (i.e., hacia el parche de vegetación y hacia el parche de suelo desnudo).

Si extendemos el análisis a los nidos presentes en el mismo lugar en 2001, suponiendo que la vegetación leñosa no se modificó en ese periodo (un supuesto algo riesgoso, ver nota arriba), ninguno de los modelos Poisson de intensidad variable implica un ajuste notoriamente mejor a los demás:

par(mfrow = c(1, 4), mar = c(2, 2, 1.5, 1))

# Modelo CSR
ModV2001_CSR <- ppm(quadscheme(V2001_ppp, nd = 50) ~ 1)

# Modelo HET categórico: vegetación
ModV2001_Veg <- ppm(quadscheme(V2001_ppp, nd = 50) ~ V2019_veg_50)
barplot(c(exp(coef(ModV2001_Veg)[1]), exp(sum(coef(ModV2001_Veg)))), names.arg = c("veg", "sd"), ylab = "Intensidad (λ)", main = "V2001:HET2019", col = "white")

# Modelo HET continuo: distancia a borde (independiente del tipo de vegetación)
ModV2001_AbsDist <- ppm(quadscheme(V2001_ppp, nd = 50) ~ abs(V2019_vegdist_50))
plot(effectfun(ModV2001_AbsDist), main = "V2001:BRD2019")

# Modelos HET interacción (distancia y vegetación)
## interacción distancia absoluta * veg
ModV2001_DistxVeg <- ppm(quadscheme(V2001_ppp, nd = 50) ~ abs(V2019_vegdist_50) / V2019_veg_50)
plot(x = rangox[rangox <= 0], y = exp(coef(ModV2001_DistxVeg)[1] - coef(ModV2001_DistxVeg)[2] * rangox[rangox <= 0]), type = "l", ylab = "Intensidad (λ)", main = "V2001:BRDHET2019", xlim = range(rangox), ylim = c(0, 0.05))
lines(x = rangox[rangox >= 0], y = exp(sum(coef(ModV2001_DistxVeg)[1]) + sum(coef(ModV2019_DistxVeg)[c(2, 3)]) * rangox[rangox >= 0]))

## cuadrática de la distancia con signo (=> máximo a cualquier valor intermedio, no necesariamente cero)
ModV2001_Dist2 <- ppm(quadscheme(V2001_ppp, nd = 50) ~ polynom(V2019_vegdist_50, 2))
plot(effectfun(ModV2001_Dist2), main = "V2001:DIST²2019")
Idem Fig. \@ref(fig:MasModelosJ2001) para la parcela V durante 2001 (usando la información de la vegetación medida en 2019).

Figura 13: Idem Fig. 9 para la parcela V durante 2001 (usando la información de la vegetación medida en 2019). 

ModsInhomV2001 <- list("V2001:HET" = ModV2001_Veg,
                       "V2001:BRD" = ModV2001_AbsDist,
                       "V2001:BRDHET" = ModV2001_DistxVeg,
                       "V2001:DIST²" = ModV2001_Dist2)

Modelos ajustados: V2001 (vegetación 2019)

Modelo Objeto Dependencia Intensidad AIC
CSR ModV2001_CSR — 0.036 773.7
HET ModV2001_Veg veg/sd λveg: 0.031 / λsd: 0.042 773.51
BRD ModV2001_AbsDist abs(dist) 0.045 × 0.842m 774.5
BRDHET ModV2001_DistxVeg abs(dist) * veg/sd λveg: 0.042 × 0.826m / λsd: 0.042 × 0.976m 775.27
DIST2 ModV2001_Dist2 dist2 0.037 × 1.093m × 0.989m2 775.57

Aunque vuelven a aparecer relaciones de intensidad decreciente con la distancia al borde (Fig. (Fig. 13), el incremento en el ajuste (reducción de devianza) no parece justificar su inclusión como variable predictora. Es probable que con las potenciales diferencias de 1-2 m en la vegetación leñosa en casi dos décadas (como las observadas al comparar la vegetación leñosa de la parcela J) las relaciones con los bordes se hayan desdibujado lo suficiente como para borrar la señal que es tan marcada en los otros tres casos.

Aún así, cualquiera de los modelos que incluyen a las distancias al borde como variable predictora acomodan mejor la baja dispersión entre las distancias al borde de los píxeles con nidos en 2019 (Fig. 14), aunque los modelos BRD o BRDHET aparecían como los mejores en términos de reducción de devianza cuando se penaliza el número de parámetros estimados. Del mismo modo, si bien los ajustes por verosimilitud de todos los modelos era similar para los nidos de 2001, la inclusión de variables predictoras asociadas con la vegetación leñosa parece necesaria para acomodar la distribución observada de distancias al borde de la vegetación (Fig. 14).

## Modelos nulos
# V2001: modelos nulos con Poisson ajustados por spatstat
V2001_vegdist_50_null <- cbind(V2001_vegdist_50_null, data.frame(BRD_mean = numeric(length = 2000), 
                                    BRD_sd = numeric(length = 2000), 
                                    BRD2_mean = numeric(length = 2000), 
                                    BRD2_sd = numeric(length = 2000), 
                                    BRDHET_mean = numeric(length = 2000), 
                                    BRDHET_sd = numeric(length = 2000), 
                                    BRDHET2_mean = numeric(length = 2000), 
                                    BRDHET2_sd = numeric(length = 2000), 
                                    DIST_mean = numeric(length = 2000), 
                                    DIST_sd = numeric(length = 2000), 
                                    DIST2_mean = numeric(length = 2000), 
                                    DIST2_sd = numeric(length = 2000)))

# la prob. de selección de c/píxel (valor del vector $v) está dada por la predicción del modelo ajustado en un grilla del mísmo número de celdas (50×50)
V2001_vegdist_50_null[c("BRD_mean", "BRD_sd")] <- t(replicate(2000, 
    do.call(function(x) c(mean(x), sd(x)), list(
    sample(V2019_vegdist_50$v, size = sum(!V2001_ants_50$v), replace = FALSE, prob = predict(ModV2001_AbsDist, ngrid = 50)$v)))))
# el tamaño de muestra es Poisson con media = observados
V2001_vegdist_50_null[c("BRD2_mean", "BRD2_sd")] <- t(replicate(2000, 
    do.call(function(x) c(mean(x), sd(x)), list(
    sample(V2019_vegdist_50$v, size = rpois(1, sum(!V2001_ants_50$v)), replace = TRUE, prob = predict(ModV2001_AbsDist, ngrid = 50)$v)))))

# idem con modelo ModV2001_DistxVeg
V2001_vegdist_50_null[c("BRDHET_mean", "BRDHET_sd")] <- t(replicate(2000, 
    do.call(function(x) c(mean(x), sd(x)), list(
    sample(V2019_vegdist_50$v, size = sum(!V2001_ants_50$v), replace = FALSE, prob = predict(ModV2001_DistxVeg, ngrid = 50)$v)))))

V2001_vegdist_50_null[c("BRDHET2_mean", "BRDHET2_sd")] <- t(replicate(2000, 
    do.call(function(x) c(mean(x), sd(x)), list(
    sample(V2019_vegdist_50$v, size = rpois(1, sum(!V2001_ants_50$v)), replace = TRUE, prob = predict(ModV2001_DistxVeg, ngrid = 50)$v)))))

# idem con modelo ModV2001_Dist2
V2001_vegdist_50_null[c("DIST_mean", "DIST_sd")] <- t(replicate(2000, 
    do.call(function(x) c(mean(x), sd(x)), list(
    sample(V2019_vegdist_50$v, size = sum(!V2001_ants_50$v), replace = FALSE, prob = predict(ModV2001_Dist2, ngrid = 50)$v)))))

V2001_vegdist_50_null[c("DIST2_mean", "DIST2_sd")] <- t(replicate(2000, 
    do.call(function(x) c(mean(x), sd(x)), list(
    sample(V2019_vegdist_50$v, size = rpois(1, sum(!V2001_ants_50$v)), replace = TRUE, prob = predict(ModV2001_Dist2, ngrid = 50)$v)))))

# incluye a los valores observados del estadístico en la fila 1 de cada vector
V2001_vegdist_50_null[1, c(9, 11, 13, 15, 17, 19)] <- mean(V2019_vegdist_50[!V2001_ants_50$v])
V2001_vegdist_50_null[1, c(10, 12, 14, 16, 18, 20)] <- sd(V2019_vegdist_50[!V2001_ants_50$v])

# V2019: modelos nulos con Poisson ajustados por spatstat
V2019_vegdist_50_null <- cbind(V2019_vegdist_50_null, data.frame(BRD_mean = numeric(length = 2000), 
                                    BRD_sd = numeric(length = 2000), 
                                    BRD2_mean = numeric(length = 2000), 
                                    BRD2_sd = numeric(length = 2000), 
                                    BRDHET_mean = numeric(length = 2000), 
                                    BRDHET_sd = numeric(length = 2000), 
                                    BRDHET2_mean = numeric(length = 2000), 
                                    BRDHET2_sd = numeric(length = 2000), 
                                    DIST_mean = numeric(length = 2000), 
                                    DIST_sd = numeric(length = 2000), 
                                    DIST2_mean = numeric(length = 2000), 
                                    DIST2_sd = numeric(length = 2000)))

# la prob. de selección de c/píxel (valor del vector $v) está dada por la predicción del modelo ajustado en un grilla del mísmo número de celdas (50×50)
V2019_vegdist_50_null[c("BRD_mean", "BRD_sd")] <- t(replicate(2000, 
    do.call(function(x) c(mean(x), sd(x)), list(
    sample(V2019_vegdist_50$v, size = sum(!V2019_ants_50$v), replace = FALSE, prob = predict(ModV2019_AbsDist, ngrid = 50)$v)))))
# el tamaño de muestra es Poisson con media = observados
V2019_vegdist_50_null[c("BRD2_mean", "BRD2_sd")] <- t(replicate(2000, 
    do.call(function(x) c(mean(x), sd(x)), list(
    sample(V2019_vegdist_50$v, size = rpois(1, sum(!V2019_ants_50$v)), replace = TRUE, prob = predict(ModV2019_AbsDist, ngrid = 50)$v)))))

# idem con modelo ModV2019_DistxVeg
V2019_vegdist_50_null[c("BRDHET_mean", "BRDHET_sd")] <- t(replicate(2000, 
    do.call(function(x) c(mean(x), sd(x)), list(
    sample(V2019_vegdist_50$v, size = sum(!V2019_ants_50$v), replace = FALSE, prob = predict(ModV2019_DistxVeg, ngrid = 50)$v)))))

V2019_vegdist_50_null[c("BRDHET2_mean", "BRDHET2_sd")] <- t(replicate(2000, 
    do.call(function(x) c(mean(x), sd(x)), list(
    sample(V2019_vegdist_50$v, size = rpois(1, sum(!V2019_ants_50$v)), replace = TRUE, prob = predict(ModV2019_DistxVeg, ngrid = 50)$v)))))

# idem con modelo ModV2019_Dist2
V2019_vegdist_50_null[c("DIST_mean", "DIST_sd")] <- t(replicate(2000, 
    do.call(function(x) c(mean(x), sd(x)), list(
    sample(V2019_vegdist_50$v, size = sum(!V2019_ants_50$v), replace = FALSE, prob = predict(ModV2019_Dist2, ngrid = 50)$v)))))

V2019_vegdist_50_null[c("DIST2_mean", "DIST2_sd")] <- t(replicate(2000, 
    do.call(function(x) c(mean(x), sd(x)), list(
    sample(V2019_vegdist_50$v, size = rpois(1, sum(!V2019_ants_50$v)), replace = TRUE, prob = predict(ModV2019_Dist2, ngrid = 50)$v)))))

# incluye a los valores observados del estadístico en la fila 1 de cada vector
V2019_vegdist_50_null[1, c(9, 11, 13, 15, 17, 19)] <- mean(V2019_vegdist_50[!V2019_ants_50$v])
V2019_vegdist_50_null[1, c(10, 12, 14, 16, 18, 20)] <- sd(V2019_vegdist_50[!V2019_ants_50$v])

## FIGURA
# matriz de 8 subfiguras (una columna para cada año)
par(mar = c(2, 4, 3, 2))
layout(matrix(c(1, 2, 3, 4, 5, 6, 7, 8), 4, 2, byrow = FALSE), heights = c(2, 1.5, 3, 3))

# columna V2001
# 1: mapa de distancias con nidos
plot(V2019_vegdist_50, col = interp.colourmap(colourmap(c("darkgreen", "yellow", "darkred"), inputs = c(-7, 0, 7)), n = 255), ribside = "left", main = "V2001: distancia a bordes")
plot(V2001_ants_50, col = c(rgb(0, 0, 0, .6), rgb(1, 1, 1, 0)), ribbon = F, add = TRUE)

# 2: histogramas como proporciones (con p-valor de test K-S)
# Prueba de K-S
V2001_vegdist_50_KS <- ks.test(V2019_vegdist_50$v, V2019_vegdist_50[as.logical(!V2001_ants_50$v)])
# distribución total
hist(V2019_vegdist_50, breaks = seq(-7.007, 7.007, 1.001), col = "darkgreen", ylab = "Proporción de pixeles", xlab = "Distancia a borde", ylim = c(0, 0.5), probability = TRUE, main = "")
# distribución observada (nidos)
hist(V2019_vegdist_50[as.logical(!V2001_ants_50$v)], breaks = seq(-7.007, 7.007, 1.001), col = rgb(0, 0, 0, .7), border = "black", density = 30, add = TRUE, probability = TRUE)
text(5, 0.5, paste("KS test:\nD =", round(V2001_vegdist_50_KS$stat, digits = 2), "\nP =", round(V2001_vegdist_50_KS$p, digits = 3)), adj = c(1, 1))

# 3: boxplot de la distribución de medias de todos los modelos
boxplot(c(V2001_vegdist_50_null$CSR_mean,
          V2001_vegdist_50_null$CSR2_mean,
          V2001_vegdist_50_null$HET_mean,
          V2001_vegdist_50_null$HET2_mean,
          V2001_vegdist_50_null$BRD_mean,
          V2001_vegdist_50_null$BRD2_mean,
          V2001_vegdist_50_null$BRDHET_mean,
          V2001_vegdist_50_null$BRDHET2_mean,
          V2001_vegdist_50_null$DIST_mean,
          V2001_vegdist_50_null$DIST2_mean)
        ~ factor(c(rep("CSR", 2000), rep("CSR+", 2000), rep("HET", 2000), rep("HET+", 2000), rep("BRD", 2000), rep("BRD+", 2000), rep("BRDHET", 2000), rep("BRDHET+", 2000), rep("DIST2", 2000), rep("DIST2+", 2000)), levels = c("CSR", "CSR+", "HET", "HET+", "BRD", "BRD+", "BRDHET", "BRDHET+", "DIST2", "DIST2+")),
        horizontal = TRUE, ylim = c(-2, 2), xlim = c(10.75, 0.25), las = 1, main = "Media", ann = FALSE)
# valor observado (primer fila en cualquiera de las columnas)
abline(v = V2001_vegdist_50_null$CSR_mean[1], col = "red", lty = 3)
# pseudo-p values from permutations (double the number of more extreme values on the smallest [<0.5] tail)
text(2, 1, paste("P =", 2 * min(mean(V2001_vegdist_50_null$CSR_mean <= V2001_vegdist_50_null$CSR_mean[1]), mean(V2001_vegdist_50_null$CSR_mean >= V2001_vegdist_50_null$CSR_mean[1]))), adj = c(1, 0.5))

text(2, 2, paste("P =", 2 * min(mean(V2001_vegdist_50_null$CSR2_mean <= V2001_vegdist_50_null$CSR2_mean[1]), mean(V2001_vegdist_50_null$CSR2_mean >= V2001_vegdist_50_null$CSR2_mean[1]))), adj = c(1, 0.5))

text(2, 3, paste("P =", 2 * min(mean(V2001_vegdist_50_null$HET_mean <= V2001_vegdist_50_null$HET_mean[1]), mean(V2001_vegdist_50_null$HET_mean >= V2001_vegdist_50_null$HET_mean[1]))), adj = c(1, 0.5))

text(2, 4, paste("P =", 2 * min(mean(V2001_vegdist_50_null$HET2_mean <= V2001_vegdist_50_null$HET2_mean[1]), mean(V2001_vegdist_50_null$HET2_mean >= V2001_vegdist_50_null$HET2_mean[1]))), adj = c(1, 0.5))

text(2, 5, paste("P =", 2 * min(mean(V2001_vegdist_50_null$BRD_mean <= V2001_vegdist_50_null$BRD_mean[1]), mean(V2001_vegdist_50_null$BRD_mean >= V2001_vegdist_50_null$BRD_mean[1]))), adj = c(1, 0.5))

text(2, 6, paste("P =", 2 * min(mean(V2001_vegdist_50_null$BRD2_mean <= V2001_vegdist_50_null$BRD2_mean[1]), mean(V2001_vegdist_50_null$BRD2_mean >= V2001_vegdist_50_null$BRD2_mean[1]))), adj = c(1, 0.5))

text(2, 7, paste("P =", 2 * min(mean(V2001_vegdist_50_null$BRDHET_mean <= V2001_vegdist_50_null$BRDHET_mean[1]), mean(V2001_vegdist_50_null$BRDHET_mean >= V2001_vegdist_50_null$BRDHET_mean[1]))), adj = c(1, 0.5))

text(2, 8, paste("P =", 2 * min(mean(V2001_vegdist_50_null$BRDHET2_mean <= V2001_vegdist_50_null$BRDHET2_mean[1]), mean(V2001_vegdist_50_null$BRDHET2_mean >= V2001_vegdist_50_null$BRDHET2_mean[1]))), adj = c(1, 0.5))

text(2, 9, paste("P =", 2 * min(mean(V2001_vegdist_50_null$DIST_mean <= V2001_vegdist_50_null$DIST_mean[1]), mean(V2001_vegdist_50_null$DIST_mean >= V2001_vegdist_50_null$DIST_mean[1]))), adj = c(1, 0.5))

text(2, 10, paste("P =", 2 * min(mean(V2001_vegdist_50_null$DIST2_mean <= V2001_vegdist_50_null$DIST2_mean[1]), mean(V2001_vegdist_50_null$DIST2_mean >= V2001_vegdist_50_null$DIST2_mean[1]))), adj = c(1, 0.5))

# 4: boxplot de la distribución de sd CSR, CSR2, HET y HET2
boxplot(c(V2001_vegdist_50_null$CSR_sd,
          V2001_vegdist_50_null$CSR2_sd,
          V2001_vegdist_50_null$HET_sd,
          V2001_vegdist_50_null$HET2_sd,
          V2001_vegdist_50_null$BRD_sd,
          V2001_vegdist_50_null$BRD2_sd,
          V2001_vegdist_50_null$BRDHET_sd,
          V2001_vegdist_50_null$BRDHET2_sd,
          V2001_vegdist_50_null$DIST_sd,
          V2001_vegdist_50_null$DIST2_sd)
        ~ factor(c(rep("CSR", 2000), rep("CSR+", 2000), rep("HET", 2000), rep("HET+", 2000), rep("BRD", 2000), rep("BRD+", 2000), rep("BRDHET", 2000), rep("BRDHET+", 2000), rep("DIST2", 2000), rep("DIST2+", 2000)), levels = c("CSR", "CSR+", "HET", "HET+", "BRD", "BRD+", "BRDHET", "BRDHET+", "DIST2", "DIST2+")),
        horizontal = TRUE, ylim = c(0, 3.5), xlim = c(10.75, 0.25), las = 1, main = "Desvío estándar", ann = FALSE)
# valor observado (primer fila en cualquiera de las columnas)
abline(v = V2001_vegdist_50_null$CSR_sd[1], col = "red", lty = 3)
# pseudo-p value from permutations (double the number of more extreme values on the smallest [<0.5] tail)
text(3.5, 1, paste("P =", 2 * min(mean(V2001_vegdist_50_null$CSR_sd <= V2001_vegdist_50_null$CSR_sd[1]), mean(V2001_vegdist_50_null$CSR_sd >= V2001_vegdist_50_null$CSR_sd[1]))), adj = c(1, 0.5))

text(3.5, 2, paste("P =", 2 * min(mean(V2001_vegdist_50_null$CSR2_sd <= V2001_vegdist_50_null$CSR2_sd[1]), mean(V2001_vegdist_50_null$CSR2_sd >= V2001_vegdist_50_null$CSR2_sd[1]))), adj = c(1, 0.5))

text(3.5, 3, paste("P =", 2 * min(mean(V2001_vegdist_50_null$HET_sd <= V2001_vegdist_50_null$HET_sd[1]), mean(V2001_vegdist_50_null$HET_sd >= V2001_vegdist_50_null$HET_sd[1]))), adj = c(1, 0.5))

text(3.5, 4, paste("P =", 2 * min(mean(V2001_vegdist_50_null$HET2_sd <= V2001_vegdist_50_null$HET2_sd[1]), mean(V2001_vegdist_50_null$HET2_sd >= V2001_vegdist_50_null$HET2_sd[1]))), adj = c(1, 0.5))

text(3.5, 5, paste("P =", 2 * min(mean(V2001_vegdist_50_null$BRD_sd <= V2001_vegdist_50_null$BRD_sd[1]), mean(V2001_vegdist_50_null$BRD_sd >= V2001_vegdist_50_null$BRD_sd[1]))), adj = c(1, 0.5))

text(3.5, 6, paste("P =", 2 * min(mean(V2001_vegdist_50_null$BRD2_sd <= V2001_vegdist_50_null$BRD2_sd[1]), mean(V2001_vegdist_50_null$BRD2_sd >= V2001_vegdist_50_null$BRD2_sd[1]))), adj = c(1, 0.5))

text(3.5, 7, paste("P =", 2 * min(mean(V2001_vegdist_50_null$BRDHET_sd <= V2001_vegdist_50_null$BRDHET_sd[1]), mean(V2001_vegdist_50_null$BRDHET_sd >= V2001_vegdist_50_null$BRDHET_sd[1]))), adj = c(1, 0.5))

text(3.5, 8, paste("P =", 2 * min(mean(V2001_vegdist_50_null$BRDHET2_sd <= V2001_vegdist_50_null$BRDHET2_sd[1]), mean(V2001_vegdist_50_null$BRDHET2_sd >= V2001_vegdist_50_null$BRDHET2_sd[1]))), adj = c(1, 0.5))

text(3.5, 9, paste("P =", 2 * min(mean(V2001_vegdist_50_null$DIST_sd <= V2001_vegdist_50_null$DIST_sd[1]), mean(V2001_vegdist_50_null$DIST_sd >= V2001_vegdist_50_null$DIST_sd[1]))), adj = c(1, 0.5))

text(3.5, 10, paste("P =", 2 * min(mean(V2001_vegdist_50_null$DIST2_sd <= V2001_vegdist_50_null$DIST2_sd[1]), mean(V2001_vegdist_50_null$DIST2_sd >= V2001_vegdist_50_null$DIST2_sd[1]))), adj = c(1, 0.5))

# columna V2019
# 5: mapa de distancias con nidos
plot(V2019_vegdist_50, col = interp.colourmap(colourmap(c("darkgreen", "yellow", "darkred"), inputs = c(-7, 0, 7)), n = 255), ribside = "left", main = "V2019: distancia a bordes")
plot(V2019_ants_50, col = c(rgb(0, 0, 0, .6), rgb(1, 1, 1, 0)), ribbon = F, add = TRUE)

# 6: histogramas como proporciones (con p-valor de test K-S)
# Prueba de K-S
V2019_vegdist_50_KS <- ks.test(V2019_vegdist_50$v, V2019_vegdist_50[as.logical(!V2019_ants_50$v)])
# distribución total
hist(V2019_vegdist_50, breaks = seq(-7.007, 7.007, 1.001), col = "darkgreen", ylab = "Proporción de pixeles", xlab = "Distancia a borde", ylim = c(0, 0.5), probability = TRUE, main = "")
# distribución observada (nidos)
hist(V2019_vegdist_50[as.logical(!V2019_ants_50$v)], breaks = seq(-7.007, 7.007, 1.001), col = rgb(0, 0, 0, .7), border = "black", density = 30, add = TRUE, probability = TRUE)
text(5, 0.5, paste("KS test:\nD =", round(V2019_vegdist_50_KS$stat, digits = 2), "\nP =", round(V2019_vegdist_50_KS$p, digits = 3)), adj = c(1, 1))

# 7: boxplot de la distribución de medias de todos los modelos
boxplot(c(V2019_vegdist_50_null$CSR_mean,
          V2019_vegdist_50_null$CSR2_mean,
          V2019_vegdist_50_null$HET_mean,
          V2019_vegdist_50_null$HET2_mean,
          V2019_vegdist_50_null$BRD_mean,
          V2019_vegdist_50_null$BRD2_mean,
          V2019_vegdist_50_null$BRDHET_mean,
          V2019_vegdist_50_null$BRDHET2_mean,
          V2019_vegdist_50_null$DIST_mean,
          V2019_vegdist_50_null$DIST2_mean)
        ~ factor(c(rep("CSR", 2000), rep("CSR+", 2000), rep("HET", 2000), rep("HET+", 2000), rep("BRD", 2000), rep("BRD+", 2000), rep("BRDHET", 2000), rep("BRDHET+", 2000), rep("DIST2", 2000), rep("DIST2+", 2000)), levels = c("CSR", "CSR+", "HET", "HET+", "BRD", "BRD+", "BRDHET", "BRDHET+", "DIST2", "DIST2+")),
        horizontal = TRUE, ylim = c(-2, 2), xlim = c(10.75, 0.25), las = 1, main = "Media", ann = FALSE)
# valor observado (primer fila en cualquiera de las columnas)
abline(v = V2019_vegdist_50_null$CSR_mean[1], col = "red", lty = 3)
# pseudo-p values from permutations (double the number of more extreme values on the smallest [<0.5] tail)
text(2, 1, paste("P =", 2 * min(mean(V2019_vegdist_50_null$CSR_mean <= V2019_vegdist_50_null$CSR_mean[1]), mean(V2019_vegdist_50_null$CSR_mean >= V2019_vegdist_50_null$CSR_mean[1]))), adj = c(1, 0.5))

text(2, 2, paste("P =", 2 * min(mean(V2019_vegdist_50_null$CSR2_mean <= V2019_vegdist_50_null$CSR2_mean[1]), mean(V2019_vegdist_50_null$CSR2_mean >= V2019_vegdist_50_null$CSR2_mean[1]))), adj = c(1, 0.5))

text(2, 3, paste("P =", 2 * min(mean(V2019_vegdist_50_null$HET_mean <= V2019_vegdist_50_null$HET_mean[1]), mean(V2019_vegdist_50_null$HET_mean >= V2019_vegdist_50_null$HET_mean[1]))), adj = c(1, 0.5))

text(2, 4, paste("P =", 2 * min(mean(V2019_vegdist_50_null$HET2_mean <= V2019_vegdist_50_null$HET2_mean[1]), mean(V2019_vegdist_50_null$HET2_mean >= V2019_vegdist_50_null$HET2_mean[1]))), adj = c(1, 0.5))

text(2, 5, paste("P =", 2 * min(mean(V2019_vegdist_50_null$BRD_mean <= V2019_vegdist_50_null$BRD_mean[1]), mean(V2019_vegdist_50_null$BRD_mean >= V2019_vegdist_50_null$BRD_mean[1]))), adj = c(1, 0.5))

text(2, 6, paste("P =", 2 * min(mean(V2019_vegdist_50_null$BRD2_mean <= V2019_vegdist_50_null$BRD2_mean[1]), mean(V2019_vegdist_50_null$BRD2_mean >= V2019_vegdist_50_null$BRD2_mean[1]))), adj = c(1, 0.5))

text(2, 7, paste("P =", 2 * min(mean(V2019_vegdist_50_null$BRDHET_mean <= V2019_vegdist_50_null$BRDHET_mean[1]), mean(V2019_vegdist_50_null$BRDHET_mean >= V2019_vegdist_50_null$BRDHET_mean[1]))), adj = c(1, 0.5))

text(2, 8, paste("P =", 2 * min(mean(V2019_vegdist_50_null$BRDHET2_mean <= V2019_vegdist_50_null$BRDHET2_mean[1]), mean(V2019_vegdist_50_null$BRDHET2_mean >= V2019_vegdist_50_null$BRDHET2_mean[1]))), adj = c(1, 0.5))

text(2, 9, paste("P =", 2 * min(mean(V2019_vegdist_50_null$DIST_mean <= V2019_vegdist_50_null$DIST_mean[1]), mean(V2019_vegdist_50_null$DIST_mean >= V2019_vegdist_50_null$DIST_mean[1]))), adj = c(1, 0.5))

text(2, 10, paste("P =", 2 * min(mean(V2019_vegdist_50_null$DIST2_mean <= V2019_vegdist_50_null$DIST2_mean[1]), mean(V2019_vegdist_50_null$DIST2_mean >= V2019_vegdist_50_null$DIST2_mean[1]))), adj = c(1, 0.5))

# 8: boxplot de la distribución de sd CSR, CSR2, HET y HET2
boxplot(c(V2019_vegdist_50_null$CSR_sd,
          V2019_vegdist_50_null$CSR2_sd,
          V2019_vegdist_50_null$HET_sd,
          V2019_vegdist_50_null$HET2_sd,
          V2019_vegdist_50_null$BRD_sd,
          V2019_vegdist_50_null$BRD2_sd,
          V2019_vegdist_50_null$BRDHET_sd,
          V2019_vegdist_50_null$BRDHET2_sd,
          V2019_vegdist_50_null$DIST_sd,
          V2019_vegdist_50_null$DIST2_sd)
        ~ factor(c(rep("CSR", 2000), rep("CSR+", 2000), rep("HET", 2000), rep("HET+", 2000), rep("BRD", 2000), rep("BRD+", 2000), rep("BRDHET", 2000), rep("BRDHET+", 2000), rep("DIST2", 2000), rep("DIST2+", 2000)), levels = c("CSR", "CSR+", "HET", "HET+", "BRD", "BRD+", "BRDHET", "BRDHET+", "DIST2", "DIST2+")),
        horizontal = TRUE, ylim = c(0, 3.5), xlim = c(10.75, 0.25), las = 1, main = "Desvío estándar", ann = FALSE)
# valor observado (primer fila en cualquiera de las columnas)
abline(v = V2019_vegdist_50_null$CSR_sd[1], col = "red", lty = 3)
# pseudo-p value from permutations (double the number of more extreme values on the smallest [<0.5] tail)
text(3.5, 1, paste("P =", 2 * min(mean(V2019_vegdist_50_null$CSR_sd <= V2019_vegdist_50_null$CSR_sd[1]), mean(V2019_vegdist_50_null$CSR_sd >= V2019_vegdist_50_null$CSR_sd[1]))), adj = c(1, 0.5))

text(3.5, 2, paste("P =", 2 * min(mean(V2019_vegdist_50_null$CSR2_sd <= V2019_vegdist_50_null$CSR2_sd[1]), mean(V2019_vegdist_50_null$CSR2_sd >= V2019_vegdist_50_null$CSR2_sd[1]))), adj = c(1, 0.5))

text(3.5, 3, paste("P =", 2 * min(mean(V2019_vegdist_50_null$HET_sd <= V2019_vegdist_50_null$HET_sd[1]), mean(V2019_vegdist_50_null$HET_sd >= V2019_vegdist_50_null$HET_sd[1]))), adj = c(1, 0.5))

text(3.5, 4, paste("P =", 2 * min(mean(V2019_vegdist_50_null$HET2_sd <= V2019_vegdist_50_null$HET2_sd[1]), mean(V2019_vegdist_50_null$HET2_sd >= V2019_vegdist_50_null$HET2_sd[1]))), adj = c(1, 0.5))

text(3.5, 5, paste("P =", 2 * min(mean(V2019_vegdist_50_null$BRD_sd <= V2019_vegdist_50_null$BRD_sd[1]), mean(V2019_vegdist_50_null$BRD_sd >= V2019_vegdist_50_null$BRD_sd[1]))), adj = c(1, 0.5))

text(3.5, 6, paste("P =", 2 * min(mean(V2019_vegdist_50_null$BRD2_sd <= V2019_vegdist_50_null$BRD2_sd[1]), mean(V2019_vegdist_50_null$BRD2_sd >= V2019_vegdist_50_null$BRD2_sd[1]))), adj = c(1, 0.5))

text(3.5, 7, paste("P =", 2 * min(mean(V2019_vegdist_50_null$BRDHET_sd <= V2019_vegdist_50_null$BRDHET_sd[1]), mean(V2019_vegdist_50_null$BRDHET_sd >= V2019_vegdist_50_null$BRDHET_sd[1]))), adj = c(1, 0.5))

text(3.5, 8, paste("P =", 2 * min(mean(V2019_vegdist_50_null$BRDHET2_sd <= V2019_vegdist_50_null$BRDHET2_sd[1]), mean(V2019_vegdist_50_null$BRDHET2_sd >= V2019_vegdist_50_null$BRDHET2_sd[1]))), adj = c(1, 0.5))

text(3.5, 9, paste("P =", 2 * min(mean(V2019_vegdist_50_null$DIST_sd <= V2019_vegdist_50_null$DIST_sd[1]), mean(V2019_vegdist_50_null$DIST_sd >= V2019_vegdist_50_null$DIST_sd[1]))), adj = c(1, 0.5))

text(3.5, 10, paste("P =", 2 * min(mean(V2019_vegdist_50_null$DIST2_sd <= V2019_vegdist_50_null$DIST2_sd[1]), mean(V2019_vegdist_50_null$DIST2_sd >= V2019_vegdist_50_null$DIST2_sd[1]))), adj = c(1, 0.5))
Idem Fig. \@ref(fig:FiguraPoissonJ) pero para la parcela V.

Figura 14: Idem Fig. 11 pero para la parcela V.

Síntesis: modelos Poisson con parcelas como réplicas

Dado que ambas parcelas fueron establecidas arbitrariamente en zonas cercanas de un mismo algarrobal y fueron muestreadas en los mismos momentos, pueden ser consideradas réplicas de un mismo proceso espacial. En ese caso, los patrones espaciales de puntos resultan agregados en 2001 (densidad alta), y con variaciones importantes a diferentes escalas pero aún sin excederse de lo esperado para un proceso aleatorio en 2019 (densidad baja; Fig. 15).

par(mfrow = c(1, 2), mar = c(2, 2, 1.5, 1))
plot(pool(envelope(J2001_ppp, Lest, nsim = 399, nrank = 10, simulate = expression(pixelcentres(owin(c(0.5, 50.5), c(0.5, 50.5), mask = matrix(sample(J2001_antmask_50$m), 50, 50), unitname = "m"))), funargs = list(r = 0:12, rmax = 12), global = FALSE, verbose = FALSE, savefuns = TRUE), envelope(V2001_ppp, Lest, nsim = 399, nrank = 10, simulate = expression(pixelcentres(owin(c(0.5, 50.5), c(0.5, 50.5), mask = matrix(sample(V2001_antmask_50$m), 50, 50), unitname = "m"))), funargs = list(r = 0:12, rmax = 12), global = FALSE, verbose = FALSE, savefuns = TRUE)), . - r ~ r, legend = FALSE, main = "2001")

plot(pool(envelope(J2019_ppp, Lest, nsim = 399, nrank = 10, simulate = expression(pixelcentres(owin(c(0.5, 50.5), c(0.5, 50.5), mask = matrix(sample(J2019_antmask_50$m), 50, 50), unitname = "m"))), funargs = list(r = 0:12, rmax = 12), global = FALSE, verbose = FALSE, savefuns = TRUE), envelope(V2019_ppp, Lest, nsim = 399, nrank = 10, simulate = expression(pixelcentres(owin(c(0.5, 50.5), c(0.5, 50.5), mask = matrix(sample(V2019_antmask_50$m), 50, 50), unitname = "m"))), funargs = list(r = 0:12, rmax = 12), global = FALSE, verbose = FALSE, savefuns = TRUE)), . - r ~ r, legend = FALSE, main = "2019")
Autocorrelación espacial según la función L de Ripley para las dos parcelas, consideradas réplicas de un mismo proceso espacial que determina la ubicación de las colonias de *P. bergi*, en 2001 y en 2019.

Figura 15: Autocorrelación espacial según la función L de Ripley para las dos parcelas, consideradas réplicas de un mismo proceso espacial que determina la ubicación de las colonias de P. bergi, en 2001 y en 2019.

El ajuste de los mismos modelos Poisson ya descriptos, pero ahora tomando como información ambas parcelas en cada año (y sujeto a la limitación ya comentada de que la información de la vegetación en V solo fue medida en 2019) produce los siguientes estimadores:

# todas las listas de spatstat (solist) relacionadas dentro de hyperframes  para los análisis como réplicas (incluyendo a los quad a partir de los ppp) [nota: para 2001 uso vegetación de V2019]
X2001_HF <- hyperframe(Nidos = solist(J2001_ppp,V2001_ppp),
                       Veg = solist(J2001_veg_50,V2019_veg_50),
                       Dist = solist(J2001_vegdist_50,V2019_vegdist_50),
                       Pix = solist(quadscheme(J2001_ppp, nd = 50), quadscheme(V2001_ppp, nd = 50)))
X2019_HF <- hyperframe(Nidos = solist(J2019_ppp,V2019_ppp),
                       Veg = solist(J2019_veg_50,V2019_veg_50),
                       Dist = solist(J2019_vegdist_50,V2019_vegdist_50), 
                       Pix = solist(quadscheme(J2019_ppp, nd = 50), quadscheme(V2019_ppp, nd = 50)))
Full_HF <- hyperframe(Nidos = solist(J2001_ppp, V2001_ppp, J2019_ppp, V2019_ppp),
                      Veg = solist(J2001_veg_50, V2019_veg_50, J2019_veg_50, V2019_veg_50),
                      Dist = solist(J2001_vegdist_50, V2019_vegdist_50, J2019_vegdist_50, V2019_vegdist_50), 
                      Pix = solist(quadscheme(J2001_ppp, nd = 50), quadscheme(V2001_ppp, nd = 50), quadscheme(J2019_ppp, nd = 50), quadscheme(V2019_ppp, nd = 50)),
                      Grid = as.factor(c("J", "V", "J", "V")),
                      Year = as.factor(c("2001", "2001", "2019", "2019")))

## 2001
# Modelo CSR
ModX2001_CSR <- mppm(Pix ~ 1, X2001_HF)

# Modelo HET categórico: vegetación
ModX2001_Veg <- mppm(Pix ~ Veg, X2001_HF)

# Modelo het continuo: BRD distancia a borde (independiente del tipo de vegetación)
ModX2001_AbsDist <- mppm(Pix ~ abs(Dist), X2001_HF)

# Modelos het interacción (distancia y vegetación)
## BRDHET: interacción distancia absoluta * veg
ModX2001_DistxVeg <- mppm(Pix ~ abs(Dist) / Veg, X2001_HF)

## DIST²: cuadrática de la distancia con signo
ModX2001_Dist2 <- mppm(Pix ~ polynom(Dist, 2), X2001_HF)

#figuras intensidad predicha por modelos
par(mfrow = c(1, 4), mar = c(2, 2, 1.5, 1))
rangox= seq(min(X2001_HF$Dist), max(X2001_HF$Dist), 0.2)

barplot(c(exp(coef(ModX2001_Veg)[1]), exp(sum(coef(ModX2001_Veg)))), names.arg = c("veg", "sd"), ylab = "Intensidad (λ)", main = "2001:HET", col = "white")

plot(x = rangox, y = c(exp(coef(ModX2001_AbsDist)[1] - coef(ModX2001_AbsDist)[2] * rangox[rangox < 0]), exp(coef(ModX2001_AbsDist)[1] + coef(ModX2001_AbsDist)[2] * rangox[rangox > 0])), type = "l", ylab = "Intensidad (λ)", main = "2001:BRD")

plot(x = rangox[rangox <= 0], y = exp(coef(ModX2001_DistxVeg)[1] - coef(ModX2001_DistxVeg)[2] * rangox[rangox <= 0]), type = "l", ylab = "Intensidad (λ)", main = "2001:BRDHET", xlim = range(rangox), ylim = c(0, 0.06))
lines(x = rangox[rangox >= 0], y = exp(sum(coef(ModX2001_DistxVeg)[1]) + sum(coef(ModX2001_DistxVeg)[c(2, 3)]) * rangox[rangox >= 0]))

plot(x = rangox, y = exp(coef(ModX2001_Dist2)[1] + coef(ModX2001_Dist2)[2] * rangox + coef(ModX2001_Dist2)[3] * rangox ^ 2), type = "l", ylab = "Intensidad (λ)", main = "2001:DIST²")
Idem Fig. \@ref(fig:MasModelosJ2001) pero para ambas parcelas, consideradas réplicas, durante 2001.

Figura 16: Idem Fig. 9 pero para ambas parcelas, consideradas réplicas, durante 2001.

Modelos ajustados: réplicas 2001

Modelo Objeto Dependencia Intensidad AIC
CSR ModX2001_CSR — 0.033 1457.71
HET ModX2001_Veg veg/sd λveg: 0.024 / λsd: 0.047 1441.26
BRD ModX2001_AbsDist abs(dist) 0.052 × 0.714m 1453.77
BRDHET ModX2001_DistxVeg abs(dist) * veg/sd λveg: 0.044 × 0.659m / λsd: 0.044 × 1.028m 1441.24
DIST2 ModX2001_Dist2 dist2 0.036 × 1.252m × 0.972m2 1442.34 
## 2019
# Modelo CSR
ModX2019_CSR <- mppm(Pix ~ 1, X2019_HF)

# Modelo HET categórico: vegetación
ModX2019_Veg <- mppm(Pix ~ Veg, X2019_HF)

# Modelo het continuo: BRD: distancia a borde (independiente del tipo de vegetación)
ModX2019_AbsDist <- mppm(Pix ~ abs(Dist), X2019_HF)

# Modelos het interacción (distancia y vegetación)
## interacción distancia absoluta * veg: BRDHET 
ModX2019_DistxVeg <- mppm(Pix ~ abs(Dist) / Veg, X2019_HF)

## DIST²: cuadrática de la distancia con signo
ModX2019_Dist2 <- mppm(Pix ~ polynom(Dist, 2), X2019_HF)

#figuras intensidad predicha por modelos
par(mfrow = c(1, 4), mar = c(2, 2, 1.5, 1))
rangox= seq(min(X2019_HF$Dist), max(X2019_HF$Dist), 0.2)

barplot(c(exp(coef(ModX2019_Veg)[1]), exp(sum(coef(ModX2019_Veg)))), names.arg = c("veg", "sd"), ylab = "Intensidad (λ)", main = "2019:HET", col = "white")

plot(x = rangox, y = c(exp(coef(ModX2019_AbsDist)[1] - coef(ModX2019_AbsDist)[2] * rangox[rangox < 0]), exp(coef(ModX2019_AbsDist)[1] + coef(ModX2019_AbsDist)[2] * rangox[rangox > 0])), type = "l", ylab = "Intensidad (λ)", main = "2019:BRD")

plot(x = rangox[rangox <= 0], y = exp(coef(ModX2019_DistxVeg)[1] - coef(ModX2019_DistxVeg)[2] * rangox[rangox <= 0]), type = "l", ylab = "Intensidad (λ)", main = "2019:BRDHET", xlim = range(rangox), ylim = c(0, 0.05))
lines(x = rangox[rangox >= 0], y = exp(sum(coef(ModX2019_DistxVeg)[1]) + sum(coef(ModX2019_DistxVeg)[c(2, 3)]) * rangox[rangox >= 0]))

plot(x = rangox, y = exp(coef(ModX2019_Dist2)[1] + coef(ModX2019_Dist2)[2] * rangox + coef(ModX2019_Dist2)[3] * rangox ^ 2), type = "l", ylab = "Intensidad (λ)", main = "2019:DIST²")
Idem Fig. \@ref(fig:MasModelosJ2001) pero para ambas parcelas, consideradas réplicas, durante 2019.

Figura 17: Idem Fig. 9 pero para ambas parcelas, consideradas réplicas, durante 2019.

Modelos ajustados: réplicas 2019

Modelo Objeto Dependencia Intensidad AIC
CSR ModX2019_CSR — 0.015 773.54
HET ModX2019_Veg veg/sd λveg: 0.012 / λsd: 0.02 769.87
BRD ModX2019_AbsDist abs(dist) 0.042 × 0.461m 762.18
BRDHET ModX2019_DistxVeg abs(dist) * veg/sd λveg: 0.036 × 0.447m / λsd: 0.036 × 0.635m 760.63
DIST2 ModX2019_Dist2 dist2 0.022 × 1.193m × 0.851m2 761.55

Nuevamente, los modelos Poisson heterogéneos que incorporan la información de la disposición de la vegetación leñosa ajustan mejor a los patrones de puntos observados, y las relaciones estimadas son relativamente consistentes.

Modelos con interacción entre puntos (Gibbs)

Como vimos, en todos los casos (dos sitios × dos épocas) los modelos Poisson que incorporan al tipo y distancia a la vegetación mejoran sustancialmente la descripción de los patrones de nidos respecto a los esperados por patrones completamente aleatorios (= Poisson homogéneos). Sin embargo, dado que ni las distancias ni los tipos de vegetación muestran patrones espaciales particulares dentro las parcelas, estos modelos no parecen ser suficientes para explicar el agregamiento o la sobredispersión observados en algunos casos.

Por ejemplo, el agregamiento espacial de los píxeles con nidos observado en J2001 (con la densidad más alta) parece exceder los límites de lo esperado por cada uno de estos modelos tomados como modelos nulos, en particular a distancias entre nidos de 4-6 m. Los p-valores de las pruebas de hipótesis globales basadas en simulaciones (MAD test en negro, Montecarlo de dos etapas en rojo) son siempre relativamente bajos, y en ocasiones menores a p = 0.05 (i.e., el patrón observado muestra valores de autocorrelación espacial más extremos que lo esperado por el modelo Poisson heterogéneo propuesto; Fig. 18).

En el extremo opuesto, en V2019, con la densidad de nidos más baja, los modelos basados en la relación entre los nidos y la vegetación tampoco parecen ser suficientes para acomodar los patrones de sobredispersión a distancias cortas (menores a 6-7 m) detectados por las funciones de autocorrelación espacial (Fig. 18).

Esta evidencia sugiere que el supuesto de independencia entre puntos de los modelos Poisson heterogéneos ajustados (i.e., que más allá de que la probabilidad de ocurrencia de nidos varíe con la presencia/distancia a la vegetación leñosa, la presencia de un nido no afecte la probabilidad de que haya otro a determinada distancia) no sería válido.

par(mfcol = c(2, 4), mar = c(2, 2, 1.5, 1))
# genera cada patrón de puntos en cada simulación (simulate =) con un muestreo aleatorio sin repetición del mismo tamaño del observado (= modelos "estrictos") a partir del vector de probabilidades de selección (intensidad) predicho por cada modelo Poisson ajustado
# los ppp respetan las coordenadas enteras (centros de pixeles) y el tamaño de ventana: owin(c(0.5, 50.5), c(0.5, 50.5))
# [si fuesen modelos sin restricciones espaciales sería directamente envelope(Modelo, …) pudiendo usar fix.n = TRUE para que el MonteCarlo sea válido (i.e., no excesivamente conservativo)]

#J2001
for(i in names(ModsInhomJ2001)) {
  temp <- envelope(J2001_ppp, fun = Linhom, nsim = 399, nrank = 10, simulate = expression(as.ppp(arrayInd(sample(2500, prob = predict(ModsInhomJ2001[[i]], ngrid = 50)$v, size = sum(!J2001_ants_50$v)), .dim = c(50, 50)), W = owin(c(0.5, 50.5), c(0.5, 50.5)))), funargs = list(r = 0:12, rmax = 12), global = FALSE, verbose = FALSE, savefuns = TRUE, savepatterns = TRUE)
  plot(temp, . - mmean ~ r, legend = FALSE, main = paste(i, "- L(inhom)"), cex.axis = 0.8)
  text(x = 12, y = 0.5, paste("p = ", round(mad.test(temp, use.theory = FALSE)$p.value, 3)), adj = c(1, 0))
  text(x = 12, y = 0.5, paste("p = ", round(bits.test(ModsInhomJ2001[[i]], fun = Linhom, exponent = Inf, rinterval = c(0, 12), nsim = 39, leaveoneout = FALSE, use.theory = FALSE, verbose = FALSE)$p.value, 3)), adj = c(1, 1.5), col = "red")
  
  temp2 <- envelope(J2001_ppp, fun = pcfinhom, nsim = 399, nrank = 10, r = 0:12, rmax = 12, verbose = FALSE, simulate = temp, savefuns = TRUE)
  plot(temp2, . - mmean + 1 ~ r, legend = FALSE, main = paste(i, "- pcf(inhom)"), cex.axis = 0.8)
  text(x = 12, y = 1.6, paste("p = ", round(mad.test(temp2, use.theory = FALSE)$p.value, 3)), adj = c(1, 0))
  text(x = 12, y = 1.6, paste("p = ", round(bits.test(ModsInhomJ2001[[i]], fun = pcfinhom, exponent = Inf, rinterval = c(0, 12), nsim = 39, leaveoneout = FALSE, use.theory = FALSE, verbose = FALSE)$p.value, 3)), adj = c(1, 1.5), col = "red")
}

#V2019
for(i in names(ModsInhomV2019)) {
  temp <- envelope(V2019_ppp, fun = Linhom, nsim = 399, nrank = 10, simulate = expression(as.ppp(arrayInd(sample(2500, prob = predict(ModsInhomV2019[[i]], ngrid = 50)$v, size = sum(!V2019_ants_50$v)), .dim = c(50, 50)), W = owin(c(0.5, 50.5), c(0.5, 50.5)))), funargs = list(r = 0:12, rmax = 12), global = FALSE, verbose = FALSE, savefuns = TRUE, savepatterns = TRUE)
  plot(temp, . - mmean ~ r, legend = FALSE, main = paste(i, "- L(inhom)"), cex.axis = 0.8)
  text(x = 12, y = -1.5, paste("p = ", round(mad.test(temp, use.theory = FALSE)$p.value, 3)), adj = c(1, 0))
  text(x = 12, y = -1.5, paste("p = ", round(bits.test(ModsInhomV2019[[i]], fun = Linhom, exponent = Inf, rinterval = c(0, 12), nsim = 39, leaveoneout = FALSE, use.theory = FALSE, verbose = FALSE)$p.value, 3)), adj = c(1, 1.5), col = "red")
  
  temp2 <- envelope(V2019_ppp, fun = pcfinhom, nsim = 399, nrank = 10, r = 0:12, rmax = 12, verbose = FALSE, simulate = temp, savefuns = TRUE)
  plot(temp2, . - mmean + 1 ~ r, legend = FALSE, main = paste(i, "- pcf(inhom)"), cex.axis = 0.8)
  text(x = 12, y = 0.6, paste("p = ", round(mad.test(temp2, use.theory = FALSE)$p.value, 3)), adj = c(1, 0))
  text(x = 12, y = 0.6, paste("p = ", round(bits.test(ModsInhomV2019[[i]], fun = pcfinhom, exponent = Inf, rinterval = c(0, 12), nsim = 39, leaveoneout = FALSE, use.theory = FALSE, verbose = FALSE)$p.value, 3)), adj = c(1, 1.5), col = "red")
}
Autocorrelación espacial residual según las distancias r entre colonias (en metros) en las parcelas J en 2001 (arriba) y V en 2019 (abajo) una vez descontada la variación de intensidad estimada por cada modelo Poisson heterogéneo (HET, BRD, BRDHET y DIST^2^, ver texto) según las funciones L de Ripley y de correlación de pares (g). Líneas: valores observados; sombreado: intervalos críticos del 95% bajo cada modelo Poisson hheterogéneo; los p-valores corresponden a pruebas globales de Máxima Desviación Absoluta (en negro, probablemente conservativas) y pruebas globales corregidas (en rojo; MonteCarlo de dos etapas).Autocorrelación espacial residual según las distancias r entre colonias (en metros) en las parcelas J en 2001 (arriba) y V en 2019 (abajo) una vez descontada la variación de intensidad estimada por cada modelo Poisson heterogéneo (HET, BRD, BRDHET y DIST^2^, ver texto) según las funciones L de Ripley y de correlación de pares (g). Líneas: valores observados; sombreado: intervalos críticos del 95% bajo cada modelo Poisson hheterogéneo; los p-valores corresponden a pruebas globales de Máxima Desviación Absoluta (en negro, probablemente conservativas) y pruebas globales corregidas (en rojo; MonteCarlo de dos etapas).

Figura 18: Autocorrelación espacial residual según las distancias r entre colonias (en metros) en las parcelas J en 2001 (arriba) y V en 2019 (abajo) una vez descontada la variación de intensidad estimada por cada modelo Poisson heterogéneo (HET, BRD, BRDHET y DIST2, ver texto) según las funciones L de Ripley y de correlación de pares (g). Líneas: valores observados; sombreado: intervalos críticos del 95% bajo cada modelo Poisson hheterogéneo; los p-valores corresponden a pruebas globales de Máxima Desviación Absoluta (en negro, probablemente conservativas) y pruebas globales corregidas (en rojo; MonteCarlo de dos etapas).

Los modelos de procesos de puntos Gibbs incorporan interacciones entre puntos de manera explícita 3 que pueden modelarse además de las covariables que generan diferencias espaciales en la intensidad del proceso. Esos patrones pueden surgir, por ejemplo, de procesos de aparición y desaparición de puntos en el tiempo, (como establecimiento y muerte de plantas o colonias de hormigas) en los que la probabilidad de esos eventos esté afectada por la distancia a los puntos ya presentes 4. Si esas interacciones son negativas se generan patrones de inhibición o repulsión, y si son positivas, de agregamiento moderado.

Existen modelos Gibbs que proponen restricciones estrictas (hardcore: no se permiten puntos a menos de cierta distancia h entre sí) o interacciones parciales dependientes de la distancia entre puntos, como el modelo de Strauss que penaliza la ocurrencia de puntos más cercanos que cierta distancia r entre sí, a partir de la cual ya no hay interacción (i.e., luego de cierta distancia unbral se comporta como un proceso Poisson). Y también existen modelos de interacción combinados (e.g., el modelo Strauss-Hardcore no permite puntos a menos de una distancia hc y los penaliza o favorece en una magnitud fija hasta una distancia r > hc). Estos últimos son particularmente apropiados para modelar los patrones de puntos de este muestreo en el que por diseño (resolución) no puede haber puntos medidos a menos de 1 m de distancia entre sí (una restricción hardcore) pero se desean explorar los patrones de agregamiento/repulsión a distancias mayores.

Los perfiles de interacción entre pares de puntos cercanos (evaluados a distancias crecientes en intervalos de 1 m) resultan similares bajo modelos que asumen una intensidad de base homogénea (CSR) o bajo los que asumen y estiman una intensidad variable relacionada de distintos modos con la distancia al borde y el tipo de vegetación (solo se muestran los del modelo DIST2, pero los otros también son similares; Fig. 19), y son compatibles con lo ya observado:
* un patrón de leves interacciones positivas (ɣ > 1) a distancias cercanas durante 2001, que generaría un leve agregamiento entre píxeles con nido, * una interacción negativa (ɣ < 1) muy marcada hasta los 6–7 m en V2019 compatible con el patrón uniforme observado, y * un perfil de interacciones negativas menos marcadas (y exceptuando al primer intervalo de 1–2 m, es decir entre píxeles adyacentes) en J2019.

#interaction plots with function PairPiece (escalonada, un gamma por valor de r); no puede haber puntos a < 1 m
par(mfrow = c(2, 2), mar = c(2, 2, 1.5, 1))
layout(matrix(c(1, 2, 0, 5, 6, 3, 4, 0, 7, 8), 2, 5, byrow = TRUE), widths = c(2, 2, 1, 2, 2))
plot(fitin(ppm(quadscheme(J2001_ppp, nd = 50) ~ 1, PairPiece(1:12), correction = "isotropic")), main = "Profile J2001:CSR", legend = FALSE)
plot(fitin(ppm(quadscheme(J2019_ppp, nd = 50) ~ 1, PairPiece(1:12), correction = "isotropic")), main = "Profile J2019:CSR", legend = FALSE)
plot(fitin(ppm(quadscheme(V2001_ppp, nd = 50) ~ 1, PairPiece(1:12), correction = "isotropic")), main = "Profile V2001:CSR", legend = FALSE)
plot(fitin(ppm(quadscheme(V2019_ppp, nd = 50) ~ 1, PairPiece(1:12), correction = "isotropic")), main = "Profile V2019:CSR", legend = FALSE)

#idem pero modelos con polinomio de vegetación
plot(fitin(ppm(quadscheme(J2001_ppp, nd = 50) ~ polynom(J2001_vegdist_50, 2), PairPiece(1:12), correction = "isotropic")), main = "Profile J2001:DIST²", legend = FALSE)
plot(fitin(ppm(quadscheme(J2019_ppp, nd = 50) ~ polynom(J2019_vegdist_50, 2), PairPiece(1:12), correction = "isotropic")), main = "Profile J2019:DIST²", legend = FALSE)
plot(fitin(ppm(quadscheme(V2001_ppp, nd = 50) ~ polynom(V2019_vegdist_50, 2), PairPiece(1:12), correction = "isotropic")), main = "Profile V2001:DIST²", legend = FALSE)
plot(fitin(ppm(quadscheme(V2019_ppp, nd = 50) ~ polynom(V2019_vegdist_50, 2), PairPiece(1:12), correction = "isotropic")), main = "Profile V2019:DIST²", legend = FALSE)
Perfiles de interacción entre pares de puntos cercanos (evaluados a distancias crecientes en intervalos de 1 m) sin contemplar los efectos de vegetación (modelos CSR, izquierda) o bien luego de descontar los efectos de la interacción entre distancia al borde y tipo de vegetación (modelados con una función cuadrática de la distancia al borde: DIST², derecha) para cada parcela y año de muestreo. Valores mayores a 1 implican agregamiento, y menores, repulsión.

Figura 19: Perfiles de interacción entre pares de puntos cercanos (evaluados a distancias crecientes en intervalos de 1 m) sin contemplar los efectos de vegetación (modelos CSR, izquierda) o bien luego de descontar los efectos de la interacción entre distancia al borde y tipo de vegetación (modelados con una función cuadrática de la distancia al borde: DIST², derecha) para cada parcela y año de muestreo. Valores mayores a 1 implican agregamiento, y menores, repulsión.

Para modelar esas interacciones negativas/positivas de la manera más simple posible (un modelo híbrido Strauss-Hard, con restricción total hasta 1 m y un único valor de penalización/favorecimiento hasta una distancia crítica r), la distancia crítica que maximiza la función de verosimilitud está en todos los casos en un rango de 6–6.5 m (Fig. 20).

# encuentra el valor del parámetro irregular r (2-12 m cada medio pixel) del Strauss+Hard que maximiza la verosimilitud [con hc = .99, los dummies restringidos por quadscheme como en los ajustes anteriores y corrección de borde isotrópica]
par(mfrow = c(2, 2), mar = c(2, 2, 1.5, 1))
plot(profilepl(data.frame(r=seq(2, 12, 0.5), hc = .99), StraussHard, quadscheme(J2001_ppp, nd = 50) ~ 1, correction = "isotropic", AIC = TRUE, fast = FALSE, verbose = FALSE), main = "Strauss+Hard J2001")
plot(profilepl(data.frame(r=seq(2, 12, 0.5), hc = .99), StraussHard, quadscheme(J2019_ppp, nd = 50) ~ 1, correction = "isotropic", AIC = TRUE, fast = FALSE, verbose = FALSE), main = "Strauss+Hard J2019")
plot(profilepl(data.frame(r=seq(2, 12, 0.5), hc = .99), StraussHard, quadscheme(V2001_ppp, nd = 50) ~ 1, correction = "isotropic", AIC = TRUE, fast = FALSE, verbose = FALSE), main = "Strauss+Hard V2001")
plot(profilepl(data.frame(r=seq(2, 12, 0.5), hc = .99), StraussHard, quadscheme(V2019_ppp, nd = 50) ~ 1, correction = "isotropic", AIC = TRUE, fast = FALSE, verbose = FALSE), main = "Strauss+Hard V2019")
Verosimilitud de modelos híbridos Gibbs (con interacción tipo Strauss-Hard) de intensidad constante en función del parámetro libre r (en m), la distancia máxima de interacción entre nidos cercanos (hc, la distancia mínima entre nidos vecinos, fue fijada en 99 cm, dado que 1 m es el tamaño de píxel o la máxima resolución del censo de colonias). La línea verde señala el valor máximo encontrado.

Figura 20: Verosimilitud de modelos híbridos Gibbs (con interacción tipo Strauss-Hard) de intensidad constante en función del parámetro libre r (en m), la distancia máxima de interacción entre nidos cercanos (hc, la distancia mínima entre nidos vecinos, fue fijada en 99 cm, dado que 1 m es el tamaño de píxel o la máxima resolución del censo de colonias). La línea verde señala el valor máximo encontrado.

Tomando esos dos valores críticos (hc < 1 m y r = 6 m) para definir una interacción potencial tipo Strauss-Hard entre píxeles con nidos cercanos, el ajuste de los modelos de intensidad heterogénea en relación con el tipo y distancia a la vegetación leñosa (ya ensayados como modelos Poisson) muestra tendencias de primer orden (i.e., la intensidad que correspondería a un patrón sin puntos) muy similares (Fig. 21).

# idem modelos Poisson pero ahora con StraussHard (r = 6, hc = .99) y corrección de borde isotrópica
par(mfrow = c(4, 4), mar = c(2, 2, 1.5, 1))

#J2001
rangox = seq(min(J2001_vegdist_50), max(J2001_vegdist_50), 0.2)
## Modelo Strauss+Hard CSR
ModJ2001_SHCSR <- ppm(quadscheme(J2001_ppp, nd = 50) ~ 1, StraussHard(6, .99), correction = "isotropic")
## Modelo Strauss+Hard HET categórico: vegetación
ModJ2001_SHVeg <- ppm(quadscheme(J2001_ppp, nd = 50) ~ J2001_veg_50, StraussHard(6, .99), correction = "isotropic")
barplot(c(exp(coef(ModJ2001_SHVeg)[1]), exp(sum(coef(ModJ2001_SHVeg)[c(1,2)]))), names.arg = c("veg", "sd"), ylab = "First order trend", main = "J2001:SH+HET", col = "white")
## Modelo Strauss+Hard HET continuo: distancia a borde (independiente del tipo de vegetación)
ModJ2001_SHAbsDist <- ppm(quadscheme(J2001_ppp, nd = 50) ~ abs(J2001_vegdist_50), StraussHard(6, .99), correction = "isotropic")
plot(effectfun(ModJ2001_SHAbsDist), main = "J2001:SH+BRD")
## Modelo Strauss+Hard HET interacción distancia absoluta / veg
ModJ2001_SHDistxVeg <- ppm(quadscheme(J2001_ppp, nd = 50) ~ abs(J2001_vegdist_50) / J2001_veg_50, StraussHard(6, .99), correction = "isotropic")
plot(x = rangox[rangox <= 0], y = exp(coef(ModJ2001_SHDistxVeg)[1] - coef(ModJ2001_SHDistxVeg)[2] * rangox[rangox <= 0]), type = "l", ylab = "First order trend", main = "J2001:SH+BRDHET", xlim = range(rangox), ylim = c(0, 0.05))
lines(x = rangox[rangox >= 0], y = exp(sum(coef(ModJ2001_SHDistxVeg)[1]) + sum(coef(ModJ2001_SHDistxVeg)[c(2, 4)]) * rangox[rangox >= 0]))
## Modelo Strauss+Hard HET cuadrática de la distancia con signo
ModJ2001_SHDist2 <- ppm(quadscheme(J2001_ppp, nd = 50) ~ polynom(J2001_vegdist_50, 2), StraussHard(6, .99), correction = "isotropic")
plot(effectfun(ModJ2001_SHDist2), main = "J2001:SH+DIST²")

# J2019
rangox = seq(min(J2019_vegdist_50), max(J2019_vegdist_50), 0.2)
## Modelo Strauss+Hard CSR
ModJ2019_SHCSR <- ppm(quadscheme(J2019_ppp, nd = 50) ~ 1, StraussHard(6, .99), correction = "isotropic")
## Modelo Strauss+Hard HET categórico: vegetación
ModJ2019_SHVeg <- ppm(quadscheme(J2019_ppp, nd = 50) ~ J2019_veg_50, StraussHard(6, .99), correction = "isotropic")
barplot(c(exp(coef(ModJ2019_SHVeg)[1]), exp(sum(coef(ModJ2019_SHVeg)[c(1,2)]))), names.arg = c("veg", "sd"), ylab = "First order trend", main = "J2019:SH+HET", col = "white")
## Modelo Strauss+Hard HET continuo: distancia a borde (independiente del tipo de vegetación)
ModJ2019_SHAbsDist <- ppm(quadscheme(J2019_ppp, nd = 50) ~ abs(J2019_vegdist_50), StraussHard(6, .99), correction = "isotropic")
plot(effectfun(ModJ2019_SHAbsDist), main = "J2019:SH+BRD")
## Modelo Strauss+Hard HET interacción distancia absoluta / veg
ModJ2019_SHDistxVeg <- ppm(quadscheme(J2019_ppp, nd = 50) ~ abs(J2019_vegdist_50) / J2019_veg_50, StraussHard(6, .99), correction = "isotropic")
plot(x = rangox[rangox <= 0], y = exp(coef(ModJ2019_SHDistxVeg)[1] - coef(ModJ2019_SHDistxVeg)[2] * rangox[rangox <= 0]), type = "l", ylab = "First order trend", main = "J2019:SH+BRDHET", xlim = range(rangox), ylim = c(0, 0.05))
lines(x = rangox[rangox >= 0], y = exp(sum(coef(ModJ2019_SHDistxVeg)[1]) + sum(coef(ModJ2019_SHDistxVeg)[c(2, 4)]) * rangox[rangox >= 0]))
## Modelo Strauss+Hard HET cuadrática de la distancia con signo
ModJ2019_SHDist2 <- ppm(quadscheme(J2019_ppp, nd = 50) ~ polynom(J2019_vegdist_50, 2), StraussHard(6, .99), correction = "isotropic")
plot(effectfun(ModJ2019_SHDist2), main = "J2019:SH+DIST²")

# V2001
rangox = seq(min(V2019_vegdist_50), max(V2019_vegdist_50), 0.2)
## Modelo Strauss+Hard CSR
ModV2001_SHCSR <- ppm(quadscheme(V2001_ppp, nd = 50) ~ 1, StraussHard(6, .99), correction = "isotropic")
## Modelo Strauss+Hard HET categórico: vegetación
ModV2001_SHVeg <- ppm(quadscheme(V2001_ppp, nd = 50) ~ V2019_veg_50, StraussHard(6, .99), correction = "isotropic")
barplot(c(exp(coef(ModV2001_SHVeg)[1]), exp(sum(coef(ModV2001_SHVeg)[c(1,2)]))), names.arg = c("veg", "sd"), ylab = "First order trend", main = "V2001:SH+HET", col = "white")
## Modelo Strauss+Hard HET continuo: distancia a borde (independiente del tipo de vegetación)
ModV2001_SHAbsDist <- ppm(quadscheme(V2001_ppp, nd = 50) ~ abs(V2019_vegdist_50), StraussHard(6, .99), correction = "isotropic")
plot(effectfun(ModV2001_SHAbsDist), main = "V2001:SH+BRD")
## Modelo Strauss+Hard HET interacción distancia absoluta / veg
ModV2001_SHDistxVeg <- ppm(quadscheme(V2001_ppp, nd = 50) ~ abs(V2019_vegdist_50) / V2019_veg_50, StraussHard(6, .99), correction = "isotropic")
plot(x = rangox[rangox <= 0], y = exp(coef(ModV2001_SHDistxVeg)[1] - coef(ModV2001_SHDistxVeg)[2] * rangox[rangox <= 0]), type = "l", ylab = "First order trend", main = "V2001:SH+BRDHET", xlim = range(rangox), ylim = c(0, 0.04))
lines(x = rangox[rangox >= 0], y = exp(sum(coef(ModV2001_SHDistxVeg)[1]) + sum(coef(ModV2001_SHDistxVeg)[c(2, 4)]) * rangox[rangox >= 0]))
## Modelo Strauss+Hard HET cuadrática de la distancia con signo
ModV2001_SHDist2 <- ppm(quadscheme(V2001_ppp, nd = 50) ~ polynom(V2019_vegdist_50, 2), StraussHard(6, .99), correction = "isotropic")
plot(effectfun(ModV2001_SHDist2), main = "V2001:SH+DIST²")

# V2019
rangox = seq(min(V2019_vegdist_50), max(V2019_vegdist_50), 0.2)
## Modelo Strauss+Hard CSR
ModV2019_SHCSR <- ppm(quadscheme(V2019_ppp, nd = 50) ~ 1, StraussHard(6, .99), correction = "isotropic")
## Modelo Strauss+Hard HET categórico: vegetación
ModV2019_SHVeg <- ppm(quadscheme(V2019_ppp, nd = 50) ~ V2019_veg_50, StraussHard(6, .99), correction = "isotropic")
barplot(c(exp(coef(ModV2019_SHVeg)[1]), exp(sum(coef(ModV2019_SHVeg)[c(1,2)]))), names.arg = c("veg", "sd"), ylab = "First order trend", main = "V2019:SH+HET", col = "white")
## Modelo Strauss+Hard HET continuo: distancia a borde (independiente del tipo de vegetación)
ModV2019_SHAbsDist <- ppm(quadscheme(V2019_ppp, nd = 50) ~ abs(V2019_vegdist_50), StraussHard(6, .99), correction = "isotropic")
plot(effectfun(ModV2019_SHAbsDist), main = "V2019:SH+BRD")
## Modelo Strauss+Hard HET interacción distancia absoluta / veg
ModV2019_SHDistxVeg <- ppm(quadscheme(V2019_ppp, nd = 50) ~ abs(V2019_vegdist_50) / V2019_veg_50, StraussHard(6, .99), correction = "isotropic")
plot(x = rangox[rangox <= 0], y = exp(coef(ModV2019_SHDistxVeg)[1] - coef(ModV2019_SHDistxVeg)[2] * rangox[rangox <= 0]), type = "l", ylab = "First order trend", main = "V2019:SH+BRDHET", xlim = range(rangox), ylim = c(0, 0.27))
lines(x = rangox[rangox >= 0], y = exp(sum(coef(ModV2019_SHDistxVeg)[1]) + sum(coef(ModV2019_SHDistxVeg)[c(2, 4)]) * rangox[rangox >= 0]))
## Modelo Strauss+Hard HET cuadrática de la distancia con signo
ModV2019_SHDist2 <- ppm(quadscheme(V2019_ppp, nd = 50) ~ polynom(V2019_vegdist_50, 2), StraussHard(6, .99), correction = "isotropic")
plot(effectfun(ModV2019_SHDist2), main = "V2019:SH+DIST²")

ModsSH_All <- list("J2001:SHCSR" = ModJ2001_SHCSR,
                   "J2001:SHHET" = ModJ2001_SHVeg,
                   "J2001:SHBRD" = ModJ2001_SHAbsDist,
                   "J2001:SHBRDHET" = ModJ2001_SHDistxVeg,
                   "J2001:SHDIST²" = ModJ2001_SHDist2,
                   "J2019:SHCSR" = ModJ2019_SHCSR,
                   "J2019:SHHET" = ModJ2019_SHVeg,
                   "J2019:SHBRD" = ModJ2019_SHAbsDist,
                   "J2019:SHBRDHET" = ModJ2019_SHDistxVeg,
                   "J2019:SHDIST²" = ModJ2019_SHDist2,
                   "V2001:SHCSR" = ModV2001_SHCSR,
                   "V2001:SHHET" = ModV2001_SHVeg,
                   "V2001:SHBRD" = ModV2001_SHAbsDist,
                   "V2001:SHBRDHET" = ModV2001_SHDistxVeg,
                   "V2001:SHDIST²" = ModV2001_SHDist2,
                   "V2019:SHCSR" = ModV2019_SHCSR,
                   "V2019:SHHET" = ModV2019_SHVeg,
                   "V2019:SHBRD" = ModV2019_SHAbsDist,
                   "V2019:SHBRDHET" = ModV2019_SHDistxVeg,
                   "V2019:SHDIST²" = ModV2019_SHDist2)
Idem Fig. 8 pero bajo modelos Gibbs híbridos de tipo Strauss-Hard, con parámetros hc = 0.99 m (distancia mínima entre píxeles con nidos) y r = 6 m (distancia máxima de inetarcción), para cada una de las parcelas y años.

Figura 21: Idem Fig. 8 pero bajo modelos Gibbs híbridos de tipo Strauss-Hard, con parámetros hc = 0.99 m (distancia mínima entre píxeles con nidos) y r = 6 m (distancia máxima de inetarcción), para cada una de las parcelas y años.

Las comparaciones de estos modelos entre sí y con los modelos Poisson previos ajustados para cada una de las parcelas y años (recordando que en el caso de V2001 la vegetación usada al ajustar los modelos heterogéneos es la medida en el mismo lugar en 2019) son:

Modelos Poisson y Strauss+Hard ajustados: J2001

Modelo Objeto Dependencia Interacción Intensidad/Tendencia AIC
CSR ModJ2001_CSR — — 0.03 684.98
ModJ2001_SHCSR 1.16 0.017 667.88
HET ModJ2001_Veg veg/sd — λveg: 0.017 / λsd: 0.051 664.9
ModJ2001_SHVeg 1.147 βveg: 0.01 / βsd: 0.029 649.26
BRD ModJ2001_AbsDist abs(dist) — 0.079 × 0.462m 678.61
ModJ2001_SHAbsDist 1.155 0.044 × 0.467m 662.28
BRDHET ModJ2001_DistxVeg abs(dist) * veg/sd — λveg: 0.059 × 0.37m / λsd: 0.059 × 0.896m 659.74
ModJ2001_SHDistxVeg 1.144 βveg: 0.038 × 0.352m / βsd: 0.038 × 0.83m 644.54
DIST2 ModJ2001_Dist2 dist2 — 0.04 × 1.553m × 0.845m2 659.79
ModJ2001_SHDist2 1.144 0.024 × 1.535m × 0.828m2 644.6

Modelos Poisson y Strauss+Hard ajustados: J2019

Modelo Objeto Dependencia Interacción Intensidad/Tendencia AIC
CSR ModJ2019_CSR — — 0.013 344.93
ModJ2019_SHCSR 0.767 0.013 347.77
HET ModJ2019_Veg veg/sd — λveg: 0.009 / λsd: 0.022 340.35
ModJ2019_SHVeg 0.715 βveg: 0.013 / βsd: 0.036 342.34
BRD ModJ2019_AbsDist abs(dist) — 0.042 × 0.423m 339.3
ModJ2019_SHAbsDist 0.711 0.072 × 0.392m 340.77
BRDHET ModJ2019_DistxVeg abs(dist) * veg/sd — λveg: 0.028 × 0.427m / λsd: 0.028 × 0.847m 335.65
ModJ2019_SHDistxVeg 0.679 βveg: 0.047 × 0.414m / βsd: 0.047 × 0.909m 336.65
DIST2 ModJ2019_Dist2 dist2 — 0.019 × 1.404m × 0.868m2 335.66
ModJ2019_SHDist2 0.678 0.033 × 1.476m × 0.872m2 336.58

Modelos Poisson y Strauss+Hard ajustados: V2001

Modelo Objeto Dependencia Interacción Intensidad/Tendencia AIC
CSR ModV2001_CSR — — 0.036 773.7
ModV2001_SHCSR 1.134 0.021 759.82
HET ModV2001_Veg veg/sd — λveg: 0.031 / λsd: 0.042 773.51
ModV2001_SHVeg 1.13 βveg: 0.019 / βsd: 0.024 760.42
BRD ModV2001_AbsDist abs(dist) — 0.045 × 0.842m 774.5
ModV2001_SHAbsDist 1.135 0.027 × 0.825m 760.1
BRDHET ModV2001_DistxVeg abs(dist) * veg/sd — λveg: 0.042 × 0.826m / λsd: 0.042 × 0.976m 775.27
ModV2001_SHDistxVeg 1.133 βveg: 0.026 × 0.814m / βsd: 0.026 × 0.912m 761.56
DIST2 ModV2001_Dist2 dist2 — 0.037 × 1.093m × 0.989m2 775.57
ModV2001_SHDist2 1.131 0.022 × 1.065m × 0.981m2 762.04

Modelos Poisson y Strauss+Hard ajustados: V2019

Modelo Objeto Dependencia Interacción Intensidad/Tendencia AIC
CSR ModV2019_CSR — — 0.017 429.26
ModV2019_SHCSR 0.34 0.018 391.31
HET ModV2019_Veg veg/sd — λveg: 0.015 / λsd: 0.019 430.64
ModV2019_SHVeg 0.328 βveg: 0.056 / βsd: 0.096 389.02
BRD ModV2019_AbsDist abs(dist) — 0.041 × 0.511m 425.95
ModV2019_SHAbsDist 0.297 0.277 × 0.394m 380.41
BRDHET ModV2019_DistxVeg abs(dist) * veg/sd — λveg: 0.04 × 0.503m / λsd: 0.04 × 0.547m 427.84
ModV2019_SHDistxVeg 0.295 βveg: 0.255 × 0.387m / βsd: 0.255 × 0.482m 380.03
DIST2 ModV2019_Dist2 dist2 — 0.023 × 1.045m × 0.859m2 428.71
ModV2019_SHDist2 0.296 0.129 × 1.118m × 0.815m2 381.21

En todos los casos los modelos Gibbs con interacción entre puntos cercanos (hasta 6 m entre sí) muestran un ajuste igual o mejor que sus contrapartes que suponen que no existe interacción (Poisson). En tres de los casos el ajuste es sustancialmente mejor, y en el otro (J2019, en el que la autocorrelación espacial es compatible con un modelo aleatorio) genera modelos equivalentes (i.e., no hay evidencia suficiente para distinguir entre los modelos con y sin interacción, al menos usando un solo valor de interacción para puntos entre 1 y 6 m de distancia entre sí 5). Ademas, en general no hay distinciones notorias entre los dos modos de modelar la relación heterogénea entre tipos de vegetación con la distancia al borde de la vegetación (BRDHET o DIST2), y las estimaciones de la variación de la intensidad con la distancia son muy similares independientemente de cómo resulte la interacción entre puntos cercanos.

En los dos casos de 2001 (densidades altas), los valores de interacción estimada entre puntos son levemente mayores a 1 (asociado a un agregamiento de nidos a distancias < 6 m), mientras que en 2019 son menores a 1 (notoriamente en V2019, que resulta en un patrón espacial uniforme).

Además, y en cierto modo crucial, estos modelos con interacción parecen acomodar mejor los patrones de autocorrelación observados que los modelos Poisson que suponían no interacción: las funciones de autocorrelación K de los residuos ya no exceden los límites de lo esperado (Fig. 22).

par(mfrow = c(4, 5), mar = c(2, 2, 1.5, 1))
#gráficos de K de Ripley de los residuos de Pearson
for(i in names(ModsSH_All)) {
  plot(Kres(ModsSH_All[[i]], r = seq(0, 12.5, 0.5), correction = "iso"), ires ~ r, shade = c("ihi", "ilo"), legend = FALSE, main = i)
  abline(h = 0, lty = 2, col = "red", lwd = 0.5)
}
Funciones de autocorrelación K de Ripley de los residuos para cada parcela y año de todos los modelos Gibbs ensayados. Cuando el modelo ajustado es exacto y no hay autocorrelación residual, el valor esperado de la función es aproximadamente cero.

Figura 22: Funciones de autocorrelación K de Ripley de los residuos para cada parcela y año de todos los modelos Gibbs ensayados. Cuando el modelo ajustado es exacto y no hay autocorrelación residual, el valor esperado de la función es aproximadamente cero.

Síntesis: modelos Gibbs con parcelas como réplicas

Como es esperable a partir de estos resultados, el modelo con interacción entre puntos también resulta en un ajuste notoriamente superior cuando se trata a las dos parcelas como réplicas de un mismo proceso espacial (variable en el tiempo): la interacción negativa entre nidos a distancias menores de 6 m durante 2019 y la positiva durante 2001 mejoran el ajuste de los modelos espacialmente heterogéneos respecto a sus equivalentes Poisson (= sin interacción).

#Ajuste a 2001 (como réplicas) de los modelos StraussHard(6, .99)
ModX2001_SHCSR <- mppm(Pix ~ 1, X2001_HF, StraussHard(6, .99), correction = "isotropic")

ModX2001_SHVeg <- mppm(Pix ~ Veg, X2001_HF, StraussHard(6, .99), correction = "isotropic")

ModX2001_SHAbsDist <- mppm(Pix ~ abs(Dist), X2001_HF, StraussHard(6, .99), correction = "isotropic")

ModX2001_SHDistxVeg <- mppm(Pix ~ abs(Dist) / Veg, X2001_HF, StraussHard(6, .99), correction = "isotropic")

ModX2001_SHDist2 <- mppm(Pix ~ Dist + I(Dist^2), X2001_HF, StraussHard(6, .99), correction = "isotropic")

#Ajuste a 2019 (como réplicas) de los modelos StraussHard(6, .99)
ModX2019_SHCSR <- mppm(Pix ~ 1, X2019_HF, StraussHard(6, .99), correction = "isotropic")

ModX2019_SHVeg <- mppm(Pix ~ Veg, X2019_HF, StraussHard(6, .99), correction = "isotropic")

ModX2019_SHAbsDist <- mppm(Pix ~ abs(Dist), X2019_HF, StraussHard(6, .99), correction = "isotropic")

ModX2019_SHDistxVeg <- mppm(Pix ~ abs(Dist) / Veg, X2019_HF, StraussHard(6, .99), correction = "isotropic")

ModX2019_SHDist2 <- mppm(Pix ~ Dist + I(Dist^2), X2019_HF, StraussHard(6, .99), correction = "isotropic")

Modelos Poisson y Strauss+Hard ajustados: réplicas 2001

Modelo Objeto Dependencia Interacción Intensidad/Tendencia AIC
CSR ModX2001_CSR — — 0.033 1457.71
ModX2001_SHCSR 1.147 0.019 1422.74
HET ModX2001_Veg veg/sd — λveg: 0.024 / λsd: 0.047 1441.26
ModX2001_SHVeg 1.139 βveg: 0.014 / βsd: 0.026 1409.12
BRD ModX2001_AbsDist abs(dist) — 0.052 × 0.714m 1453.77
ModX2001_SHAbsDist 1.148 0.029 × 0.707m 1418.55
BRDHET ModX2001_DistxVeg abs(dist) * veg/sd — λveg: 0.044 × 0.659m / λsd: 0.044 × 1.028m 1441.24
ModX2001_SHDistxVeg 1.14 βveg: 0.027 × 0.647m / βsd: 0.027 × 0.961m 1409.1
DIST2 ModX2001_Dist2 dist2 — 0.036 × 1.252m × 0.972m2 1442.34
ModX2001_SHDist2 1.139 0.021 × 1.22m × 0.963m2 1410.62

Modelos Poisson y Strauss+Hard ajustados: réplicas 2019

Modelo Objeto Dependencia Interacción Intensidad/Tendencia AIC
CSR ModX2019_CSR — — 0.015 773.54
ModX2019_SHCSR 0.53 0.033 742.24
HET ModX2019_Veg veg/sd — λveg: 0.012 / λsd: 0.02 769.87
ModX2019_SHVeg 0.497 βveg: 0.026 / βsd: 0.057 732.9
BRD ModX2019_AbsDist abs(dist) — 0.042 × 0.461m 762.18
ModX2019_SHAbsDist 0.477 0.137 × 0.388m 721.58
BRDHET ModX2019_DistxVeg abs(dist) * veg/sd — λveg: 0.036 × 0.447m / λsd: 0.036 × 0.635m 760.63
ModX2019_SHDistxVeg 0.464 βveg: 0.108 × 0.389m / βsd: 0.108 × 0.641m 716.89
DIST2 ModX2019_Dist2 dist2 — 0.022 × 1.193m × 0.851m2 761.55
ModX2019_SHDist2 0.464 0.063 × 1.285m × 0.836m2 717.74

Modelo global espacio-temporal (diseño completo)

Habiendo revisado ya las opciones de ajustes de modelos bajo distintos supuestos y niveles de agregación (muestras independientes o réplicas espaciales), encontramos evidencia de que:

  • los nidos no se disponen aleatoriamente: hay una relación positiva con la menor distancia a la interfase entre parches (vegetación leñosa–suelo desnudo), y
  • aún contemplando esa relación, existe una interacción espacial entre nidos cercanos (a menos de 6–7 m entre sí) que varía de signo (= atractiva o repulsiva) y magnitud en el tiempo (o con la densidad, que cambió notoriamente entre los periodos muestreados).

En consecuencia, se pueden ajustar modelos al diseño general de muestreo (pseudo-experimental), considerando a los patrones de nidos como provenientes de un par de réplicas de un mismo proceso espacial que es variable en el tiempo [recordando la limitación dada por que uno de los mapas de vegetación no está disponible y se infiere a partir de la relevada en el mismo lugar en el otro periodo]. Por ejemplo, a partir de las inferencias anteriores, se pueden ajustar y comparar modelos Gibbs con la tendencia del proceso espacial (= la intensidad de un patrón sin nidos) afectada por una relación común (= única para todas las parcelas y tiempos) con la distancia al borde y el tipo de vegetación, y un factor binario “Año” que solo afecta a la ordenada (= la intensidad de un proceso sin puntos cuando Dist = 0) y al parámetro ɣ de interacción Strauss entre nidos a menos de 6 m de distancia entre sí:

#Ajuste a todas las muestras de modelos Poisson con Year como factor (afectando la tendencia del proceso) para comparar con los Strauss+Hard
ModFull_CSR <- mppm(Pix ~ Year, Full_HF, correction = "isotropic")

ModFull_Veg <- mppm(Pix ~ Year + Veg, Full_HF, correction = "isotropic")

ModFull_AbsDist <- mppm(Pix ~ Year + abs(Dist), Full_HF, correction = "isotropic")

ModFull_DistxVeg <- mppm(Pix ~ Year + abs(Dist) / Veg, Full_HF, correction = "isotropic")

ModFull_Dist2 <- mppm(Pix ~ Year + Dist + I(Dist^2), Full_HF, correction = "isotropic")

#Ajuste a todas las muestras como modelos StraussHard(6, .99) con Year como factor (afectando la tendencia del proceso y la interacción entre nidos)
ModFull_SHCSR <- mppm(Pix ~ Year, Full_HF, StraussHard(6, .99), iformula = ~ Interaction:Year, correction = "isotropic")

ModFull_SHVeg <- mppm(Pix ~ Year + Veg, Full_HF, StraussHard(6, .99), iformula = ~ Interaction:Year, correction = "isotropic")

ModFull_SHAbsDist <- mppm(Pix ~ Year + abs(Dist), Full_HF, StraussHard(6, .99), iformula = ~ Interaction:Year, correction = "isotropic")

ModFull_SHDistxVeg <- mppm(Pix ~ Year + abs(Dist) / Veg, Full_HF, StraussHard(6, .99), iformula = ~ Interaction:Year, correction = "isotropic")

ModFull_SHDist2 <- mppm(Pix ~ Year + Dist + I(Dist^2), Full_HF, StraussHard(6, .99), iformula = ~ Interaction:Year, correction = "isotropic")

Modelos Poisson y Strauss+Hard: un solo diseño (con dos réplicas espaciales por año)

Modelo Objeto Dependencia Interacción Intensidad/Tendencia AIC
CSR ModFull_CSR año — β2001: 0.033
β2019: 0.015		 2231.25
HOM+ ModFull_SHCSR año ɣ2001: 1.147
ɣ2019: 0.53		 β2001: 0.019
β2019: 0.033		 2166.97
HET ModFull_Veg año + veg/sd — 2001: βveg: 0.024 / βsd: 0.046
2019: βveg: 0.011 / βsd: 0.021		 2209.3
HET+ ModFull_SHVeg año + veg/sd ɣ2001: 1.138
ɣ2019: 0.503		 2001: βveg: 0.014 / βsd: 0.027
2019: βveg: 0.027 / βsd: 0.053		 2142.3
BRD ModFull_AbsDist año + abs(dist) — 2001: 0.061 × 0.622m
2019: 0.029 × 0.622m		 2216.34
BRD+ ModFull_SHAbsDist año + abs(dist) ɣ2001: 1.149
ɣ2019: 0.49		 2001: 0.038 × 0.571m
2019: 0.082 × 0.571m		 2144.63
BRDHET ModFull_DistxVeg año + abs(dist) * veg/sd — 2001: βveg: 0.052 × 0.583m / βsd: 0.052 × 0.891m
2019: βveg: 0.025 × 0.583m / βsd: 0.025 × 0.891m		 2200.03
BRDHET+ ModFull_SHDistxVeg año + abs(dist) * veg/sd ɣ2001: 1.14
ɣ2019: 0.473		 2001: βveg: 0.033 × 0.537m / βsd: 0.033 × 0.839m
2019: βveg: 0.073 × 0.537m / βsd: 0.073 × 0.839m		 2126.89
DIST2 ModFull_Dist2 año + dist2 — 2001: 0.039 × 1.238m × 0.936m2
2019: 0.018 × 1.238m × 0.936m2		 2202.6
DIST2+ ModFull_SHDist2 año + dist2 ɣ2001: 1.139
ɣ2019: 0.474		 2001: 0.023 × 1.251m × 0.918m2
2019: 0.052 × 1.251m × 0.918m2		 2130.04

(misma tabla, traducida y formateada para el manuscrito)
Table 1:

Model Structure Interactions Intensities/Trends AIC
CSR year — β2001: 0.033
β2019: 0.015		 2231.25
HOM+ year ɣ2001: 1.147
ɣ2019: 0.53		 β2001: 0.019
β2019: 0.033		 2166.97
HET year + veg/open — 2001: βveg: 0.024 / βopen: 0.046
2019: βveg: 0.011 / βopen: 0.021		 2209.3
HET+ year + veg/open ɣ2001: 1.138
ɣ2019: 0.503		 2001: βveg: 0.014 / βopen: 0.027
2019: βveg: 0.027 / βopen: 0.053		 2142.3
BRD year + abs(dist) — 2001: 0.061 × 0.622m
2019: 0.029 × 0.622m		 2216.34
BRD+ year + abs(dist) ɣ2001: 1.149
ɣ2019: 0.49		 2001: 0.038 × 0.571m
2019: 0.082 × 0.571m		 2144.63
BRDHET year + abs(dist) * veg/open — 2001: βveg: 0.052 × 0.583m / βopen: 0.052 × 0.891m
2019: βveg: 0.025 × 0.583m / βopen: 0.025 × 0.891m		 2200.03
BRDHET+ year + abs(dist) * veg/open ɣ2001: 1.14
ɣ2019: 0.473		 2001: βveg: 0.033 × 0.537m / βopen: 0.033 × 0.839m
2019: βveg: 0.073 × 0.537m / βopen: 0.073 × 0.839m		 2126.89
DIST2 year + dist2 — 2001: 0.039 × 1.238m × 0.936m2
2019: 0.018 × 1.238m × 0.936m2		 2202.6
DIST2+ year + dist2 ɣ2001: 1.139
ɣ2019: 0.474		 2001: 0.023 × 1.251m × 0.918m2
2019: 0.052 × 1.251m × 0.918m2		 2130.04

Nuevamente, todos los modelos con interacción entre colonias son más adecuados que su versión con la misma estructura pero asumiendo que no hay interacción (= Poisson). El modelo que plantea una relación con la distancia al borde de la vegetación de distinta magnitud dependiendo del tipo de parche (BRDHET) es el que mejor se ajusta a los datos entre los que plantean diferentes relaciones con la vegetación leñosa.

Un modelo que plantee una única magnitud del proceso en ambos periodos y que las diferencias en la cantidad de nidos observados solo se deba a un cambio temporal en la intensidad de la interacción espacial entre nidos cercanos muestra un peor ajuste (abs(dist) * veg/sd + StraussHardaño, AIC: 2137.8), al igual que uno que plantee una intensidad variable con la distancia y el tiempo pero con una interacción entre nidos constante (año + abs(dist) * veg/sd + Strauss+Hard, AIC: 2183.7). Por cierto, tampoco hay evidencias importantes que justifiquen un modelado más complejo, como que la intensidad variable con la distancia al límite entre parches cambie entre años (año * abs(dist) * veg/sd + StraussHardaño, AIC: 2127.78; año * dist² + StraussHardaño, AIC: 2130.11).

Bajo el modelo elegido, los parámetros y relaciones estimadas son las que se fueron deduciendo de los modelados parciales previos: una probabilidad de ocurrencia de píxeles con nidos que desciende con la distancia al borde de la vegetación leñosa (con un descenso ~56% más rápido dentro del parche vegetado) y una interacción entre nidos cercanos (<6 m) que pasa de ligeramente positiva en 2001 (ɣ>1, causando un leve agregamiento de los nidos) a claramente negativa en 2019 (ɣ<<1, resultando en un marcado patrón inhibitorio en una de las dos parcelas) cuando la densidad de nidos se redujo aproximadamente a la mitad.

par(mar = c(4, 4, 4, 2))
layout(matrix(c(1, 2), 1, 2))
rangox = c(seq(min(V2019_vegdist_50), -.99, .1), seq(1, max(V2019_vegdist_50), 0.1))
#2001
plot(x = rangox[rangox <= 0], y = exp(coef(ModFull_SHDistxVeg)[1] - coef(ModFull_SHDistxVeg)[3] * rangox[rangox <= 0]), type = "l", xlab = "Distance to inter-patch border (m)", ylab = "First order term (“trend”, β)", xlim = range(rangox), ylim = c(0, 0.07))
lines(x = rangox[rangox >= 0], y = exp(sum(coef(ModFull_SHDistxVeg)[1]) + sum(coef(ModFull_SHDistxVeg)[c(3,4)]) * rangox[rangox >= 0]))
#2019
lines(x = rangox[rangox <= 0], y = exp(sum(coef(ModFull_SHDistxVeg)[c(1,2)]) - coef(ModFull_SHDistxVeg)[3] * rangox[rangox <= 0]), lty = 2, lwd = 1.5)
lines(x = rangox[rangox >= 0], y = exp(sum(coef(ModFull_SHDistxVeg)[c(1,2)]) + sum(coef(ModFull_SHDistxVeg)[c(3,4)]) * rangox[rangox >= 0]), lty = 2, lwd = 1.5)
#interactions
plot(x = c(0, 0.99, 1:6, 6.01, 7.5), y = c(0, 0, rep(exp(coef(ModFull_SHDistxVeg)[5]), 6), 1, 1), type = "l", xlab = "Distance between colonies (m)", ylab = "Second-order term (interaction, ɣ)", ylim = c(0, 1.5))
lines(x = c(0, 0.99, 1:6, 6.01, 7.5), y = c(0, 0, rep(exp(coef(ModFull_SHDistxVeg)[6]), 6), 1, 1), lwd = 1.5, lty = 2)
abline(h = 1, col = "red", lty = 3)
mtext("BRDHET+", line = -1.5, outer = TRUE, cex = 1.5)
mtext(text = expression(paste("( ~ Year + abs(dist) * veg/open + ", "Strauss-Hard"[Year], " )")), line = -3, outer = TRUE, cex = 1)
Relación de la intensidad condicional con la distancia al borde entre tipos de parches (distancias negativas hacia el interior del parche de vegetación leñosa alta; izquierda) y de las interacciones de tipo Strauss-Hard con la distancia entre colonias (derecha: ɣ > 1: interacciones positivas o de atracción, ɣ < 1: negativas o de repulsión) bajo el modelo espacial de patrones de puntos seleccionado (BRDHET+, ver texto)para explicar la distribución espacial de colonias de *Pheidole bergi* en el algarrobal protegido de la Reserva de la Biosfera de Ñacuñan. Los valores estimados para ambos componentes del modelo dependen del año (2001: líneas continuas; 2019: lineas interrumpidas).

Figura 23: Relación de la intensidad condicional con la distancia al borde entre tipos de parches (distancias negativas hacia el interior del parche de vegetación leñosa alta; izquierda) y de las interacciones de tipo Strauss-Hard con la distancia entre colonias (derecha: ɣ > 1: interacciones positivas o de atracción, ɣ < 1: negativas o de repulsión) bajo el modelo espacial de patrones de puntos seleccionado (BRDHET+, ver texto)para explicar la distribución espacial de colonias de Pheidole bergi en el algarrobal protegido de la Reserva de la Biosfera de Ñacuñan. Los valores estimados para ambos componentes del modelo dependen del año (2001: líneas continuas; 2019: lineas interrumpidas).

Al igual que los modelos ajustados de manera individual, el modelo completo que incluye interacciones entre nidos cercanos de magnitud variable dependiendo del tiempo (o densidad) parece acomodar más adecuadamente los patrones de autocorrelación observados que los modelos sin interacciones (las funciones de autocorrelación K de los residuos ya no exceden los límites de lo esperado; Fig. 24).

par(mfrow = c(2, 2), mar = c(2, 2, 1.5, 1))
#mismo gŕafico con one-liner elegante, pero problema para títulos
#with(hyperframe(Model = subfits(ModFull_SHDistxVeg)), plot(Kres(Model, r = seq(0, 12.5, 0.5), correction = "iso"), ires ~ r, shade = c("ihi", "ilo"), legend = FALSE, main = "FULL:BRDHET"))
for (i in 1:4) {
  plot(Kres(subfits(ModFull_SHDistxVeg)[[i]], r = seq(0, 12.5, 0.5), correction = "iso"), ires ~ r, shade = c("ihi", "ilo"), legend = FALSE, main = c("FULL:J2001", "FULL:J2019", "FULL:V2001", "FULL:V2019")[i])
}
Idem Fig. \@ref(fig:KresidualStraussHard) pero estimando la autocorrelación residual para cada parcela y año de acuerdo al modelo completo Strauss-Hard (BRDHET) con los parámetros estimados que se presentan en la Fig. \@ref(fig:FullmodelFig).

Figura 24: Idem Fig. 22 pero estimando la autocorrelación residual para cada parcela y año de acuerdo al modelo completo Strauss-Hard (BRDHET) con los parámetros estimados que se presentan en la Fig. 23.

Una observación relevante es que practicamente no existen en estas parcelas distancias a los bordes entre parches de vegetación leñosa y suelo desnudo que excedan las distancias en que se observan interacciones (leves y positivas en 2001, más importantes y negativas en 2019) entre nidos cercanos (= 6 m), por los que los factores detectados (cambios de intensidad con distancia a vegetación e interacciones entre nidos) están parcialmente confundidos.

Conclusiones

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

  1. ver sección 3.4.5 en Baddeley et al sobre los problemas de baja precisión (en este caso 1 m) para los análisis de patrones de puntos por estar basados en distancias: “The effects of rounding can substantially change the results of some statistical techniques, particularly those which deal with distances between neighbouring points.” (p.61).↩︎

  2. sobre las magnitudes de las distancias estimadas, ver help(“distancemap”): “If X is a binary pixel mask, the distance values computed are not the usual Euclidean distances. Instead the distance between two pixels is measured by the length of the shortest path connecting the two pixels. A path is a series of steps between neighbouring pixels (each pixel has 8 neighbours). This is the standard ‘distance transform’ algorithm of image processing (Rosenfeld and Kak, 1968; Borgefors, 1986).”↩︎

  3. Ver Cap. 13 en Baddeley et al. 2016↩︎

  4. Ver sección 13.3.5 en Baddeley et al. 2016↩︎

  5. El valor de AIC puede mejorar algo más si se utiliza una distancia de corte para la penalización de 7 m (AIC = 335.08, aunque no lo suficiente para que sea distinguible de un modelo Poisson. El ajuste sí resulta notoriamente mejor (AIC = 322.89) si se plantea un modelo multiescala que, además de la restricción de >1 m impuesta por el diseño de muestreo, favorezca los nidos inmediatamente vecinos (1–1.42 m de distancia) y luego penalice al resto hasta los 7 m. Esto, sin embargo, aumenta el riesgo de estar “sobreajustando”, generando un modelo demasiado específico a la muestra particular (i.e., una vez observado el patrón de interacción óptimo)↩︎