simulacao_pontos_no_espaco.Rmd

---
title: "Modelos de distribuição espacial"
subtitle: "Simulação de pontos no espaço"
author: "Exercício prático da disciplina Ecologia de Populações de Comunidades Vegetais, IB-USP"
output:
  html_document:
    toc: true
    theme: united
  pdf_document:
    toc: true
    highlight: zenburn
---

```{r settings, echo=FALSE,  message=FALSE}
library(knitr)
knitr::opts_chunk$set(echo = TRUE, eval = FALSE)
library(spatstat)
## descomente para gerar com todos os outputs
##knitr::opts_chunk$set(echo = TRUE, eval = TRUE)
```

# Introdução

Um modelo é uma descrição de um fenômeno, que preserva alguns dos
aspectos relevantes deste fenômeno. Um dos tipos de modelos que
usamos para descrever a estrutura espacial de populações de plantas
são algorítimos que distribuem pontos num plano, segundo um conjunto
de regras. Estes são os modelos de **distribuição de pontos no
espaço**, ou **modelos de padrões de pontos no espaço** [^4].

Aqui vamos entender as regras dos quatro modelos de padrões de pontos
no espaço usados no artigo de Linn *et al.* (2011), usando a
linguagem de programação R.

## Preparação para o exercício

Abra o programa R, clicando no ícone que está na área de trabalho do
seu computador: ![](images/logor.png)

Se tudo deu certo até aqui, abrirá uma janela do R como essa:

![](images/tela_inicial_r.png)

O símbolo ">", circundado em azul na imagem, indica o início da linha
de comando ou **prompt**, onde você deve escrever comandos para o R.


Copie e cole o comando abaixo na linha de comando do R, para carregar
o pacote (biblioteca de funções) que vamos usar:

```{r }
library(spatstat)
```

Para realizar este roteiro, basta seguir as instruções de cada
seção. Na maioria dos casos, você terá apenas que copiar o colar
comandos, e interpretar os resultados.  Em alguns casos, você terá que
modificar um pouco alguns parâmetros dos comandos. Isto estará
explicado em cada uma das seções seguintes.


# Alguns modelos de distribuição de pontos em um plano

## Poisson homogêneo

O primeiro modelo que vamos usar é chamado **Poisson homogêneo**, e
descreve pontos distribuídos ao acaso em uma área. Execute os comandos
abaixo no R, para simular um padrão de pontos gerados por este
processo Poisson Homogêneo, em uma área quadrada unitária (ou seja, de
dimensões 1x1) [^2] 

```{r  Poisson homogeneo}
## Gera o padrão de um processo PH, com intensidade = 100
PH <- rpoispp(lambda = 100, win = unit.square())
## Plota o padrão de pontos
plot(PH, main = "Poisson homogêneo")
```

A função que simula o processo chama-se `rpoispp`. Ela tem o argumento
`lambda`, que é a **intensidade** do processo. A intensidade aqui
corresponde à densidade média de pontos por unidade de área. Depois de
gerar os pontos, o código acima usa a função `plot` para fazer um
gráfico do padrão de pontos. Experimente executar os comandos acima
com outros valores de intensidade, para verificar o que muda no padrão
de pontos gerados pelo modelo.


Como a intensidade do padrão que geramos inicialmente é 100 e a área é
um quadrado de 1x1, o valor **esperado** de pontos nesta área é
de 100. Isso não quer dizer que em toda simulação teremos exatamente
esta quantidade de pontos, e sim que o **número médio** de pontos
é 100. Ou seja, se repetirmos várias vezes a simulação deste processo,
a média das densidades de todas as simulações corresponde à
intensidade do processo.

A figura abaixo mostra seis repetições do padrão Poisson Homogêneo com
intensidade de cem. O número de pontos em cada simulação está escrito
acima de cada gráfico.

```{r  Poisson homogeneo media, echo = FALSE,  eval=TRUE}
## seis  simulações do modelo
set.seed(3)
HP1 <- rpoispp(lambda = 100, win = unit.square(), nsim = 6)
names(HP1) <- sapply(HP1, function(x) paste0("N = ", x$n))
plot(HP1, multiplot=TRUE, main = "")
```

