¿Determinar si los árboles dentro de las brechas forestales se agrupan usando R?

14

El conjunto de datos adjunto muestra aproximadamente 6000 árboles jóvenes en aproximadamente 50 huecos forestales de tamaño variable. Estoy interesado en aprender cómo crecen estos retoños dentro de sus respectivos huecos (es decir, agrupados, aleatorios, dispersos). Como saben, un enfoque tradicional sería ejecutar Global Moran's I. Sin embargo, las agregaciones de árboles dentro de agregaciones de huecos parecen ser un uso inapropiado de Moran's I. Realicé algunas estadísticas de prueba con Moran's I usando una distancia de umbral de 50 metros, que produjo resultados sin sentido (es decir, valor p = 0.0000000 ...). Es probable que la interacción entre las agregaciones de brechas produzca estos resultados. He considerado crear una secuencia de comandos para recorrer las brechas de dosel individuales y determinar la agrupación dentro de cada brecha, aunque mostrar estos resultados al público sería problemático.

¿Cuál es el mejor enfoque para cuantificar la agrupación dentro de los grupos?

ingrese la descripción de la imagen aquí

Aaron
fuente
1
Aaron, dices que intentaste ejecutar el I de Moran, ¿estás interesado en medir cómo se compara el atributo de un retoño con los atributos de los retoños vecinos (es decir, estás tratando con un patrón de puntos marcado )? El título parece implicar que solo está interesado en la ubicación de los retoños entre sí y no en sus atributos.
MannyG
@MannyG Sí, solo estoy interesado en determinar si los árboles jóvenes están agrupados en relación con las ubicaciones de otros árboles jóvenes dentro de cualquier brecha forestal dada. Solo hay una especie de interés y el tamaño de los retoños no es de interés.
Aaron

Respuestas:

7

No tiene un campo aleatorio uniforme, por lo que intentar analizar todos sus datos a la vez violará los supuestos de cualquier estadística que elija lanzar al problema. No está claro en su publicación si sus datos son un proceso de puntos marcados (es decir, diámetro o altura asociados con la ubicación de cada árbol). Si estos datos no representan un proceso de puntos marcados, no tengo idea de cómo aplicó un Moran's-I. Si los datos solo representan ubicaciones espaciales, recomendaría usar un Ripley's-K con la transformación Besag-L para estandarizar la expectativa nula en cero. Esto permite una evaluación multiescala del agrupamiento. Si sus datos tienen un valor asociado, entonces su mejor opción es un Moran's-I (LISA) local. Realmente lo miraría con ambas estadísticas. Independientemente de su elección, aún deberá recorrer cada sitio individualmente para producir resultados válidos. Aquí hay un ejemplo de código R para una simulación de Monte Carlo de Ripley's-K / Besag's-L utilizando el conjunto de datos de árbol de redwood incorporado. Debería ser bastante sencillo modificar esto para recorrer sus sitios y generar un gráfico para cada uno.

# ADD REQUIRED PACKAGES
require(sp)
require(spatstat)
options(scipen=5)

# USE REDWOOD SAPLING DATASET
spp <- SpatialPoints(coords(redwood))

###################################################
###### START BESAG'S-L MONTE CARLO  ANALYSUS ######
###################################################
# CREATE CONVEX HULL FOR ANALYSIS WINDOW                       
W=ripras(coordinates(spp)) 

# COERCE TO spatstat ppp OBJECT
spp.ppp=as.ppp(coordinates(spp), W)                     
  plot(spp.ppp) 

# ESTIMATE BANDWIDTH
area <- area.owin(W)
lambda <- spp.ppp$n/area
 ripley <- min(diff(W$xrange), diff(W$yrange))/4
   rlarge <- sqrt(1000/(pi * lambda))
     rmax <- min(rlarge, ripley)
bw <- seq(0, rmax, by=rmax/10)  

# CALCULATE PERMUTED CROSS-K AND PLOT RESULTS       
Lenv <- envelope(spp.ppp, fun="Kest", r=bw, i="1", j="2", nsim=99, nrank=5, 
                 transform=expression(sqrt(./pi)-bw), global=TRUE)            
plot(Lenv, main="Besag's-L", xlab="Distance", ylab="L(r)", legend=F, col=c("white","black","grey","grey"), 
    lty=c(1,2,2,2), lwd=c(2,1,1,1) )
     polygon( c(Lenv$r, rev(Lenv$r)), c(Lenv$lo, rev(Lenv$hi)), col="lightgrey", border="grey")
       lines(supsmu(bw, Lenv$obs), lwd=2)
       lines(bw, Lenv$theo, lwd=1, lty=2)
         legend("topleft", c(expression(hat(L)(r)), "Simulation Envelope", "theo"), pch=c(-32,22),
                col=c("black","grey"), lty=c(1,0,2), lwd=c(2,0,2), pt.bg=c("white","grey"))
