Há algum tempo venho analizando um grupo de dados coletados em uma lagoa do Litoral do Rio Grande do Sul com o intuito de modelar a distribuição de moluscos limnicos. Venho usando o software R para estas análises. Um dia posto aqui os prós e os contras de usar o programa R, apesar de muita gente achar que o R é impossível.
Como uma primeira abordagem decidi investigar a colinearidade das variáveis que disponho para o ajuste de um modelo com o menor número possível de variáveis. Dando uma buscada no blog R-Bloggers descrobri uma função criada por Stephen Turner (http://gettinggeneticsdone.blogspot.com.br/2012/08/more-on-exploring-correlations-in-r.html) que coloca informações relevantes a cerca dos dados que disponho como o coeficiente de correlação e o valor de p(asteriscos) em um formato de matriz com todas as variáveis que quero analizar. Bem útil para ter uma visão da colinearidade dos preditores e o uso das funções são simples, sem problemas. Esse é um exemplo que o ambiente R e um ferramenta para a criatividade. Para mais detalhes acesse o post original. Abaixo reproduzo o código das funções para serem rodadas no R, usadas e aprovadas.
Função que cria a matriz de correlação com o valor de p.
cor.prob <- function(X, dfr = nrow(X) - 2) {
R <- cor(X, use = "pairwise.complete.obs")
above <- row(R) < col(R)
r2 <- R[above]^2
Fstat <- r2 * dfr/(1 - r2)
R[above] <- 1 - pf(Fstat, 1, dfr)
R[row(R) == col(R)] <- NA
R
}
Essa função “achata” a matriz criando uma matriz com quatro colunas: uma coluna para a variável i, uma para a variável j, uma para a correlação das variáveis e uma última coluna para o valor de p.
flattenSquareMatrix <- function(m) {
if ((class(m) != "matrix") | (nrow(m) != ncol(m)))
stop("Must be a square matrix.")
if (!identical(rownames(m), colnames(m)))
stop("Row and column names must be equal.")
ut <- upper.tri(m)
data.frame(i = rownames(m)[row(m)[ut]], j = rownames(m)[col(m)[ut]], cor = t(m)[ut],
p = m[ut])
}
Agora selecione uma matriz de dados, aqui o autor usou uma matriz de dados mtcars do R.
mydata <- mtcars[, c(1, 3, 4, 5, 6)]
Gerando a matriz de correlação com a função cor().
cor(mydata)
## mpg disp hp drat wt
## mpg 1.0000 -0.8476 -0.7762 0.6812 -0.8677
## disp -0.8476 1.0000 0.7909 -0.7102 0.8880
## hp -0.7762 0.7909 1.0000 -0.4488 0.6587
## drat 0.6812 -0.7102 -0.4488 1.0000 -0.7124
## wt -0.8677 0.8880 0.6587 -0.7124 1.0000
Matriz de correlação com o valor de p.
cor.prob(mydata)
## mpg disp hp drat wt
## mpg NA 9.380e-10 1.788e-07 1.776e-05 1.294e-10
## disp -0.8476 NA 7.143e-08 5.282e-06 1.222e-11
## hp -0.7762 7.909e-01 NA 9.989e-03 4.146e-05
## drat 0.6812 -7.102e-01 -4.488e-01 NA 4.784e-06
## wt -0.8677 8.880e-01 6.587e-01 -7.124e-01 NA
“Achatando” a matriz.
flattenSquareMatrix(cor.prob(mydata))
## i j cor p
## 1 mpg disp -0.8476 9.380e-10
## 2 mpg hp -0.7762 1.788e-07
## 3 disp hp 0.7909 7.143e-08
## 4 mpg drat 0.6812 1.776e-05
## 5 disp drat -0.7102 5.282e-06
## 6 hp drat -0.4488 9.989e-03
## 7 mpg wt -0.8677 1.294e-10
## 8 disp wt 0.8880 1.222e-11
## 9 hp wt 0.6587 4.146e-05
## 10 drat wt -0.7124 4.784e-06
Gerando o lindo gráfico, genial.
library(PerformanceAnalytics)
## Loading required package: zoo
## Attaching package: 'zoo'
## The following object(s) are masked from 'package:base':
##
## as.Date, as.Date.numeric
## Loading required package: xts
## Attaching package: 'PerformanceAnalytics'
## The following object(s) are masked from 'package:graphics':
##
## legend
chart.Correlation(mydata)