Curso de R no IPÊ
Script do dia 4.1: Matrizes e Cluster
Marcos V. C. Vital e participantes da turma de 2014
Para mais material, acesse http://marcosvital.wordpress.com/
Este documento é a 1a parte do quarto de uma série que nasceu durante as atividades do “Curso de introdução ao uso do software R e suas aplicações estatísticas nas ciências biológicas”, realizado no IPÊ em dezembro de 2014.
Nesta parte do guia nós vamos criar matrizes de distância e representá-las em dendogramas hierárquicos.
Como sempre, todo o documento foi criado utilizando o software R, o programa RStudio e a linguagem Markdown. Para saber mais, acesse:
Antes de começar, carregue os dados necessários:
setwd("C:/R/Curso IPÊ")
dados<-read.table("insetos.txt", h=T, row.names=1)
attach(dados)
#Aqui um summary vai ser muito carregado, então vou usar um str:
str(dados)## 'data.frame': 30 obs. of 40 variables:
## $ Ambiente : Factor w/ 2 levels "Mata_Primária",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Gado : Factor w/ 2 levels "Ausente","Presente": 2 1 2 2 2 1 1 1 2 2 ...
## $ Temperatura: num 21.8 21.3 23.8 23.2 21.5 21.7 20.3 19.7 24.7 20.2 ...
## $ Cobertura : int 73 62 57 49 59 78 64 74 48 75 ...
## $ Luz : num 53.5 61 71.5 69.5 61.5 53 58 56 71 60.5 ...
## $ Flores : int 14 28 17 29 15 24 21 29 20 5 ...
## $ sp.1 : int 4 5 1 3 5 1 1 0 6 3 ...
## $ sp.2 : int 0 0 0 0 0 0 3 1 4 0 ...
## $ sp.3 : int 7 3 2 6 7 2 2 1 8 3 ...
## $ sp.4 : int 2 7 0 4 2 4 5 3 3 1 ...
## $ sp.5 : int 0 2 2 1 1 0 3 4 0 0 ...
## $ sp.6 : int 0 5 2 14 4 3 3 9 8 0 ...
## $ sp.7 : int 3 2 3 2 1 4 3 2 1 3 ...
## $ sp.8 : int 2 4 3 5 2 3 0 0 1 0 ...
## $ sp.9 : int 1 2 0 0 0 0 2 0 1 0 ...
## $ sp.10 : int 1 3 5 2 1 1 0 0 1 1 ...
## $ sp.11 : int 0 0 0 4 0 2 0 2 0 0 ...
## $ sp.12 : int 1 5 0 2 0 2 4 5 0 0 ...
## $ sp.13 : int 3 6 0 3 3 2 3 1 4 2 ...
## $ sp.14 : int 3 4 2 10 5 2 2 5 8 0 ...
## $ sp.15 : int 0 0 0 0 1 0 0 0 1 0 ...
## $ sp.16 : int 3 3 0 2 0 0 3 0 1 0 ...
## $ sp.17 : int 0 0 0 1 0 0 0 0 0 0 ...
## $ sp.18 : int 0 1 2 3 1 3 0 0 2 0 ...
## $ sp.19 : int 0 0 2 2 3 0 0 2 0 1 ...
## $ sp.20 : int 0 3 0 0 0 4 0 2 0 0 ...
## $ sp.21 : int 0 1 3 4 2 1 3 3 3 2 ...
## $ sp.22 : int 1 2 1 2 2 2 1 2 2 2 ...
## $ sp.23 : int 1 0 1 2 1 1 2 1 2 1 ...
## $ sp.24 : int 1 1 0 2 1 2 0 1 1 1 ...
## $ sp.25 : int 0 1 0 2 0 1 1 0 0 1 ...
## $ sp.26 : int 2 2 1 0 1 1 2 1 1 1 ...
## $ sp.27 : int 0 0 0 0 0 0 3 0 0 0 ...
## $ sp.28 : int 1 0 0 0 0 0 0 0 0 2 ...
## $ sp.29 : int 0 0 0 0 1 0 0 0 0 0 ...
## $ sp.30 : int 0 0 0 0 0 0 0 0 0 0 ...
## $ sp.31 : int 0 5 0 0 0 0 1 0 0 0 ...
## $ sp.32 : int 1 1 1 0 1 1 1 0 1 1 ...
## $ sp.33 : int 0 0 0 0 1 0 0 0 1 0 ...
## $ sp.34 : int 0 0 0 0 1 0 0 0 0 0 ...Notou uma coisa diferente?
Nós usamos o argumento row.names para informar o R os nomes das linhas!
Passo a passo para criar um cluster:
Etapa 0: se necessário, padronize os seus dados
dados.amb<-dados[,3:6] #Separei os dados de interesse
dados.amb.pad<-scale(dados.amb) #Padronizei os dadosA novidade aqui é a função scale que padroniza os dados (subtrai da média e divide pelo desvio padrão).
Etapa 1: crie sua matriz de distância
amb.dist<-dist(dados.amb.pad, method="euclid")Aqui usamos a função dist, e escolhemos a distância euclidiana com o argumento method
Etapa 2: crie o seu cluster
amb.cluster<-hclust(amb.dist, method="average")Usamos a função hclust, e definimos o método de UPGMA pelo argumento method
Etapa 3: calcule o coeficiente de correlação cofenética
matriz.cof<-cophenetic(amb.cluster) #Cria a matriz coefenética
cor(amb.dist, matriz.cof) #Calcula o coeficiente## [1] 0.7388642Aqui usamos a função cophenetic para extrair a matriz cofenética, e a função cor para calcular a correlação de pearson.
O coeficiente apresentou um valor adequado, então vamos lá plotar o nosso cluster!
Plotando o nosso cluster
Existem diversas maneiras de se plotar um cluster no R. Vamos explorar apenas algumas, ok?
Rápidão!
O mais rápido e simples é:
plot(amb.cluster)
Não é la muito bonito, mas quebra galho para uma visualização rápida.
Você pode usar o argumento hang para que ele estique os ramos e deixe os nomes lado a lado, assim:
plot(amb.cluster, hang=-1)
De lado
Você pode colocar o cluster na horizontal convertendo ele em outra categoria de objeto. São duas alternativas: você pode usar a função as.phylo do pacote ape; ou pode usar a as.dendrogram.
Usando a primeira:
library(ape)
plot(as.phylo(amb.cluster))
E a segunda (neste caso precisamos do argumento horiz para deixá-lo na horizontal:
plot(as.dendrogram(amb.cluster), horiz=T) #horiz=T faz na horizontal
Mas nos dois casos uma coisa não ficou legal: os nomes ficaram espremidos…
Isto pode ser resolvido de uma maneira simples: aumentando o tamanho da janela do gráfico. Vamos lá:
windows(width=8, height=8, rescale="fixed")
plot(as.dendrogram(amb.cluster), horiz=T)
Aqui temos novidades:
- a função
windowsabre uma janela para gráfico com o tamanho que determinarmos; - os argumentos
widtheheighdeterminam a largura e altura da janela, respectivamente; - o argumento
rescale="fixed"faz com que a figura não mude de tamanho quando mexemos na janela;
Alterando nomes
Uma das coisas legais do R é que podemos mexer aqui e ali e encontrar maneiras de solucionar problemas, representar coisas e tudo mais sem ter que voltar atrás, mexer nos dados originais e coisas do tipo.
Imagine, por exemplo, que você queira alterar os nomes das unidades amostrais para identificá-las como sendo pertencentes aos dois tipos de mata descritos na variável ambiente. Como fazer isso sem sair do R?
Vamos lá:
códigos<-ifelse(dados$Ambiente=="Mata_Primária", "MP", "MS")
novos.nomes<-paste(row.names(dados), códigos, sep="-")
amb.cluster$labels<-novos.nomes
plot(as.phylo(amb.cluster))
O que fizemos aqui?
ifelsevocê já conhece, e foi usada para criar códigos para serem acrescentados aos nomes;pastepermite que eu junte nomes, expressões e coisas assim, e foi usada para criarmos os nomes novos;$labelsfoi usado para acessarmos os nomes usados no cluster.
Bacana, né?
Mas peraí, como você sabia desse esquema do $labels?
Existe uma função bem básica mas extremamente útil no R, chamada ls.
Esta função serve, inicialmente, para listarmos os objetos na memória do R. Por exemplo:
ls()## [1] "amb.cluster" "amb.dist" "códigos" "dados"
## [5] "dados.amb" "dados.amb.pad" "matriz.cof" "novos.nomes"Mas ela também pode ser usada para listarmos o que está dentro de um objeto, o que pode nos ajudar a entender os seus componentes. No caso do cluster, eu soube que os nomes deveriam estar ali, por olhei para o resultado do ls, assim:
ls(amb.cluster)## [1] "call" "dist.method" "height" "labels" "merge"
## [6] "method" "order"
Criando um mapa de calor
Este pedacinho só vai funcionar se você tiber o arquivo ‘Coldiss.R’ que acompanha o livro Numerical Ecology with R.
Um mapa de calor é uma representação visual da matriz de distância, na qual os números são substituídos por cores em uma escala. O resultado pode ser bastante intuitivo:
source("Coldiss.R")
coldiss(amb.dist)## Loading required package: gclus
## Loading required package: cluster
Aqui nós usamos a função source para acessar o arquivo e criar a função coldiss
Escolhendo outras métricas de distância
Ao fazer o nosso cluster la em cima, trabalhamos com a distância euclidiana. Mas, como vimos, outras métricas podem ser necessárias, dependendo de coisas como a natureza dos meus dados e as minhas perguntas.
Vamos, então, olhar rapidamente para alguns exemplos bem simples, ok? Você irá precisar do pacote vegan.
Coeficiente de Bray-Curtis:
library(vegan)## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.0-10dist.bray<-vegdist(dados[,7:40], method="bray")
cluster.bray<-hclust(dist.bray, method="average")
plot(cluster.bray, hang=-1)
Aqui usamos a função vegdist, que tem um monte de métricas acessíveis pelo argumento method;
Coeficiente de Jaccard:
dist.jac<-vegdist(dados[,7:40], method="jac", binary=T)
cluster.jac<-hclust(dist.jac, method="average")
plot(cluster.jac, hang=-1)
Aqui adicionamos o argumento binary para indicar que queremos usar apenas a presença e ausência das espécies.
Outras métricas?
O pacote vegan possui uma função bem bacana, que permite criar uma métrica desejada. Vale à pena dar uma espiada nela:
?designdist
Ufa, terminamos mais um. Ainda temos mais coisas do mesmo dia, então vamos continuar!