## Modelo Poisson não-homogêneo

Este é um modelo de distribuição de pontos ao acaso, mas com uma
intensidade que varia no espaço. Para exemplificar, vamos imaginar que
em metade da nossa área quadrada de 1x1 a intensidade do processo é
$\lambda = 20$, e na outra metade a intensidade é de $\lambda =
180$. Execute os comandos abaixo para simular um padrão Poisson não
homogêneo com estas intensidades:

```{r  Poisson nao homogeneo}
## Aqui definimos uma função que descreve como a intensidade muda na área
## No caso, a intensidade muda de 180 para 20 a partir do valor da coordenada X = 0.5
padrao1 <- function(x,y) ifelse(test = x > 0.5, yes = 20, no = 180)
## Gera o padrão de pontos
IP <- rpoispp(lambda = padrao1, win = unit.square())
## Plota o gráfico
plot(IP, main = "Poisson não-homogêneo")
## Uma linha para marcar as duas regiões
abline(v = 0.5, lty =2 , col = "blue")
```

Um dos usos do modelo Poisson não-homogêneo é descrever a variação na
densidade de plantas sob diferentes condições ambientais. Um exemplo
pode ser a maior densidade de uma espécie de planta em áreas de menor
altitude. Para ilustrar, vamos usar o mapa de altitude do terreno na
[parcela permanente de Barro Colorado](http://ctfs.si.edu/webatlas/datasets/bci/), 
no Panamá. Esta é uma área de 500 x 1000m, na qual todas as árvores acima de 1 cm de
diâmetro são monitoradas desde 1983.

O comando abaixo mostra o mapa de altitude da parcela. A escala de cor
mostra a altitude acima do nível do mar, em metros.

```{r  BCI altitude}
image(bei.extra$elev, main = "Elevação do terreno em BCI")
```

Agora vamos criar uma função que descreve um decréscimo da intensidade
do processo Poisson com altitude, de forma exponencial. Os números
foram escolhidos para obtermos uma quantidade razoável de pontos neste
exemplo. Os comandos abaixo calculam os valores de intensidade para a
altitude de cada região da parcela:

```{r  BCI padrao ficticio}
## Valores de intensidade em função da elevação na parcela
padrao2 <- exp(20 - bei.extra$elev/10)*5e-6
```

Ou seja, criamos um padrão fictício de variação da intensidade, que
decai exponencialmente emfunção da altitude. A relação entre
intensidade do processo e altitude segundo este padrão fictício está
no gráfico abaixo:

```{r  BCI  relacao intensidade x altitude, eval=TRUE, echo=FALSE}
padrao2 <- exp(20 - bei.extra$elev/10)*5e-6
plot(bei.extra$elev$v, I(padrao2$v*25),
     main = "Padrão fictício de intensidade em BCI",
     xlab = "Altitude em relação ao nível do mar (m)",
     ylab = "Intensidade por m2")
```


Agora vamos simular um padrão de pontos com o modelo Poisson
não-homogêneo, com esta função de intensidade. Com isso, simulamos uma
espécie fictícia que tem maiores densidades nas áreas mais baixas:

```{r  BCI planta ficticia}
## Gera padrão de pontos com a intensidade definida pelo padrão
bci.pp <- rpoispp(lambda = padrao2)
## Plota o mapa do terreno e os pontos
image(bei.extra$elev, main = "Elevação do terreno em BCI")
points(bci.pp)
```

## Modelo Thomas Homogêneo

No modelo de Thomas Homogêneo, os pontos são distribuídos em
agregados, ou *clusters*. A posição de cada agregado é definida por uma
coordenada central, em torno da qual se distribuem pontos, até uma
certa distância do centro, e com uma certa intensidade. As coordenadas
centrais dos agregados são distribuídas por um processo Poisson
Homogêneo. Os pontos em torno de cada coordenada central são
distribuídos com uma intensidade que decresce com a distância ao
centro do agregado[^3].

Execute o comando a seguir para simular um processo de Thomas homogêneo:

```{r  Thomas homogeneo}
TH <- rThomas(kappa = 10, scale = 0.01, mu = 10, win = unit.square())
plot(TH, main = "Thomas homogêneo")
```

O parâmetro `mu` da função `rThomas` define a intensidade de pontos em
cada agrupamento. Veja o que acontece se triplicamos esta intensidade:

```{r  Thomas homogeneo varia mu}
## Nova simulacao com mu = 30
TH1b <- rThomas(kappa = 10, scale = 0.01, mu = 30, win = unit.square())
## Plota o grafico anterior junto com o grafico do novo padrao
par(mfrow=c(1,2)) ## abre janela para dois graficos juntos
plot(TH, main = "mu original")
plot(TH1b, main = "mu alterado")
par(mfrow = c(1,1))
```

Já o parâmetro `kappa` da função define a intensidade de agrupamentos,
ou seja, a densidade esperada de agregados (ou a densidade esperada de
pontos imaginários que definem o centro de cada agrupamento). Veja o
que acontece se triplicamos esta intensidade:

```{r  Thomas homogeneo varia kappa}
## Nova simulacao com kappa = 20
TH2 <- rThomas(kappa = 30, scale = 0.01, mu = 10, win = unit.square())
## Plota o grafico anterior junto com o grafico do novo padrao
par(mfrow=c(1,2)) ## abre janela para dois graficos juntos
plot(TH, main = "Kappa original")
plot(TH2, main = "Kappa alterado")
par(mfrow = c(1,1))
```

O terceiro parâmetro do processo de Thomas homogêneo é chamado `scale` no R. 
Use os comandos abaixo para umentar o valor de `scale` em
relação à primeira simulação, e descobrir o que isso causa:

```{r  Thomas homogeneo varia scale}
## Nova simulacao com scale = 0.03
TH3 <- rThomas(kappa = 10, scale = 0.03, mu = 10, win = unit.square())
## Plota o grafico anterior junto com o grafico do novo padrao
par(mfrow=c(1,2)) ## abre janela para dois graficos juntos
plot(TH, main = "scale original")
plot(TH3, main = "scale alterado")
par(mfrow = c(1,1))
```

## Modelo de Thomas não-homogêneo

No modelo de Thomas não-homogêneo, a intensidade de pontos **nos
agrupamentos** varia no espaço. Podemos simular este processo no R com
a mesma função que usamos acima, mas indicando um padrão de variação
para o argumento `mu`. 

Os comandos a seguir criam um padrão com duas
zonas de diferentes intensidades, em uma área quadrada unitária. Em
seguida, simula um processo de Thomas não homogêneo, em que a
intensidade de pontos nos aglomerados segue este padrão de
intensidade:


```{r  Thomas  nao homogeneo}
## Cria o padrão de intensidade de agregados: duas zonas com intensidades 5 e 20
padrao3 <- function(x,y) ifelse(x > 0.5, 5, 20)
## Simula o processo de Thomas com intensidade de pontos no cluster igual ao padrão acima
IT <- rThomas(kappa = 20, scale = 0.01, mu = padrao3, win = unit.square())
## Plota os pontos e a linha que divide as duas zonas de intensidade
plot(IT)
abline(v = 0.5, lty =2 , col = "blue")
```

# Referências bibliográficas

Lin, Y. C., Chang, L. W., Yang, K. C., Wang, H. H., & Sun,
    I. F. (2011). Point patterns of tree distribution determined by
    habitat heterogeneity and dispersal limitation. Oecologia, 165(1),
    175-184.


[^2]: As linhas precedidas de `##` são comentários para
explicar o que cada comando faz. 

[^3]: Mais precisamente, a distribuição dos pontos em cada agregado
    segue uma normal bivariada, centrada na coordenada
    central. Imagine um sino sobre cada coordenada central, que
    descreve a intensidade em torno desta coordenada. É justamente a
    menor probabilidade de pontos distantes do centro que criam o
    agregado.

[^4]: Em inglês *point process models* e *point pattern models*.