Determinando se as árvores dentro das lacunas da floresta estão agrupadas usando R?

14

O conjunto de dados em anexo mostra aproximadamente 6000 mudas em aproximadamente 50 aberturas florestais de tamanho variável. Estou interessado em aprender como essas mudas estão crescendo dentro de suas respectivas lacunas (ou seja, agrupadas, aleatórias, dispersas). Como você sabe, uma abordagem tradicional seria executar o I. global de Moran. No entanto, agregações de árvores dentro de agregados de lacunas parecem ser um uso inadequado do I. de Moran. Fiz algumas estatísticas de teste com o I de Moran, usando uma distância de 50 metros, que produziu resultados sem sentido (p-valor = 0,0000000 ...). A interação entre as agregações de gap provavelmente está produzindo esses resultados. Eu considerei criar um script para percorrer lacunas individuais do dossel e determinar o agrupamento dentro de cada lacuna, embora exibir esses resultados ao público seja problemático.

Qual é a melhor abordagem para quantificar o cluster dentro de clusters?

insira a descrição da imagem aqui

Aaron
fonte
1
Aaron, você diz que tentou executar o I de Moran, está interessado em medir como o atributo de um broto se compara aos atributos dos brotos vizinhos (ou seja, você está lidando com um padrão de pontos marcado )? O título parece sugerir que você está interessado apenas na localização das mudas em relação uma à outra e não em seus atributos.
MannyG
@MannyG Sim, só estou interessado em determinar se as mudas estão agrupadas em relação aos locais de outras mudas em qualquer lacuna florestal. Existe apenas uma espécie de interesse e o tamanho das mudas não é de interesse.
Aaron

Respostas:

7

Você não possui um campo aleatório uniforme, portanto, a tentativa de analisar todos os seus dados ao mesmo tempo violará as suposições de qualquer estatística que você escolher lançar no problema. Não está claro em sua postagem se os dados são um processo pontual marcado (ou seja, diâmetro ou altura associado a cada local da árvore). Se esses dados não estão representando um processo pontual marcado, não faço ideia de como você aplicou um I de Moran. Se os dados representarem apenas locais espaciais, eu recomendaria o uso de um Ripley's-K com a transformação Besag-L para padronizar a expectativa nula em zero. Isso permite uma avaliação em várias escalas do cluster. Se seus dados tiverem um valor associado, sua melhor opção será um Moran's-I (LISA) local. Na verdade, eu olhava para isso com as duas estatísticas. Independentemente da sua escolha, você ainda precisará percorrer cada site individual para produzir resultados válidos. Aqui está um exemplo de código R para uma simulação de Monte Carlo do Ripley's-K / Besag's-L usando o conjunto de dados embutido de mudas de sequóias. Deve ser bastante simples modificá-lo para percorrer seus sites e produzir um gráfico para cada um.

# 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
fonte
1
Mas você não pode simplesmente usar o casco convexo como janela para o seu padrão de pontos! Lembre-se de que a janela é a área na qual o padrão que produz os pontos opera. Você sabe a priori que as árvores apenas crescem nessas regiões definidas e você precisa definir sua janela para refletir isso. Você pode atenuar isso definindo o intervalo de K (r) para algo muito pequeno, da ordem do tamanho de 0,3x das suas clareiras, mas obterá resultados tendenciosos devido à falta de correções do efeito de borda. Jeffrey está usando o tamanho de toda a área de estudo para definir seu rmax.
Spacedman
1
No meu exemplo, sim, estou usando a região inteira. É exatamente por isso que eu recomendei percorrer cada local de amostra (lacuna). Cada vez que você define um subconjunto de uma área de amostra específica, executa novamente a análise. Você não pode tratar toda a área de estudo como seu campo aleatório, porque você não possui amostragem contínua. Tendo apenas lacunas amostradas, na verdade você tem gráficos independentes. A função Kest que estou chamando, por padrão, usa uma correção de borda "borda". Existem outras opções de correção de borda disponíveis. Eu diria que a sua unidade experimental é a lacuna do dossel e deve ser analisada como tal.
Jeffrey Evans
1
Ao pensar nisso um pouco mais. Você realmente deve usar polígonos que representam cada espaço como sua janela. Se você configurar seu problema para refletir a unidade experimental, o CSR e K serão enviesados ​​porque a área não reflete o tamanho da lacuna do dossel. Esse é um problema nas recomendações de minhas e de @ Spacedman.
Jeffrey Evans
2
Observe que meu exemplo estendido usou apenas uma grade grossa, porque era uma maneira bastante simples de criar algo com a estrutura certa. Sua máscara deve parecer um mapa de suas áreas de floresta aberta. É tecnicamente errado tentar definir a máscara a partir dos dados!
Spacedman 29/11
1
@Spacedman. Gosto da sua abordagem e é certamente eficiente. Minha preocupação específica é que as lacunas do dossel sejam a unidade experimental. Na sua abordagem, se duas lacunas forem proximais, a largura de banda poderá incluir observações de diferentes unidades de amostragem. Além disso, a estatística resultante não deve refletir o "conjunto" de unidades experimentais, mas deve ser representativa de cada unidade e inferência no processo espacial extraído de padrões comuns entre as unidades experimentais. Se tratado globalmente, representa um processo não-estacionário de intensidade que viola suposições estatísticas.
11139 Jeffrey Evans
4

O que você tem é um padrão de pontos com uma janela que é um número de pequenas regiões poligonais desconectadas.

Você deve poder usar qualquer um dos testes package:spatstatpara o CSR desde que o alimente com uma janela correta. Pode ser um número de conjuntos de pares (x, y) que definem cada limpeza ou uma matriz binária de (0,1) valores no espaço.

Primeiro vamos definir algo que se parece um pouco com seus dados:

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)

e vamos fingir que nossas clareiras são células quadradas que são as seguintes:

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)

Então, podemos traçar a função K desses pontos nessa janela. Esperamos que isso não seja CSR porque os pontos parecem agrupados dentro das células. Observe que eu tenho que alterar o intervalo de distâncias para ser pequeno - da ordem do tamanho da célula - caso contrário, a função K é avaliada por distâncias do tamanho de todo o padrão.

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

Se gerarmos alguns pontos de CSR nas mesmas células, poderemos comparar os gráficos da função K. Este deve ser mais parecido com CSR:

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

dois padrões de pontos em janelas irregulares

Você não pode realmente ver os pontos agrupados nas células no primeiro padrão, mas se você os desenhar por conta própria em uma janela gráfica, ficará claro. Os pontos no segundo padrão são uniformes dentro das células (e não existem na região preta) e a função K é claramente diferente da Kpois(r)função CSR K para os dados em cluster e semelhante para os dados uniformes.

Spacedman
fonte
2

Além do post de Andy:

O que você deseja calcular é uma medida de homogeneidade espacial (logo, a hipótese: "Seus pontos estão agrupados?"), Como as funções L e K de Ripley .

Esta postagem no blog explica muito bem o como fazer em R. Com base no código descrito, eu primeiro rotularia cada cluster no seu conjunto de dados e, em seguida, calcularia em um loop para cada cluster o envelope crítico através do K da Ripley

Maçarico
fonte
No momento, excluí minha resposta. Algumas breves análises sugeriram que identificar oportunisticamente os gráficos com base em médias K fazia com que as estatísticas locais fossem mais agrupadas do que seria sugerido por acaso. Essa resposta ainda se aplica a +1, (apenas a criação das janelas a partir dos dados é mais problemática do que minha resposta original sugeriria).
22712 Andy