A armadilha da colineridade: como não cair nela usando o VIF (Parte 1).
Análises no R por Luis Macedo-Soares.
Neste post vou falar sobre a colinearidade entre variáveis explicativas e como evitar cair nesta armadilha da busca pelo máximo de explicação da variável resposta, usando o variance inflation factor ou VIF, um teste que pode ser usado para verificar e eliminar a colinearidade da matriz de variáveis explicativas.
Quando vamos realizar uma análise, geralmente procuramos obter o máximo de explicação possível para a hipótese formulada acerca da variável resposta com que estamos trabalhando. Quando possível, medimos e montamos um conjunto de variáveis explicativas com muitas variáveis, pois queremos obter o máximo de explicação. No entanto, algumas vezes, caímos em uma armadilha chamada colinearidade, que é quando as variáveis explicativas possuem uma alta correlação entre si. Isto pode atrapalhar a sua análise e trazer complicações para a interpretação dos resultados.
Antes de realizar a análise, seja uma regressão múltipla, modelos lineares generalizados, ou uma ordenação canônica, é recomendado que você faça um teste com a sua matriz de variáveis explicativas. Neste post, eu vou te apresentar o VIF – variance inflation factor, um teste rápido e de fácil interpretação, para evitar que você caia na armadilha da colinearidade.
Quando vamos realizar uma análise, geralmente procuramos obter o máximo de explicação possível para a hipótese formulada acerca da variável resposta com que estamos trabalhando. Quando possível, medimos e montamos um conjunto de variáveis explicativas com muitas variáveis, pois queremos obter o máximo de explicação. No entanto, algumas vezes, caímos em uma armadilha chamada colinearidade, que é quando as variáveis explicativas possuem uma alta correlação entre si. Isto pode atrapalhar a sua análise e trazer complicações para a interpretação dos resultados.
Antes de realizar a análise, seja uma regressão múltipla, modelos lineares generalizados, ou uma ordenação canônica, é recomendado que você faça um teste com a sua matriz de variáveis explicativas. Neste post, eu vou te apresentar o VIF – variance inflation factor, um teste rápido e de fácil interpretação, para evitar que você caia na armadilha da colinearidade.
Vou apresentar duas funções para você calcular o VIF no R, ambas possuem vantagens e desvantagens, então escolha a que melhor atender as suas necessidades. A primeira delas é a vif.cca() do pacote ‘vegan’, que vamos ver neste post. Para utilizá-la você precisa ter o pacote ‘vegan’ instalado e carregado:
install.packages("vegan")
library(vegan)
Como exemplo vou usar a base de dados Doubs utilizada no livro Numerical Ecology with R de Daniel Borcard e colaboradores (2011). Para mais detalhes sobre a base de dados consulte o livro. Vou importar a planilha de variáveis explicativas ‘env’ e a das espécies ‘spe’. As duas são necessárias na aplicação da função vif.cca().
env <- read.csv("DoubsEnv.csv", header=TRUE, sep=",", dec=".", row.names=1)
spe <- read.csv("DoubsSpe.csv", header=TRUE, sep=",", dec=".", row.names=1)
library(vegan)
Como exemplo vou usar a base de dados Doubs utilizada no livro Numerical Ecology with R de Daniel Borcard e colaboradores (2011). Para mais detalhes sobre a base de dados consulte o livro. Vou importar a planilha de variáveis explicativas ‘env’ e a das espécies ‘spe’. As duas são necessárias na aplicação da função vif.cca().
env <- read.csv("DoubsEnv.csv", header=TRUE, sep=",", dec=".", row.names=1)
spe <- read.csv("DoubsSpe.csv", header=TRUE, sep=",", dec=".", row.names=1)
Seguindo a recomendação do livro, retirei a variável ‘das’ da matriz de variáveis explicativas, salvando-a em um novo objeto. Também criei uma variável categórica a partir da variável declividade ‘pen’, criando um novo objeto ‘env2’.
das <- env[ ,1] env <- env[ ,-1]
pen2 <- rep("very_steep", nrow(env))
pen2[env$pen <= quantile(env$pen)[4]] = "steep"
pen2[env$pen <= quantile(env$pen)[3]] = "moderate"
pen2[env$pen <= quantile(env$pen)[2]] = "low"
pen2 <- factor(pen2, levels=c("low", "moderate", "steep", "very_steep"))
table(pen2)
env2 <- env
env2$pen <- pen2
Com a função str() podemos ver que o objeto ‘env’ possui somente variáveis numéricas, enquanto que em ‘env2’ a variável ‘pen’ é um fator (variável categórica com níveis).
str(env) str(env2)
Para usar a função vif.cca() temos que primeiro realizar uma análise de ordenação canônica, CCA ou RDA por exemplo, pois a função é aplicada em um objeto resultado das funções cca(), rda() ou capscale(). Veja o help da função com ?vif.cca. DESVANTAGEM da função vif.cca(): ela só pode ser aplicada nestes tipos de objetos. Vou fazer a RDA com ‘env’ (onde todas as variáveis são numéricas) e com ‘env2’ (onde ‘pen’ é categórica):
vif.cca(rda1)
alt pen deb pH dur pho nit amm oxy dbo
16.51 1.86 7.02 1.18 2.91 24.13 16.46 30.73 7.10 17.53
rda2 <- rda(spe ~ ., data=env2)
vif.cca(rda2)
alt penmoderate pensteep penvery_steep deb pH dur pho nit
18.84 2.10 2.97 5.50 7.27 1.35 3.05 29.61 17.15
amm oxy dbo
36.18 6.55 15.59
O output da função retorna os valores de VIF para cada variável, e para cada um dos casos temos valores diferentes, uma vez que no primeiro caso temos um valor de VIF para a variável ‘pen’, que é numérica, e no segundo caso temos um valor de VIF para os níveis da variável ‘pen’ que é um fator. VANTAGEM da função vif.cca(): ela pode ser usada com variáveis categóricas.
Segundo os autores do livro Numerical Ecology with R o VIF é uma medida da proporção em que a variância de um coeficiente de regressão é inflacionado pela presença de outra variável explicativa. Para os autores, valores de VIF acima de 20 indicam forte colinearidade e a variável deve ser removida da análise, enquanto valores de VIF acima de 10 devem ser avaliados e evitados se possível. Para Alain Zuur e colaboradores, autores do livro Mixed Effects Models and Extensions in Ecolgy with R (2009) o valor de VIF para remoção da variável é de 3 a 5 (valor usado nos exemplos do livro VIF = 3). A remoção da variável com maior VIF deve ser feita uma de cada vez, e a cada remoção o teste deve ser refeito para verificar se os novos valores de VIF então dentro do limite indicado, escolhido por você.
Voltando aos exemplos, nos dois casos a primeira variável a ser removida, pois apresentou maior VIF, foi ‘amm’ com valores acima de 30. Vamos seguir com o primeiro caso, onde usamos a matriz de variáveis explicativas ‘env’. DESVANTAGEM da função vif.cca(): temos que realizar a análise de ordenação canônica (RDA) novamente, excluindo a variável ‘amm’ para depois aplicar o VIF:
rda1 <- rda(spe ~ alt + pen + deb + pH + dur + pho + nit + oxy + dbo, data=env)
vif.cca(rda1)
alt pen deb pH dur pho nit oxy dbo
14.47 1.77 7.00 1.14 2.88 14.32 11.43 6.50 14.40
Podemos observar a redução do valor do VIF de várias variáveis que provavelmente tinha uma alta correlação com a variável ‘amm’. A variável que mais reduziu o valor do VIF foi ‘pho’, que foi de 24.13 para 14.32. Dependendo do valor de corte adotado, temos que continuar retirando variáveis, e a candidata com maior VIF é ‘alt’.
rda1 <- rda(spe ~ pen + deb + pH + dur + pho + nit + oxy + dbo, data=env)
vif.cca(rda1)
pen deb pH dur pho nit oxy dbo
1.59 2.70 1.14 2.65 9.34 4.32 4.55 10.58
Obtivemos uma redução considerável nos valores de VIF, mas se seguirmos a recomendação de Zuur e colaboradores (2009) temos que continuar com a remoção de variáveis, retirando a variável ‘dbo’.
rda1 <- rda(spe ~ pen + deb + pH + dur + pho + nit + oxy, data=env)
vif.cca(rda1)
pen deb pH dur pho nit oxy
1.54 2.43 1.14 2.42 3.86 3.99 2.39
Chegamos a um resultado que eu considero satisfatório, pois obtivemos valores baixos de VIF mantendo um bom conjunto de variáveis explicativas. Caso você queira obter valores de VIF < 3, continue retirando variáveis uma a uma. Faça isso como um exercício.
No próximo post vamos ver como executar este procedimento usando a função vif() do pacote ‘HH’.
Gostou do conteúdo? Achou relevante? Curta e compartilhe com seus amigos e colegas.
rda1 <- rda(spe ~ pen + deb + pH + dur + pho + nit + oxy + dbo, data=env)
vif.cca(rda1)
pen deb pH dur pho nit oxy dbo
1.59 2.70 1.14 2.65 9.34 4.32 4.55 10.58
Obtivemos uma redução considerável nos valores de VIF, mas se seguirmos a recomendação de Zuur e colaboradores (2009) temos que continuar com a remoção de variáveis, retirando a variável ‘dbo’.
rda1 <- rda(spe ~ pen + deb + pH + dur + pho + nit + oxy, data=env)
vif.cca(rda1)
pen deb pH dur pho nit oxy
1.54 2.43 1.14 2.42 3.86 3.99 2.39
Chegamos a um resultado que eu considero satisfatório, pois obtivemos valores baixos de VIF mantendo um bom conjunto de variáveis explicativas. Caso você queira obter valores de VIF < 3, continue retirando variáveis uma a uma. Faça isso como um exercício.
No próximo post vamos ver como executar este procedimento usando a função vif() do pacote ‘HH’.
Gostou do conteúdo? Achou relevante? Curta e compartilhe com seus amigos e colegas.
Luis Macedo-Soares
Fundador do easy R
Entre o doutorado em Ecologia e o pós-doutorado em Oceanografia Biológica, aprendeu e ensinou a linguagem R, utilizando-a nas diferentes atividades que realizou. Luis é um entusiasta pelo ensino e uso da linguagem R, buscando sempre as formas de aprendizado mais efetivas possíveis para os alunos.