É comum em biologia nos deparamos com conjuntos de dados com muitas variáveis. Por exemplo, se investigamos a resposta de estresse de um organismo perante a um agente estressor, há muitas respostas metabólicas possíveis e muitas vezes analisar cada uma separadamente pode trazer problemas depois para interpretar os resultados.
Nesta aula aprenderemos um método muito utilizado neste contexto. Chama-se Análise de Componentes Principais (PCA). Vou ilustrar esta abordagem utilizando um conjunto de dados gerado em meu laboratório sobre as tolerâncias ambientais em espécies de primatas neotropicais.
Primeiramente é necessário baixar um pacote para poder ler um arquivo disponível online:
install.packages("RCurl")
## Error: trying to use CRAN without setting a mirror
library(RCurl)
## Loading required package: bitops
drop <- getURL("https://dl.dropbox.com/u/472794/dados.csv")
dat <- read.csv(textConnection(drop))
head(dat)
## alt bio1 bio2 bio3 bio4 bio5 bio6 bio7 bio8
## Alouatta_arctoidea 223.5 262.8 100.07 75.70 672.5 332.2 200.5 131.7 258.1
## Alouatta_belzebul 181.1 258.8 108.68 79.06 432.2 329.2 191.4 137.8 255.4
## Alouatta_caraya 397.8 238.0 123.72 66.10 1797.1 325.2 137.7 187.5 249.7
## Alouatta_discolor 228.6 256.8 114.65 72.83 492.4 337.8 179.0 158.8 252.2
## Alouatta_guariba 534.9 199.9 110.46 59.76 2492.7 288.5 104.1 184.4 217.0
## Alouatta_juara 387.5 251.5 97.83 82.54 497.7 312.0 193.8 118.2 249.0
## bio9 bio10 bio11 bio12 bio13 bio14 bio15 bio16 bio17
## Alouatta_arctoidea 262.6 271.2 254.3 1413 244.8 10.70 70.31 678.1 50.07
## Alouatta_belzebul 260.0 264.0 253.3 1946 335.2 23.14 69.60 930.8 91.46
## Alouatta_caraya 216.4 256.7 212.3 1430 239.5 22.45 64.79 653.8 84.01
## Alouatta_discolor 255.9 262.8 251.3 2199 353.6 24.34 64.05 985.8 100.87
## Alouatta_guariba 178.5 229.1 166.3 1423 196.1 62.06 40.33 534.2 214.67
## Alouatta_juara 251.5 256.9 244.3 2543 314.9 104.37 34.24 874.4 365.32
## bio18 bio19
## Alouatta_arctoidea 217.7 492.7
## Alouatta_belzebul 230.5 690.6
## Alouatta_caraya 424.6 112.6
## Alouatta_discolor 349.9 740.8
## Alouatta_guariba 477.2 236.8
## Alouatta_juara 516.7 684.5
O comando head é útil para visualizar o começo de um conjunto de dados grande (como esse). Você verá na primeira coluna o nome de cada espécie, seguido de 20 colunas. A primeira é a altitude média ao longo da distribuição geográfica de cada espécie, enquanto as outras variáveis são:
BIO1 = Annual Mean Temperature
BIO2 = Mean Diurnal Range (Mean of monthly (max temp - min temp))
BIO3 = Isothermality (BIO2/BIO7) (* 100)
BIO4 = Temperature Seasonality (standard deviation *100)
BIO5 = Max Temperature of Warmest Month
BIO6 = Min Temperature of Coldest Month
BIO7 = Temperature Annual Range (BIO5-BIO6)
BIO8 = Mean Temperature of Wettest Quarter
BIO9 = Mean Temperature of Driest Quarter
BIO10 = Mean Temperature of Warmest Quarter
BIO11 = Mean Temperature of Coldest Quarter
BIO12 = Annual Precipitation
BIO13 = Precipitation of Wettest Month
BIO14 = Precipitation of Driest Month
BIO15 = Precipitation Seasonality (Coefficient of Variation)
BIO16 = Precipitation of Wettest Quarter
BIO17 = Precipitation of Driest Quarter
BIO18 = Precipitation of Warmest Quarter
BIO19 = Precipitation of Coldest Quarter
Todos estes dados foram obtidos a partir de um banco de dados mundial disponível aqui (http://www.worldclim.org/) e descritos nesse artigo:
Hijmans, R.J., S.E. Cameron, J.L. Parra, P.G. Jones and A. Jarvis, 2005. Very high resolution interpolated climate surfaces for global land areas. International Journal of Climatology 25: 1965-1978.
Note que diferentes variáveis tem escalas distintas. Variáveis que tenham escalas maiores em termos absolutos podem influenciar disproporcionalmente mais as análises. Para evitar isso, vamos escalonar as variáveis. Fazemos isso utilizando os seguintes passos: (1) calculamos a média e o desvio-padrão de cada variável (coluna); (2) subtraímos de cada valor a média da coluna; (3) dividimos de cada valor o desvio-padrão da coluna.
dat <- scale(dat)
Agora estamos prontos para realizar a PCA, baseada na matriz de covariância dos dados. A PCA funciona da seguinte forma. Tomando os dados originais, avalia-se em qual direção há mais variação no espaço multivariado. Este eixo é chamado de eixo principal 1 (PC1). Busca-se então um segundo eixo (PC2), onde exista mais variação, desde que seja ortogonal ao PC1, e assim por diante. Vejamos nos nossos dados o que podemos observar.
PCA <- prcomp(dat, cor = F)
summary(PCA)
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 3.197 2.211 1.3454 1.2653 0.8879 0.529 0.40201
## Proportion of Variance 0.511 0.244 0.0905 0.0801 0.0394 0.014 0.00808
## Cumulative Proportion 0.511 0.755 0.8459 0.9260 0.9654 0.979 0.98750
## PC8 PC9 PC10 PC11 PC12 PC13
## Standard deviation 0.3066 0.24470 0.22298 0.14794 0.10778 0.0631
## Proportion of Variance 0.0047 0.00299 0.00249 0.00109 0.00058 0.0002
## Cumulative Proportion 0.9922 0.99520 0.99768 0.99878 0.99936 0.9996
## PC14 PC15 PC16 PC17 PC18 PC19
## Standard deviation 0.05791 0.0454 0.04283 0.03285 0.02167 0.00787
## Proportion of Variance 0.00017 0.0001 0.00009 0.00005 0.00002 0.00000
## Cumulative Proportion 0.99972 0.9998 0.99992 0.99997 1.00000 1.00000
## PC20
## Standard deviation 4.08e-05
## Proportion of Variance 0.00e+00
## Cumulative Proportion 1.00e+00
Veja a segunda linha. Ela indica que, por mais que estejamos tratando de 20 variáveis, as mudanças ao longo do PC1 explicam pouco mais do que metade de toda a variância nos dados (51.11%). O PC2 explica 24.4% da variância e assim por diante. Normalmente interpreta-se os PCs com standard deviation acima de 1, indicando que o PC explica mais do que qualquer variável isoladamente, mas não há uma regra rígida sobre isso. Uma vez definidos os eixos que serão interpretados, precisamos agora compreender o que eles querem dizer. O primeiro passo é compreender como as variáveis originais se relacionam aos PCs, os chamados loadings. No código a seguir, podemos ver os loadings dos primeiros 4 eixos. Para facilitar, vamos arredondar para duas casas decimais
round(PCA$rotation[, 1:4], digits = 2)
## PC1 PC2 PC3 PC4
## alt 0.19 0.26 0.14 0.35
## bio1 -0.27 -0.22 0.01 -0.06
## bio2 0.17 -0.13 -0.32 0.32
## bio3 -0.23 0.17 0.23 0.25
## bio4 0.23 -0.04 -0.24 -0.36
## bio5 -0.17 -0.35 -0.18 -0.01
## bio6 -0.30 -0.08 0.19 -0.03
## bio7 0.24 -0.17 -0.37 0.03
## bio8 -0.22 -0.27 -0.14 -0.25
## bio9 -0.29 -0.16 0.12 0.06
## bio10 -0.23 -0.28 -0.08 -0.21
## bio11 -0.28 -0.16 0.09 0.07
## bio12 -0.25 0.21 -0.25 0.09
## bio13 -0.23 0.08 -0.33 0.31
## bio14 -0.18 0.33 -0.09 -0.22
## bio15 0.11 -0.32 -0.14 0.37
## bio16 -0.24 0.08 -0.33 0.28
## bio17 -0.18 0.33 -0.10 -0.20
## bio18 -0.05 0.30 -0.45 -0.20
## bio19 -0.25 0.14 0.01 0.16
Quais variáveis apresentam maiores valores (posivos ou negativos) em cada eixo? O que isso quer dizer biologicamente?
É possível agora avaliar a posição de cada espécie neste novo espaço multivariado. Como nosso conjunto de dados é grande demais, a visualização fica comprometida, mas é possível visualizar as espécies nos extremos. Neste caso, olharemos apenas os PCs 1 e 2.
biplot(PCA)
Procure no google pelo nome das espécies em posições mais extremas do biplot. O local onde elas ocorrem fazem sentido com a sua interpretação sobre as variáveis climáticas? Para ter todos os scores de cada espécie, você pode usar o seguinte comando:
PCA$x[, 1:4]
## PC1 PC2 PC3 PC4
## Alouatta_arctoidea -0.012163 -2.89214 1.517588 0.3124266
## Alouatta_belzebul -0.684546 -2.03377 0.363862 1.6015252
## Alouatta_caraya 4.104055 -1.78026 -1.653711 -0.2202022
## Alouatta_discolor -0.303649 -1.85462 -0.967789 1.5878130
## Alouatta_guariba 6.984879 2.07520 -1.114884 -1.5188075
## Alouatta_juara -1.394705 1.60545 0.546209 -0.2032221
## Alouatta_macconnelli -2.005652 -0.32077 0.798883 0.5543846
## Alouatta_nigerrima -1.988746 -1.42069 0.524614 -0.0090665
## Alouatta_palliata -0.213621 0.48931 -0.533092 0.4867187
## Alouatta_pigra 1.652336 -1.96327 -1.271996 -1.2358434
## Alouatta_puruensis -0.064673 -1.23091 -0.377935 0.5324445
## Alouatta_sara 1.705498 -1.18161 -1.776121 -0.6929199
## Alouatta_seniculus -1.293117 1.73918 0.574015 -0.1600517
## Alouatta_ululata -0.333035 -4.15632 1.016839 1.7191600
## Aotus_azarai 1.860016 -2.15567 -1.106777 0.7021896
## Aotus_brumbacki -1.921751 -1.39344 0.287981 0.9090895
## Aotus_griseimembra -0.718813 -0.60075 1.002886 0.6320735
## Aotus_jorgehernandezi -8.209435 10.16635 -6.298630 0.8137207
## Aotus_lemurinus 8.495049 9.90347 3.548576 4.3381461
## Aotus_miconax 10.048898 6.57605 2.557232 4.1259451
## Aotus_nancymaae -1.002236 1.13371 0.532645 -0.8999271
## Aotus_nigriceps -0.216116 -1.07681 -0.294481 0.3547316
## Aotus_trivirgatus -2.614748 -0.11118 0.383245 0.7795588
## Aotus_vociferans -3.050716 2.29087 0.242127 -1.1228151
## Aotus_zonalis -3.361224 3.54433 -0.942331 0.8756663
## Ateles_belzebuth -2.328102 1.55280 0.115104 0.0008087
## Ateles_chamek 0.119852 -1.10588 -0.593671 0.1435482
## Ateles_fusciceps -2.268132 3.63560 -0.657584 0.8080913
## Ateles_geoffroyi 1.877337 -0.85911 -1.167185 0.1072401
## Ateles_hybridus -0.360200 0.30722 0.745505 0.8927051
## Ateles_marginatus 0.221748 -1.98306 -1.337935 1.7692916
## Ateles_paniscus -2.136902 -0.72805 1.084295 0.2576182
## Brachyteles_arachnoides 7.382206 2.84855 -1.822980 -1.1814906
## Brachyteles_hypoxanthus 5.906895 -0.15300 -1.184201 -0.1277527
## Cacajao_ayresi -3.523090 0.24952 0.899248 -0.2524900
## Cacajao_calvus -1.446157 0.37612 0.195534 -1.0749625
## Cacajao_hosomi -3.437430 2.30821 -0.207764 -0.1920172
## Cacajao_melanocephalus -3.644465 1.37256 0.038385 -0.6032477
## Callibella_humilis -2.649333 -0.46831 -0.037850 -0.1358816
## Callicebus_aureipalatii 0.754226 -0.48975 -1.287874 -0.6286253
## Callicebus_baptista -3.238748 -1.28296 1.555818 -0.5322219
## Callicebus_barbarabrownae 5.188133 -0.57554 2.266693 -0.9276390
## Callicebus_bernhardi 0.045154 -1.71900 -0.811232 1.1283258
## Callicebus_brunneus 0.783688 -1.73132 -0.867755 0.4618041
## Callicebus_caligatus -2.981172 -0.48526 1.173109 -0.4990416
## Callicebus_cinerascens 0.386225 -1.89450 -0.940682 1.0475789
## Callicebus_coimbrai 2.393742 -0.96801 2.417967 -1.8310971
## Callicebus_cupreus -1.268753 0.20409 0.100648 -0.5392483
## Callicebus_discolor -2.391829 3.01986 0.288127 -1.5297388
## Callicebus_donacophilus 1.941942 -1.91872 -2.019363 -0.5171093
## Callicebus_dubius 0.006492 -1.71171 -0.086031 0.0883626
## Callicebus_hoffmannsi -2.487517 -1.32116 0.923776 0.0730421
## Callicebus_lucifer -3.332027 2.68085 0.133621 -1.6734300
## Callicebus_lugens -2.559451 0.61200 0.151724 0.4117408
## Callicebus_medemi -3.222312 2.96936 -0.356494 -0.6416672
## Callicebus_melanochir 3.257753 0.64373 2.661814 -2.5903764
## Callicebus_modestus 0.890361 -2.01828 -1.046596 -0.8736060
## Callicebus_moloch 0.303032 -2.04034 -1.074493 1.8221243
## Callicebus_nigrifrons 6.547612 0.96738 -1.917540 0.6400270
## Callicebus_oenanthe 3.431152 1.66584 2.471983 0.0827987
## Callicebus_olallae 0.874152 -2.22646 -1.011762 -0.9042244
## Callicebus_ornatus -2.262147 -0.65392 0.094882 0.8808604
## Callicebus_pallescens 4.605308 -3.78536 -1.804178 -2.6558353
## Callicebus_personatus 5.172583 -0.65452 -0.458786 -0.3695050
## Callicebus_purinus -2.535763 0.25697 0.609336 -0.6031178
## Callicebus_regulus -1.804632 1.16167 0.267112 -0.9190372
## Callicebus_stephennashi -1.423828 -1.45336 0.275685 0.1483968
## Callicebus_torquatus -3.761064 1.80285 0.616858 -1.1541968
## Callimico_goeldii -1.010886 0.76269 -0.133035 -0.8394084
## Callithrix_aurita 7.229806 1.85485 -1.683587 0.0491308
## Callithrix_flaviceps 5.954317 -0.02611 -1.048193 0.1224137
## Callithrix_geoffroyi 4.260594 -0.75150 0.427617 -1.0621329
## Callithrix_jacchus 1.983116 -3.51072 1.301097 0.7669140
## Callithrix_kuhlii 3.976948 0.42791 2.504022 -2.1161045
## Callithrix_penicillata 4.130233 -1.48579 -0.851905 1.2366751
## Cebuella_pygmaea -1.647575 0.70736 0.110598 -0.7811786
## Cebus_albifrons -1.418174 0.44769 0.074086 -0.2075483
## Cebus_apella -0.691561 -1.48232 -0.015165 0.8099213
## Cebus_capucinus -1.183091 1.66814 -0.481823 0.5071922
## Cebus_cay 4.719619 -1.62009 -1.886014 -2.2631898
## Cebus_flavius 1.376397 -0.38454 2.046576 -1.9103067
## Cebus_kaapori -2.455563 -2.12212 0.939006 1.8385009
## Cebus_libidinosus 2.720324 -2.69842 -0.471543 1.1082170
## Cebus_macrocephalus -1.758221 0.63805 0.150194 -0.4816511
## Cebus_nigritus 7.370040 2.16836 -1.838923 -1.0285354
## Cebus_olivaceus -1.440468 -0.97327 0.909406 0.6645169
## Cebus_robustus 4.342905 -0.72628 0.399263 -0.9836157
## Cebus_xanthosternos 4.802540 -1.03052 1.556051 -0.4327836
## Chiropotes_albinasus 0.193731 -1.83180 -1.008593 1.3413888
## Chiropotes_chiropotes -2.091513 -0.43119 0.724331 0.7164827
## Chiropotes_satanas -1.879060 -2.47453 0.851467 1.7558455
## Chiropotes_utahickae -0.435381 -1.95213 -0.031576 1.7148703
## Lagothrix_cana 0.490617 -0.28812 -0.203185 0.7344889
## Lagothrix_lagotricha -3.305608 1.93866 -0.131049 -0.8519223
## Lagothrix_lugens 0.534238 2.79332 1.168306 1.8294002
## Lagothrix_poeppigii -1.710003 1.49674 0.320210 -1.1724628
## Leontopithecus_caissara 3.071836 2.00650 -3.817857 -3.9152967
## Leontopithecus_chrysomelas 2.582780 0.63722 2.817585 -3.0652634
## Leontopithecus_chrysopygus 7.325879 1.18921 -1.772373 -1.0762183
## Leontopithecus_rosalia 3.942450 -0.30553 1.529163 -3.3652676
## Mico_acariensis -1.161034 -1.64279 -0.318144 0.2125959
## Mico_argentatus -1.835353 -1.40443 1.275199 1.0697040
## Mico_chrysoleucus -2.424987 -0.97219 1.068217 -0.6338693
## Mico_emiliae 1.661135 -2.43779 -3.038964 2.5015896
## Mico_humeralifer -2.542359 -1.27747 0.897673 0.1040111
## Mico_intermedius 1.762988 -2.10772 -1.489818 1.8238556
## Mico_leucippe -1.922239 -0.95566 0.226445 0.3465376
## Mico_manicorensis -2.394714 -0.58053 -0.361133 0.1180926
## Mico_marcai -0.951640 -1.73491 -0.891665 0.9256567
## Mico_mauesi -2.060264 -1.34092 0.458454 -0.0466520
## Mico_melanurus 3.107643 -2.37878 -1.645102 -0.1666119
## Mico_nigriceps -0.885744 -1.90057 -0.447024 0.6437885
## Mico_rondoni 0.297572 -2.21226 -0.552034 0.7902344
## Mico_saterei -2.365838 -1.08054 0.816539 -0.5016042
## Oreonax_flavicauda 7.545197 5.22496 3.015596 2.2401590
## Pithecia_aequatorialis -2.843565 3.28123 -0.033807 -1.5293959
## Pithecia_albicans -2.670853 -0.04207 1.224190 -0.6847136
## Pithecia_irrorata -0.320708 -1.28410 -0.312569 0.3593995
## Pithecia_monachus -2.168877 1.37701 0.263758 -1.2676517
## Pithecia_pithecia -2.318692 -0.82814 1.127161 0.1002385
## Saguinus_bicolor -3.317712 -0.69109 1.682149 -0.9320489
## Saguinus_fuscicollis -1.005353 0.51837 0.008022 -0.5247549
## Saguinus_geoffroyi -5.213306 3.74385 -1.818336 0.4896348
## Saguinus_imperator 0.755817 -0.62558 -0.639814 -0.0825823
## Saguinus_inustus -3.671739 1.92962 -0.150678 -0.8471067
## Saguinus_labiatus -1.420480 -0.85275 0.299008 -0.2851307
## Saguinus_leucopus -0.660095 3.08872 0.180883 2.0434701
## Saguinus_martinsi -3.452137 -1.46325 1.067042 0.0814697
## Saguinus_melanoleucus -0.072319 -0.42355 0.041955 -0.3692087
## Saguinus_midas -2.273084 -0.91905 0.996897 0.3517931
## Saguinus_mystax -1.812167 0.47885 0.228659 -0.7503766
## Saguinus_niger -1.475994 -2.12623 0.741859 1.6375595
## Saguinus_nigricollis -2.605019 2.48408 0.255124 -0.9269966
## Saguinus_oedipus -1.861085 -1.89874 1.404704 0.2875480
## Saguinus_tripartitus -2.699233 3.10827 0.182971 -1.2237926
## Saimiri_boliviensis 1.530966 0.31141 -0.517834 0.4061818
## Saimiri_oerstedii -2.434640 0.46152 -3.039621 2.6793873
## Saimiri_sciureus -1.681150 -0.22059 0.202383 0.4223447
## Saimiri_ustus -0.010033 -1.59866 -0.376266 0.7552005
## Saimiri_vanzolinii -4.598494 2.59138 -0.440462 -1.3305725