Jeffrey Evans
fuente
1
¡Pero no puede simplemente usar el casco convexo como la ventana para su patrón de puntos! Recuerde, la ventana es el área en la que opera el patrón que produce los puntos. Usted sabe a priori que los árboles solo crecen en estas regiones establecidas, y debe configurar su ventana para reflejar eso. Puede mitigar esto estableciendo el rango de K (r) en algo muy pequeño, del orden de 0.3x del tamaño de sus claros, pero obtendrá resultados sesgados debido a la falta de correcciones del efecto de borde. Jeffrey está utilizando el tamaño de toda el área de estudio para definir su rmax.
Spacedman
1
En mi ejemplo, Sí, estoy usando toda la región. Por esto es exactamente por lo que recomendé recorrer cada sitio de muestra (brecha). Cada vez que se subdivide en un área de muestra específica, volvería a ejecutar el análisis. No puede tratar toda el área de estudio como su campo aleatorio porque no tiene muestreo continuo. Teniendo solo huecos muestreados, en efecto, tiene parcelas independientes. La función Kest que estoy llamando, por defecto, utiliza una corrección de borde "borde". Hay otras opciones de corrección de bordes disponibles. Yo diría que su unidad experimental es la brecha del dosel y debe analizarse como tal.
Jeffrey Evans el
1
Al pensar en esto un poco más. Realmente deberías usar polígonos que representen cada espacio como tu ventana. Si subconjusta su problema para reflejar la unidad experimental, entonces la CSR y K estarán sesgadas porque el área no refleja el tamaño de la brecha del dosel real. Este es un problema tanto en mi como en las recomendaciones de @ Spacedman.
Jeffrey Evans el
2
Tenga en cuenta que mi ejemplo extendido solo usaba una grilla gruesa porque era una forma bastante simple de crear algo con la estructura más o menos correcta. Su máscara debe verse como un mapa de sus áreas de bosque abierto. ¡Es técnicamente incorrecto intentar definir la máscara a partir de los datos!
Spacedman
1
@Spacedman. Me gusta su enfoque y ciertamente es eficiente. Mi preocupación específica es que las brechas del dosel son la unidad experimental. En su enfoque, si dos espacios son proximales, el ancho de banda podría incluir observaciones de diferentes unidades de muestreo. Además, la estadística resultante no debe reflejar el "conjunto" de unidades experimentales, sino que debe ser representativa de cada unidad e inferencia sobre el proceso espacial extraído de patrones comunes a través de las unidades experimentales. Si se trata globalmente, representa un proceso de intensidad no estacionario que viola los supuestos estadísticos.
Jeffrey Evans el
4

Lo que tiene es un patrón de puntos con una ventana que es un número de pequeñas regiones poligonales desconectadas.

Debería poder utilizar cualquiera de las pruebas package:spatstatpara CSR siempre que la alimente con una ventana correcta. Esto puede ser una serie de conjuntos de pares (x, y) que definen cada compensación o una matriz binaria de valores (0,1) sobre el espacio.

Primero definamos algo que se parezca un poco a sus datos:

set.seed(310366)
nclust <- function(x0, y0, radius, n) {
               return(runifdisc(n, radius, centre=c(x0, y0)))
             }
c = rPoissonCluster(15, 0.04, nclust, radius=0.02, n=5)
plot(c)

y pretendamos que nuestros claros son celdas cuadradas que resultan ser esto:

m = matrix(0,20,20)
m[1+20*cbind(c$x,c$y)]=1
imask = owin(c(0,1),c(0,1),mask = t(m)==1 )
pp1 = ppp(x=c$x,y=c$y,window=imask)
plot(pp1)

Entonces podemos trazar la función K de esos puntos en esa ventana. Esperamos que esto no sea CSR porque los puntos parecen agrupados dentro de las celdas. Tenga en cuenta que tengo que cambiar el rango de distancias para que sea pequeño, del orden del tamaño de la celda, de lo contrario, la función K se evalúa a distancias del tamaño de todo el patrón.

plot(Kest(pp1,r=seq(0,.02,len=20)))

Si generamos algunos puntos CSR en las mismas celdas, podemos comparar las gráficas de la función K. Este debería ser más como CSR:

ppSim = rpoispp(73/(24/400),win=imask)
plot(ppSim)
plot(Kest(ppSim,r=seq(0,.02,len=20)))

patrones de dos puntos en ventanas irregulares

Realmente no puede ver los puntos agrupados en las celdas en el primer patrón, pero si lo traza solo en una ventana gráfica, está claro. Los puntos en el segundo patrón son uniformes dentro de las celdas (y no existen en la región negra) y la función K es claramente diferente de Kpois(r)la función K de CSR para los datos agrupados y similar para los datos uniformes.

Hombre espacial
fuente
2

Además de la publicación de Andy:

Lo que desea calcular es una medida de homogeneidad espacial (por ejemplo, la hipótesis: "¿Están agrupados sus puntos?") Como la función L y K de Ripley .

Esta publicación de blog explica bastante bien cómo hacerlo en R. Según el código descrito, primero etiquetaría cada grupo en su conjunto de datos y luego calcularía en un bucle para cada grupo el sobre crítico a través de K de Ripley

Zarapito
fuente
Actualmente he eliminado mi respuesta. Algunos análisis breves sugirieron que la identificación oportunista de las parcelas basadas en K-medias sesgó las estadísticas locales para estar más agrupadas de lo que sugeriría por casualidad. Esta respuesta aún se aplica aunque +1, (solo la creación de las ventanas a partir de los datos es más problemática de lo que sugeriría mi respuesta original).
Andy W