Curso de R no IPÊ
Script do dia 4.2: PCA
Marcos V. C. Vital e participantes da turma de 2014
Para mais material, acesse http://marcosvital.wordpress.com/
Este documento é a 2a 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 realizar uma análise de componentes principais, e vamos explorar seus resultados utlizando-os de várias maneiras.
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 ...
A tal da PCA
Existem diversas funções diferentes no R para calcular uma PCA. Aqui nós vamos usar a função rda do pacote vegan.
Não há muito mistério: basta executar a função rda com o argumento scale=T, que faz com que os dados sejam padronizados (uma característica da PCA). Depois, é só chamar o summary. Note que o resultado será um pouco longo, ok?
Aí está:
dados.amb<-dados[,3:6]
library(vegan)## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.0-10resultado.pca<-rda(dados.amb, scale=T)
summary(resultado.pca)##
## Call:
## rda(X = dados.amb, scale = T)
##
## Partitioning of correlations:
## Inertia Proportion
## Total 4 1
## Unconstrained 4 1
##
## Eigenvalues, and their contribution to the correlations
##
## Importance of components:
## PC1 PC2 PC3 PC4
## Eigenvalue 2.9765 0.8225 0.15353 0.04753
## Proportion Explained 0.7441 0.2056 0.03838 0.01188
## Cumulative Proportion 0.7441 0.9497 0.98812 1.00000
##
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## * General scaling constant of scores: 3.281818
##
##
## Species scores
##
## PC1 PC2 PC3 PC4
## Temperatura 1.5353 -0.2536 -0.516797 0.062751
## Cobertura -1.5698 0.2432 -0.348431 -0.218239
## Luz 1.5912 -0.2436 0.157733 -0.276448
## Flores 0.8129 1.4254 -0.005553 0.001159
##
##
## Site scores (weighted sums of species scores)
##
## PC1 PC2 PC3 PC4
## área1 -0.84248 -0.328739 -0.70072 0.78732
## área2 -0.43311 0.333202 0.39672 0.20386
## área3 -0.10240 -0.554527 0.23828 -1.37674
## área4 0.03636 0.130896 0.80606 -0.06254
## área5 -0.49144 -0.460239 0.53132 0.44512
## área6 -0.83250 0.301032 -0.98034 0.32923
## área7 -0.64405 0.008396 0.52433 0.54177
## área8 -0.77822 0.592733 0.08087 -0.20293
## área9 0.08170 -0.469072 0.39966 -0.10143
## área10 -0.88643 -0.871534 0.06352 -1.33527
## área11 -0.63711 0.443519 -1.28259 -0.43492
## área12 0.05798 -1.096462 -0.78626 0.03589
## área13 -0.02232 -1.170334 1.32489 0.84332
## área14 -0.52560 0.005283 0.54777 -0.35114
## área15 -0.81449 0.103509 -0.05313 0.34057
## área16 0.04581 -0.644038 -0.66399 0.43237
## área17 -0.52146 0.255348 -0.04659 -0.68786
## área18 0.91741 -0.628082 0.34388 0.24362
## área19 -0.01556 1.303795 0.37912 0.15126
## área20 0.98503 0.029162 -0.77138 -0.18492
## área21 0.81241 -0.443408 -0.40845 -0.35618
## área22 0.63011 0.488145 -0.07909 0.54670
## área23 0.12729 1.033580 0.96023 0.35001
## área24 1.07618 0.234964 0.04026 -0.37270
## área25 0.26441 -0.004436 -0.29826 1.13692
## área26 0.51127 0.426610 -0.64239 0.39092
## área27 0.29284 0.616231 -0.27960 0.36667
## área28 0.46060 0.777148 -0.03278 -1.05485
## área29 0.48418 0.256979 0.61396 -0.68838
## área30 0.76358 -0.669661 -0.22528 0.06431Plotando de maneira rápida
Não tem mistério aqui, vamos usar a função biplot:
biplot(resultado.pca)
Se quiser plotar outros componentes, você pode usar o argumento choices. Experimente, por exemplo:
biplot(resultado.pca, choices=c(1,3)) #Plotando os componentes 1 e 3
Incluindo uma variável categórica no grpafico:
Podemos delimitar os pontos de mata primária e de mata secundária no nosso gráfico assim:
biplot(resultado.pca, display="sites")
ordihull(resultado.pca, group=Ambiente,
show="Mata_Primária", col="green")
ordihull(resultado.pca, group=Ambiente,
show="Mata_Secundária", col="red")
Os detalhes:
- usamos o argumento
display="sites"para mostrar apenas os locais de coleta, sem as variáveis; - usamos a função
ordihullpara criar polígonos;
Se quiser, confira também o mesmo procedimento acima com a função ordiellipse.
Brincando com os resultados da PCA
DE novo, vamos explorar os dados sem muito compromisso, apenas para ver como podemos lidar com este tipo de análise no R.
Extraindo os valores dos componentes principais
Esta é uma etapa essencial, se quisermos usar os resultados da PCA em outras análises. É bem simples:
pc1<-resultado.pca$CA$u[,1]
pc2<-resultado.pca$CA$u[,2]Exemplinho rápido com teste T:
t.test(pc1~Ambiente, var.equal=T)##
## Two Sample t-test
##
## data: pc1 by Ambiente
## t = -7.599, df = 28, p-value = 2.811e-08
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.3838587 -0.2208520
## sample estimates:
## mean in group Mata_Primária mean in group Mata_Secundária
## -0.1310207 0.1713347boxplot(pc1~Ambiente)
E exemplinho rápido com regressão:
abundâncias<-rowSums(dados[,7:40])
modelo<-lm(abundâncias~pc1)
summary(modelo)##
## Call:
## lm(formula = abundâncias ~ pc1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -28.485 -9.075 1.509 10.140 44.424
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 60.933 2.931 20.788 < 2e-16 ***
## pc1 75.435 16.055 4.699 6.32e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 16.06 on 28 degrees of freedom
## Multiple R-squared: 0.4409, Adjusted R-squared: 0.4209
## F-statistic: 22.08 on 1 and 28 DF, p-value: 6.323e-05plot(abundâncias~pc1, pch=16)
abline(modelo, lty=2, col="red")
Aqui notem que eu usei a função rowSums para calcular as abundâncias totais das espécies.
E é isso aí, dois do dia já foram, ainda falta um. Nos vemos lá!