A autocorrelação espacial pode trazer alguns problemas para nossos modelos de nicho, inflando de forma indesejada os modelos e incluindo o ruído proveniente do viés associado.
Nestes casos, podemos utilizar algumas estratégias para escapar desse problema. São basicamentes duas categorias, filtragem e informação de viés. Vamos explorar alguns exemplos e comparações dessas idéias nas próximas páginas
Preparando o ambiente
Carregando pacotes
Seguindo a lógica de programação, é sempre importante carregar antes, pacotes e variáveis que iremos utilizar. Nenhum objeto será reconhecido a não ser que tenha sido definido antes, então:
Setando o diretório
Crie o hábito de setar diretórios onde irá trabalhar. Tudo fica mais fácil de ser encontrado depois.
Carregando dados de ocorrência
Neste exercício usaremos os dados disponíveis no pacote spThin que acabamos de carregar, com occorrências de Heteromys anomalus, um “camundongo de bolso” que também ocorre no norte da América do Sul. Para isso, basta digitar o seguinte:
#Carregando os dados na meória:
#veja na aba environment do seu Rstudio, ou digite o comando ls()
data("Heteromys_anomalus_South_America")
#agora atribuiremos esse dado à um objeto chamado `df`, apenas para ficar mais fácil de trabalhar
df <- Heteromys_anomalus_South_America
#Repare que chamamos o objeto agora sem aspas pois ele já está carregado na memória global do ambiente.
#Com os dois próximos comandos veros e cabeçalho e um resumo das variáveis.
head(df)## SPEC LAT LONG REGION
## 1 anomalus 7.883333 -75.20000 mainland
## 2 anomalus 8.000000 -76.73333 mainland
## 3 anomalus 10.616667 -75.03333 mainland
## 4 anomalus 8.633333 -74.06667 mainland
## 5 anomalus 9.966667 -75.06667 mainland
## 6 anomalus 10.216667 -73.38333 mainland## SPEC LAT LONG REGION
## anomalus:201 Min. : 7.717 Min. :-76.73 mainland:174
## 1st Qu.: 9.500 1st Qu.:-71.70 mar : 2
## Median :10.300 Median :-67.92 tobago : 4
## Mean :10.027 Mean :-68.02 trin : 21
## 3rd Qu.:10.600 3rd Qu.:-63.92
## Max. :11.317 Max. :-60.53Repare que na tabela de ocorrências de Heteromys anomalus há quatro colunas, a primeira com o epíteto específico da espécie, o segundo com a latitude, terceiro longitude e quarto com um atributo chamado região
Note também que o sumário mostra um resumo de cada coluna, de onde podemos rapidamente checar o número de linhas pela coluna SPEC (201), máximos mínimos e medianas, médias e quartis.
### Plotando os dados Carregaremos o “mapa-mundi” do pacote maptools, chamado de wrld_simpl, e depois plotaremos este mapa mais os pointos de ocorrência salvos no objeto df.
## integer(0)Neste plot podemos ver que os pontos estão todos concentrados na porção superior da América do Sul.
Dica: você pode alterar o zoom do plot utilizando os parâmetros
xlimeylimna funçãoplot, ou plotando primeiro os pontos e depois o mapa mundi, assim:
plot(df$LONG, df$LAT)+ plot(wrld_simpl, add=T)
## integer(0)Podemos deixar esse plot mais apresentável com o código abaixo, via ggplot.
shape_df<-fortify(wrld_simpl)#preparamos os dados para o ggplot2
#plotando...
#nessa linha criamos um objeto chamado fig1 o qual sempre chamara esse plot que estamos fazendo agora.
fig1<-ggplot()+
#plotamos o mapa
geom_polygon(data = shape_df , aes(x = long, y = lat, group = group), fill = '#33FF44', color = '#333333')+
#agora os pontos
geom_point(data = df, aes(x = LONG, y = LAT), color="blue")+
#usamos um tema sem linhas ou grades
theme_void()+
#adicionanmos o "mar" como retângulo de fundo azul
theme(panel.background = element_rect('light blue'))
#chamamos fig1
fig1
Para observa a área mais de perto com ggplot, criamos os limítes do plot de acordo com os pontos de plotagem. Note que será útil que a extenção seja um pouco maior que as coordenadas precisas de cada ponto, para garantir que todos os pontos aparecerão no plot sem serem cortados
#objeto para estocar os valores de limite no seguinte
#formato: xmin, xmax, ymin, ymax.
ext<-c(range(df$LONG)+c(-1,+1), range(df$LAT)+c(-1,+1))
#Note que `range` retorna o maximo e mínimo de cada coluna, mas como queremos ampliar os limites um pouco, pediremos ao R para somar -1 ao mínimo e +1 ao máximo de cada coluna.
#Para plotar agora basta o seguinte:
fig1+
coord_equal(xlim=ext[1:2], ylim=ext[3:4])+
theme_minimal()
Agora que temos os pontos e a área devidamante verificadas, vamos as estratégias para eliminar a autocorrelação espacial:
Estratégia 1: arquivo bias (viés)
Nesse estratégia criaremos uma camada dizendo onde o viés se concentra, mas para isso precisaremos de um raster de referência para criar alguns parâmetros. Lembre-se que esse raster deve ser as mesmas camadas que você vai utilizar para a modelagem.
Carregando dados de referência:
Vamos utilizar os dados ambientais como referência, assim copiaremos desse raster a extenção, resolução e projeção.
Caso você tenha os dados na sua máquina
#lista de var ambientais que possuímos
#repare que a função list.files pede("dir", "padrão no nome do arquivo", "caminho comple ou não")
#Aqui utilizei o caminho para o diretório onde estão minhas variáveis, o padrão '.tif' e full.names = verdadeiro.
env.list<-list.files("/media/wilson/MIDGARD/Documents/GIS DataBase/bioclim/2.1/cur/", ".tif", full.names = TRUE)
#aqui pedimos para criar um raster com a primeira variável da lista com a função `raster` do pacote `raster`
#Note que está diferente da aula, pois achei que o processamento na máquina de vocês seria mais eficiente com apenas uma camada e não todas.
env <- stack(env.list)Caso não tenha os dados salvos na sua máquina:
#Note que selecionamos as 12 variáveis bioclimáticas (var ="bio") do worldclim
#com resolução de 10 min (res="10") e já foi feito o download na sua máquina
#(download = TRUE). Agora elas devem estar no seu diretório de trabalho dentro
#da pasta wc10
env <- getData("worldclim", var="bio", res="10", download = TRUE)
###Vamos olhar o objeto
### Cortando as camadas
Este passo pode ser opcional. Aqui neste exercício, note que estou trabalhando com a cama a nível mundial, pois não havia cortado ela para alguma área específica. Então vamos agora cortar essa camada para a área de interesse, neste caso, onde estão nossos pontos.
#vamos utilizar a função `crop` do pacote `raster`
#note que ext é o objeto que criamos a pouco
#que tem a exteção dos pontos de ocorrencia + 1°
env<-crop(env, ext)
#plotando para checar o raster e os pontos
plot(env[[1]])+points(df$LONG, df$LAT)
## integer(0)IMPORTANTE: a extenção que devemos cortar deve ser a mesma das camadas que iremos utilizar para a modelagem. Nesse ponto, vou utilizar um recore quadrado um pouco maior que os pontos de ocorrência, mas há outros métodos mais adequados para delimitar a área de modelagem. Apenas tenha em mente!
transformando os pontos em raster
Nosso arquivo de bias é será um raster, logo precisamos transformar nossos pontos em raster para ter um ponto de saída.
Nas próximas linhas vamos então:
1. criar um objeto apenas com lat e long 2. criar um raster com a função rasterize 3. plotar para checar se está tudo Ok.
#1. criando o objeto occ
occ<- df[,c("LONG", "LAT")]
#2. usando a funçõa rasterize. Note que a função pede:
#a) os pontos; b) um objeto de referência; c) um valor para dar ao pontos
occur.ras <- rasterize(occ, env, 1)
#3 plotando
plot(occur.ras, col="blue", legend=F)
Agora vamos extrair alguns dados dessa camada que acabamos de criar:
#nessa primeira linha, pegamos as células que tem
#valor = 1, ou seja as células que correspondem aos
#pontos de amostragem
presences <- which(values(occur.ras)==1)
#Aqui pegamos as coordenadas x e y de cada célula
pres.locs <- coordinates(occur.ras)[presences,]
#veja como eles se parecem
head(presences)## [1] 675 676 678 703 709 757## x y
## [1,] -74.25000 11.25
## [2,] -74.08333 11.25
## [3,] -73.75000 11.25
## [4,] -69.58333 11.25
## [5,] -68.58333 11.25
## [6,] -60.58333 11.25Criando mapa de densidade
Agora vamos criar um raster com densidade de pontos a partir dos dados que extraímos, utilizando a função kde2d do pacote MASS.
dens <- kde2d(
pres.locs[,1],
pres.locs[,2],
n= c( nrow(occur.ras), ncol(occur.ras) ),
lims = c(extent(env)[1], extent(env)[2],
extent(env)[3], extent(env)[4])
)
str(dens)## List of 3
## $ x: num [1:34] -77.7 -77.1 -76.6 -76 -75.5 ...
## $ y: num [1:109] 6.67 6.72 6.77 6.82 6.88 ...
## $ z: num [1:34, 1:109] 1.19e-05 1.64e-05 2.18e-05 2.83e-05 3.61e-05 ...Note que essa função cria um objeto com três listas de tamanho n. As listas x e y gardam os valores das coordenadas de cada ponto estimado, enquanto a lista z guarda os valores de densidade de cada pixel
Agora vamos criar o raster dessa densidade
#criamos o raster usando os valores do objeto dens
#e os valores de env
dens.ras <- raster(dens, env)
#Agoras colocamos dens.ras na mesma resolução das camadas ambientais
dens.ras2 <- resample(dens.ras, env)
plot(dens.ras2) 
Note que o raster de densidade é varia de acordo com a concentraçõ de pontos que tinhamos de Heteromys anomalus no iníncio. Essa camada pode então ser utilizada como informação para o MaxEnt ou outros algortmos para amostrar mais o background nesse região do que em outras.
Para salvar o raster podemos utilizar a seguinte linha:
## class : RasterLayer
## dimensions : 34, 109, 3706 (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667 (x, y)
## extent : -77.66667, -59.5, 6.666667, 12.33333 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : /media/wilson/MIDGARD/Documents/Academico/AARE - MNE/Aula 8 - Autocorrelacao/biasfile.grd
## names : layer
## values : 1.173139e-09, 0.0403418 (min, max)Estratégia 2: Filtrando por distância mínima
Nessa estratégia, ao invés de informarmos o viés da amostragem, tentaremos removê-la. Para isso utilizaremos a aleatorização do pacote spThin, utiliando os mesmos dados que já temos. Vamos realizar a filtrangem tanto no espaço geográfico quanto no espaço ambiental.
Filtro Geográfico
Nessa estratégia, vamos setar um valor mínimo de distância em quilômetros os quais nossos pontos devem estar uns dos outros.
thinned_data <-
thin(
loc.data = df, #dataframe com os dados de ocorrência
lat.col = "LAT", #coluna com info de latitude
long.col = "LONG", #coluna com info de longitude
spec.col = "SPEC", #coluna com nome da espécie
thin.par = 10, #distância mínima entre pontos (em Km)
reps = 100, # numero de repetições
locs.thinned.list.return = TRUE, #retornar a lista com os dataframes selecionados
write.files = FALSE,
write.log.file = FALSE
)## **********************************************
## Beginning Spatial Thinning.
## Script Started at: Mon May 3 20:39:46 2021
## lat.long.thin.count
## 122 123 124
## 11 46 43
## [1] "Maximum number of records after thinning: 124"
## [1] "Number of data.frames with max records: 43"
## [1] "No files written for this run."Note que o pacote retornou uma lista com três soluções possíves, selecionando 13 datasets com 122 ocorrências, 47 datasets com 123 e 40 com 124 ocorrências.
Podemos olhar a tabela da primeira solução selecionada e plotá-la para ver como ficou:
## Longitude Latitude
## 1 -75.20000 7.883333
## 2 -76.73333 8.000000
## 3 -75.03333 10.616667
## 4 -74.06667 8.633333
## 5 -75.06667 9.966667
## 6 -73.38333 10.216667plot(thinned_data[[1]], col="blue", pch=16)
plot(wrld_simpl, add=T)
points(df[,3:2], col="red", pch=3)
Os pontos azuis preenchidos (pch=16, col ="blue") correnspondem aos dados selecionados, enquanto os pontos representados por cruzes vermelhas (pch=3, col="red") representam os pontos excluídos.
Nesta escala é difícil de ver os pontos excluídos. É importante lembrar que próximo ao equador cada grau decimal (e.g. -65, -66) equivale à cerca de 111 km, 60 arc minutos ou 60 milhas náuticas. Ou seja, como estamos observando um mapa com cerca de 16,2 graus de longitude entre as latitudes 7 e 11, logo devemos estar observando uma distância de ~1700 km. Isso torna bastante difícil visualizar as diferenças de apenas 10 km, parâmetro que utilizamos para selecionar os pontos.
Então vamos olhar mais de perto entre as longitudes -68.5 e -67.5 e latitudes 10.1, 10.6, onde a extenção total no eixo x deve corresponder a cerca de 100 km:
plot(thinned_data[[1]],
xlim=c(-68.5,-67.5),
ylim=c(10.1,10.6))
plot(wrld_simpl,
col="#FFFFAA",
add=T)
points(df[,3:2], pch=3, col="red", cex=1.1)
points(thinned_data[[1]], col="blue", pch=16)
Note agora que fica muito mais fácil observar quais pontos foram excluídos e como eles são próximos uns aos outros.
Estratégia 2.1: filtrar no espaço ambiental (bonus)
Como vimos na aula é possível filtrar os dados no espaço ambiental. Neste caso, assumiríamos que a proxímidade dos pontos naquele ambiente corresponderia a não independência dos dados, e de certo modo, uma repetição não aleatória da amostragem do nicho fundamental.
Neste caso, podemos utilizar a mesma lógica de aleatorização da função thin do pacote spThin com a diferença de fornecer um dataset com valores de coordenadas ambientais ao invés de latitudes e longitudes.
Nesta etapa do exercício iremos basicamente:
1. Extrair dados ambientais a partir dos pontos de ocorrência 2. Reduzir o espaço via PCA 3. Avaliar um threshold 4. filtrar dados ambientais 5. observar os dados no espaço geográfico a partir do filtro ambiental
1. Extrair dados ambientais
O primeiro passo a aqui é extrair os dados ambientais a partir dos nossos pontos de ocorrência
#extraindo dados ambientais e criando um data.frame para a função `thin`
env.values<-extract(env, occ)
head(env.values)## bio1 bio2 bio3 bio4 bio5 bio6 bio7 bio8 bio9 bio10 bio11 bio12 bio13 bio14
## [1,] 279 101 84 397 344 224 120 276 283 284 274 2560 364 26
## [2,] 264 78 89 218 308 221 87 261 265 266 261 2414 298 80
## [3,] 278 97 82 405 336 218 118 277 276 282 273 1114 198 4
## [4,] 276 104 84 450 341 218 123 272 275 281 270 2020 347 15
## [5,] 264 113 77 463 337 191 146 262 260 269 257 1336 200 25
## [6,] 281 118 80 593 356 210 146 276 279 287 272 1308 204 15
## bio15 bio16 bio17 bio18 bio19
## [1,] 58 1055 135 299 1009
## [2,] 34 770 276 422 770
## [3,] 71 517 21 188 168
## [4,] 62 859 82 270 821
## [5,] 52 503 103 355 374
## [6,] 58 517 59 187 357Podemos plotar esses dados:
Neste gráfico, visualizamos a relação dos dados extraídos a partir de todos os nossos pontos de ocorrência em um espaço ambiental composto por apenas duas variáveis bio1: temperatura e bio12: precipitação
Note que há um valor de alpha na cor preta nos pontos do gráfico, ou seja, pontos mais escuros são sobreposições nesse espaço, logo são pontos que gostaríamos de remover.
Um aspecto muito importante é que estamos olhando apenas duas variáveis nesse plot, e os pontos de ocorrência podem ter outra relação quando considerarmos outras variáveis. Para tentar lidar com esse problema, iremos construir um espaço ambiental via PCA. Assim, estaremos olhando não mais para cada variável por vez, mas para os componentes que reduzem a variância entre elas. Vejamos isso de forma mais concreta a seguir.
2. Reduzindo as dimensões com uma PCA
Gerar PCAs no R é simples, basta utilizar a função prcomp. Há vários pacotes dedicados à isso, mas para o nosso propósito aqui, essa função será mais que suficiente.
Na linha a seguir, processamos um PCA a partir dos dados ambientais (env.values) corrigindo as coordenadas dos pontos no centro de variância (center = T) e bem como sua escala (scale. = T). Além disso, sumarizamos os resultados com a função summary()
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 2.7069 2.3057 1.5555 1.36112 1.00070 0.77038 0.45284
## Proportion of Variance 0.3856 0.2798 0.1273 0.09751 0.05271 0.03124 0.01079
## Cumulative Proportion 0.3856 0.6655 0.7928 0.89031 0.94301 0.97425 0.98504
## PC8 PC9 PC10 PC11 PC12 PC13 PC14
## Standard deviation 0.40720 0.25090 0.15881 0.10816 0.09217 0.07164 0.04535
## Proportion of Variance 0.00873 0.00331 0.00133 0.00062 0.00045 0.00027 0.00011
## Cumulative Proportion 0.99377 0.99708 0.99841 0.99903 0.99947 0.99974 0.99985
## PC15 PC16 PC17 PC18 PC19
## Standard deviation 0.03766 0.03310 0.01349 0.01093 5.882e-16
## Proportion of Variance 0.00007 0.00006 0.00001 0.00001 0.000e+00
## Cumulative Proportion 0.99993 0.99998 0.99999 1.00000 1.000e+00Note que a variância acumulada (Cumulative Proportion) aumenta a medida que os componentes são adicionados. Neste exemplo, vemos que os cinco primeiros componentes já retém mais de 94% da variância. Com isso, sabemos que é complemente possível continuar a análise considerando apenas esses primeiros cinco componentes.
Agora vamos plotar as coordenadas dos nossos pontos nos dois primeiros componentes (“66% da variância”) com: plot(pca$x). Note ainda que mudei apenas as cor para preto rgb(0,0,0) com transparência (“alfa”) de 30% (o último 0 do rgb). A forma é um círculo preenchido (pch=20) e o tamanho é uma vez e meia (cex = 1.5)
Interessantemente temos mais uma vez pontos mais escuros que significam sobreposição. Um outro detalhe é que, olhando a disposição dos pontos, fica difícil dizer qual a distância mínima que um ponto deve estar de outro para decidirmos como filtrá-los. Neste caso, precismos identificar um threshold.
3. Econtrando um threshold
Nas proíxmas linhas vamos: 1. criar um data frame apenas com os componentes que iremos utilizar (no nosso cas de 1 a 5); 2. vamos calcular a distância euclidiana entre os pontos 3. vamos criar um vetor apenas com as menores distâncias de cada ponto 4. vamos observar quantís das menores distâncias para saber quais são os thresholds e quanto perderemos em cada escolha
#1. criando o data frame para os cinco primeiros componentes
espacoAmbiental<-as.data.frame(pca$x[,c(1:2)])
#2. nessa linha criamos uma matríz das distâncias euclidianas
pcaDistMatrix<-as.matrix(dist(espacoAmbiental, method = "euclidian"))
#2.1 Aqui apenas removemos os zeros da comparação de um ponto com ele mesmo
pcaDistMatrix[pcaDistMatrix==0]<-NA
#3. Neste passo utilizamos a função `apply` no objeto pcaDistMatrix
#apenas nas linhas (sungundo argumento = 1), com a função min, e removendo quaisquer `NA`
DistMinimaPorAmostra<-
apply(as.data.frame(pcaDistMatrix),1, min, na.rm=T)
#4. Aqui criamos um vetor que armazena essas distâncias separadas em quantis de 1% saindo de 1 até 100% a cada 1%
filter.par<-quantile(DistMinimaPorAmostra, probs = seq(0.01, 1, by = 0.01))
#4.1 Veja como ele se parecem os primeiros 10 elementos
head(filter.par, 10)## 1% 2% 3% 4% 5% 6% 7%
## 0.05116354 0.05116354 0.05116354 0.05116354 0.07655885 0.07751782 0.09751995
## 8% 9% 10%
## 0.09751995 0.11870406 0.12237304Esse vetor nos dá a distribuição das frequências das menores distâncias euclidianas entre os pontos. Em outras palavras, podemos dizer que 1% das menores distâncias são menores ou iguais à 0.05116354 ou 10% das distâncias são menores ou iguais à 0.0.12237304. Na prática, “em situações normais de temperatura e pressão”, isso significa que filtraremos ~10% das amostras se escolhermos um threshold de ~0.12.
4. filtrar dados ambientais
Neste passo vamos enfim filtrar os dados. Para isso utilizaremos uma adaptação da função thin do pacote spThin que filtrará não mais por distância geográfica, mas por distância euclidiana entre as coordenadas dos pontos da PCA.
Para importarmos a função, vamos utilizar a função source e a sintaxe thin.env para chamar a função:
#importando o arquivo com a função customizada
source(file = "https://raw.githubusercontent.com/wilsonfrantine/ENM101/main/R/thin.env.r")Agora a função thin.env deve estar na memória local (veja a aba environment do Rstudio ou digite ls() no terminal) podemos utilizá-la:
#executando a função e salvando no objeto env.filtered.
#Note que passamos para a função apenas o data.frame com os valores da PCA (espacoAmbiental);
#o quinto elemento dos nossos quantis ~ 0.6
#e pedimos para repetir 10 vezes
env.filtered<-thin.env(espacoAmbiental, filter.par[[50]], 10)
#veja quantos pontos haviam e quantos restaram
cat(paste("N no dataset original:",nrow(espacoAmbiental),"\n"))## N no dataset original: 201## N no dataset filtrado: 94Vemos uma redução importante no número de pontos utilizando um threshold relativamente baixo.
SUGESTÃO DE EXERCÍCIO:
experimente filtrar os dados com diferentes thresholdsthin.par[[1]],thin.par[[2]],thin.par[[50]]… Note quantos pontos são perdidos à medida que se aumenta o threshold.
Vamos plotar os dados para ver como eles ficam. O código abaixo gera um plot para os dois primeiros componentes. Lembre-se que sempre é possível observar mais componentes em uma PCA. No nosso caso, lembre-se de voltar atrás e definir a variável espaçoAmbiental como pca$x[,c(1:3)] ou pca$x[,c(1:4)] para pegar os três ou quatro primeiros componentes, respectivamente.
plot(espacoAmbiental[,c(1,2)], pch=20, cex=1.5, col="grey70")
points(env.filtered[[1]][,c(1,2)], pch=20, cex=1.5, col="blue")
legend("bottomright",
legend = c("original","filtrado"),
pch = c(20,20),
cex = c(1.5,1.5),
col = c("grey70","blue") )
Uma vez que estamos satisfeitos com a filtragem, podemos salvar as coordenadas lat e long dos pontos de entrada. Para isso, iremos pegar as linhas que restaram após a filtragem e utilizá-las como índice para selecionar os pontos de ocorrência:
#Primeiro criamos um objeto que guarda quais linhas foram selecionadas no primeiro
#dataset do objeto env.filtered
selected<-as.numeric(row.names(env.filtered[[1]]))
#agora simplesmente dizesmo que queremos as ocorrências que batem com essa lista
occ.env.filtered<-occ[selected,]
#vamos ver como ficou?
head(occ.env.filtered)## LONG LAT
## 1 -75.20000 7.883333
## 2 -76.73333 8.000000
## 4 -74.06667 8.633333
## 5 -75.06667 9.966667
## 6 -73.38333 10.216667
## 8 -73.48333 10.3666675. Observando o filtro ambiental no espaço geográfico
Uma vez que nossos dados foram selecionados independentemente do espaço geográfico, podemos ver como eles se comportam.
Vamos plotar o objeto occ.env.filtered:
base<- ggplot()+
geom_polygon(data = shape_df , aes(x = long, y = lat, group = group), fill = '#33FF44', color = '#333333')+
theme_void()+
theme(panel.background = element_rect('light blue'))+
coord_equal(xlim=ext[1:2], ylim=ext[3:4])
colors<-c("Sem filtro"="grey40","Ambiental"="blue", "Geográfico"="red")
base+
geom_point(data = occ, aes(x = LONG, y = LAT, col="Sem filtro"), size=2)+
geom_point(data = occ.env.filtered, aes(x = LONG, y = LAT, col="Ambiental"), size=2)+
scale_color_manual(values=colors, name="Filtros")
ggplot()+
geom_point(data=espacoAmbiental, aes(PC1, PC2, col="Sem filtro"), size=2)+
geom_point(data=env.filtered[[1]], aes(PC1, PC2, col="Ambiental"))+
scale_color_manual(values=colors, name="Filtros")
#geom_point(data = thinned_data[[1]], aes(x = Longitude, y = Latitude), color=rgb(1,0,0,0.3), cex=2)+Por fim, podemos comparar filtros distintos. Primeiro no espaço geográfico:
#Plot geral
base+
geom_point(data = occ, aes(x = LONG, y = LAT, col="Sem filtro"), size=2)+
geom_point(data = thinned_data[[1]],
aes(x = Longitude, y = Latitude, col="Geográfico"), cex=2)+
geom_point(data = occ.env.filtered,
aes(x = LONG, y = LAT, col="Ambiental"), size=2)+
scale_color_manual(values=colors, name="Filtros")
Vemos uma boa disperssão dos dois filtros. Como há vários pontos sobrepostos para ambos, vamos ver mais de perto:
#Plot de um zoom
base+
geom_point(data = occ, aes(x = LONG, y = LAT, col="Sem filtro"), size=2)+
geom_point(data = thinned_data[[1]],
aes(x = Longitude, y = Latitude, col="Geográfico"), cex=2)+
geom_point(data = occ.env.filtered,
aes(x = LONG, y = LAT, col="Ambiental"), size=2)+
scale_color_manual(values=colors, name="Filtros")+
coord_equal(xlim=c(-68.5,-67.5), ylim=c(10.1,10.6))
Aqui vemos como os dois métodos selecionaram de modo distintos os ambientes. Neste quadro, as duas abordagens parecem ter amostrado bem o espaço geográfico.
Vamos checar os filtros no Espaço Ambiental:
#Plot do espaço ambiental e seus filtros
ggplot()+
geom_point(data=espacoAmbiental, aes(PC1, PC2, col="Sem filtro"))+
geom_point(data=env.filtered[[1]], aes(PC1, PC2, col="Ambiental"))+
geom_point(data=espacoAmbiental[row.names(thinned_data[[1]]),], aes(PC1, PC2, col="Geográfico"))+
scale_color_manual(values=colors, name="Filtros")