Introdução
Neste exercício veremos como tratar a colinearidade entre preditores de um modelo. Lembrando que nos modelos de nicho ecológico, os preditores são geralmente as variáveis bioclimáticas, altitude, biomas, solo, radiação solar, entre outras…
Preparando o ambiente
Neste exercício usaremos os pacotes abaixo. Caso não tenha algum deles basta utilizar o comando: install.packages("nome_do_pacote")
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 ambientais:
Vamos utilizar os dados ambientais do pacote dismo como referência, assim copiaremos desse raster a extenção, resolução e projeção.
#primeiro setamos uma varíavel para captura o caminho dos arquivos no pacote dismo
path <- file.path(system.file(package="dismo"), 'ex')
#depois criamos uma lista neste diretório que tem o padrão grd no final (por isso o sinal de $)
files <- list.files(path, pattern='grd$', full.names=TRUE )
#agora criamos o objeto raster do tipo stack
predictors <- stack(files)Isso deve resultar em um stack de 9 variáveis, com os seguintes nomes: bio1, bio12, bio16, bio17, bio5, bio6, bio7, bio8, biome
Alguns dados de ocorrência para extrairmos os dados
#mais uma vez utilizamos os dados do pacote dismo
file <- paste(system.file(package="dismo"), "/ex/bradypus.csv", sep="")
#agora criamos uma tabela com os dados da preguiça (Bradypus variegatus)
bradypus <- read.table(file, header=TRUE, sep=',')
#Veja a tabela abaixo
head(bradypus)## species lon lat
## 1 Bradypus variegatus -65.4000 -10.3833
## 2 Bradypus variegatus -65.3833 -10.3833
## 3 Bradypus variegatus -65.1333 -16.8000
## 4 Bradypus variegatus -63.6667 -17.4500
## 5 Bradypus variegatus -63.8500 -17.4000
## 6 Bradypus variegatus -64.4167 -16.0000Neste caso, não precisaremos da coluna 1, com o nome da espécie então executaremos o seguinte código:
Vamos olhar em um plot simples como os dados ficam sobre o mapa. Plotando os dados
#carregando vetor com os limites dos continentes
data(wrld_simpl)
#Plotando a primeira variavel para termos um fundo (Opcional)
plot(predictors, 1)
#Plotando os limites
plot(wrld_simpl, add=TRUE)
#Plotando os pontos
points(bradypus, col='blue')
Agora que temos esses dados, vamos extrair os valores de cada ponto
#extraindo os dados ambientais
#Criamos uma variável chamada prevals para abrigar os dados extraídos dos pontos de presença
presvals <- extract(predictors, bradypus)
#Setamos uma valor fixo para a raíz aleatória do R
set.seed(1)
#criamos um vetor de pontos aleatórios sobre a área das variáveis preditoras
backgr <- randomPoints(predictors, 500)
#extraímos os valores a partir desses pontos
absvals <- extract(predictors, backgr)
#criamos um vetor de 1 ou 0, para dizer se um determinado ponto é de presença ou ausência
pb <- c(rep(1, nrow(presvals)), rep(0, nrow(absvals)))
#Agora criamos um dataframe final, combinando os valores dos dados de ocorrência
sdmdata <- data.frame(cbind(pb, rbind(presvals, absvals)))
#Apenas transformamos a coluna bioma em fator
sdmdata[,'biome'] = as.factor(sdmdata[,'biome'])
#Veja o resultado abaixo
head(sdmdata)## pb bio1 bio12 bio16 bio17 bio5 bio6 bio7 bio8 biome
## 1 1 263 1639 724 62 338 191 147 261 1
## 2 1 263 1639 724 62 338 191 147 261 1
## 3 1 253 3624 1547 373 329 150 179 271 1
## 4 1 243 1693 775 186 318 150 168 264 1
## 5 1 243 1693 775 186 318 150 168 264 1
## 6 1 252 2501 1081 280 326 154 172 270 1Olhando a colinearidade
Vamos observar agora a colinearidade desses preditores na área que vamos estudar.
Note que os valores que observaremos dizem respeito tanto aos dados de presença quanto aos dados de background, ou seja mostra de forma geral a colinearidade entre os preditores que pretendemos utilizar no modelo.
Para visualizar a relação é bastante simples. Basta plotar os dados com a função pairs do pacote graphics, nativo do R.
As variáveis que exibem colinearidade são aquelas que apresentam dispersão próximas da linearidade. Ou seja, uma aumenta ou diminui em função da outra. No nosso exemplo, bio12 e bio16, são bons exemplos.
Vejamos essa relação mais de perto.
Note que no caso de bio12 vs bio16 é muito fácil observar a relação entre esses preditores. Contudo, o que dizer da relação entre bio16 vs bio17. Seriam esses correlacionados? qual o limíte aceitável?
Tudo fica mais fácil se quantificarmos a diferença. Índices de correlação de Person podem nos dar uma quantificação da relação entre os preditores. Para isso, faremos o seguinte:
##Medindo a colinearidade
Por conveniência, vamos utilizar um pacote chamado corrplot, que apresenta uma função otimizada para o nosso propósito. Você pode instalar esse pacote, correndo o código comentado.
Uma vez carregado o pacote, vamos à função:
# Para criar a matrix de correlação basta utilizar a função
cors <- cor(presvals)
# checking the table
cors## bio1 bio12 bio16 bio17 bio5 bio6
## bio1 1.00000000 0.05364955 0.06743062 0.04331060 0.94250240 0.91386031
## bio12 0.05364955 1.00000000 0.92643374 0.86053523 -0.08991102 0.18148337
## bio16 0.06743062 0.92643374 1.00000000 0.66297006 -0.04081337 0.17998963
## bio17 0.04331060 0.86053523 0.66297006 1.00000000 -0.08168529 0.11384821
## bio5 0.94250240 -0.08991102 -0.04081337 -0.08168529 1.00000000 0.74778656
## bio6 0.91386031 0.18148337 0.17998963 0.11384821 0.74778656 1.00000000
## bio7 -0.17453449 -0.38418535 -0.32087945 -0.27297842 0.14512737 -0.54821686
## bio8 0.96584412 0.01921980 0.04562234 0.04203550 0.94225305 0.81035244
## biome 0.04027417 -0.25170955 -0.22475186 -0.19086026 0.06770097 -0.02505668
## bio7 bio8 biome
## bio1 -0.17453449 0.96584412 0.04027417
## bio12 -0.38418535 0.01921980 -0.25170955
## bio16 -0.32087945 0.04562234 -0.22475186
## bio17 -0.27297842 0.04203550 -0.19086026
## bio5 0.14512737 0.94225305 0.06770097
## bio6 -0.54821686 0.81035244 -0.02505668
## bio7 1.00000000 -0.02055398 0.12329956
## bio8 -0.02055398 1.00000000 0.06424123
## biome 0.12329956 0.06424123 1.00000000Veja que a matrix é espelhada, ou seja a diagonal superior é a mesma da inferior, bem como a diagonal sempre terá o valor de 1, já que compara a variável com ela mesma.
Caso queira salvar a matrix de correlação use a linha abaixo:
Plotando a correlação
para plotar a correlação vamos criar uma função bastante simples. Ela receberá o argumento cor que será nossa matrix de correlação, e o argumento threshold que dirá qual valor queremos usar como “corte” para salientar a signficância
# Declaramos o nome da função e seus argumentos. Repare que além dos argumentos cor e threshold, os ... nos deixa passar qualquer outro parâmetro para a função corrplot
correlation.plot <- function(cor, threshold, ...){
#aqui dizemos que qualquer valor menor que o threshold será equivalente à 1 no plot
cor[ cor < -threshold] <- -1
#aqui dizemos que qualquer valor maior que o threshold será equivalente à 1 no plot
cor[ cor > threshold] <- 1
#Aqui retornamos à função que gera o plot
return(corrplot(cor, type = "upper", ...))
}
# uma vez criada a função, basta executá-la:
correlation.plot(cors, 0.7)
Compare a matrix dos dados que imprimimos acima com esse plot. Veja que os valores acima de 0.6 foram igualados à 1. Logo, fica nítido quais variáveis têm correlação acima de 0.6 (círculos cheios) e que, portanto, podemos excluir.
Selecionando o set de variáveis.
Selecionar variáveis é sempre crítico. Além dos critérios que comentamos na aula sobre variáveis e seu potencial informativo, é importante manter aquelas com menor colinearidade. Devemos então balancear o potencial informativo das nossas variáveis contra suas autocorrelações e colinearidades.
As linhas abaixo salvam as variáveis ambientais em um novo diretório.
#Capturando o nome das variáveis
vars<-names(predictors)
#Aqui podemos entrar com o nome das variáveis que achamos mais adequadas de acordo com a colinearidade
selected <- c("bio1", "bio12", "bio7", "biome")
# Vamos criar uma pasta para salvar essas variáveis "uncorrelated"
dir.create("uncor") ## Warning in dir.create("uncor"): 'uncor' already exists# Criando um nome para essas variáveis
variable_names <- paste0("uncor/", selected, ".asc")
## escrevendo as variáveis:
#criamos aqui um laço de repetição
#laços tipo `for` repetem-se mudando o valor do primeiro parâmetro (i)
#até que ele bata o tamanho do parâmetro 2, no nosso caso uma sequência de 1 até o tamanho das variáveis selecionadas
#Logo se selecionarmos quatro variáveis o laço rodará quatro vezes e em cada vez, "i" valerá um dos numeros de 1 à 4.
#Se escolhermos cinco variáveis,
for (i in 1:length(selected)) {
#checando se os nomes batem
#se a variável seleciona está relacionada em "vars" então escreve o preditor na pasta "uncor/"
if(selected[[i]] %in% vars){
writeRaster(predictors[[ selected[[i]] ]], filename = variable_names[i], format = "ascii", overwrite = TRUE)
}else{
#Senão avisa qual varível não foi corretamente carregada
cat(paste0("var '",i,"' não encontrada entre as preditoras abaixo: \n "))
cat(vars, "\n Verifique a ortografia \n\n")
}
}Se tudo estiver correto, você deve ter uma nova pasta no seu diretório de trabalho com o nome “uncor” que abrigará as variáveis não correlacionadas.
Espero que o exercício tenha sido útil. Bons estudos.
Prof. Dr. Wilson Frantine-Silva
Pesquisador PNPD - ORCID
Universidade Estadual do Norte Fluminense Darcy Ribeiro - LCA.
Programa de Pós-Graduação em Ecologia e Recursos Naturais