A couple of days I received the following inquiry about phylogenetic PCA as implemented
in the phytools function phyl.pca:
“Also, another less pressing question I had was what to do about multiple samples from one species. I think the phylo.pca doesn't like 0 branch lengths (I think because of the inverted matrix error), so I've been putting in "placeholder” small branch lengths and/or averaging individuals of the same species. I'm not sure if this was more acceptable or if there was a more standard practice, though.“
Firstly, one should virtually never include interspecific data by adding terminal branches of zero (or very short) length to the tree. I won't get into why here - just don't do it please.
It is possible, however, to compute PC scores for one set of data based on a rotation computed for another set. With respect to the phylogenetic PCA, what I would thus recommend it to compute the rotation using the species means data & then if one wants to obtain scores for individuals, use this rotation to compute those scores.
I previously documented how to do this
here,
but today I also added a new S3 method called scores( ) that automates this
calculation for a "phyl.pca" option. I just
pushed
this update to GitHub so it can be installed using devtools in the typical way.
Here's a quick demo of how it works:
library(phytools)
## latest phytools version from GitHub
packageVersion("phytools")
## [1] '0.7.3'
tree
##
## Phylogenetic tree with 26 tips and 25 internal nodes.
##
## Tip labels:
## Y, K, R, O, J, W, ...
##
## Rooted; includes branch lengths.
## individual data
Xi
## [,1] [,2] [,3] [,4] [,5]
## A1 1.20787154 0.47914789 0.8623047 1.2495477 -1.42947991
## A2 2.12320487 0.47402360 0.9550751 1.6972576 -1.91713104
## B1 -3.36948268 3.77947415 -0.8104234 5.5267960 -4.22373306
## B2 -3.89240291 3.54207211 -1.3588802 5.1915551 -3.65014197
## B3 -3.63879657 3.37713991 -1.9685069 4.3689452 -4.14255235
## B4 -4.17859299 3.73085941 -1.3349099 5.3970606 -3.86395132
## B5 -2.77632771 3.40620475 -2.0374127 5.4595384 -3.94879047
## C1 2.66252966 -1.69501406 5.0620859 1.4238734 -1.00437382
## C2 2.27931674 -2.71852958 4.8953938 1.0881464 -1.39655915
## C3 3.19271894 -1.67392645 4.3500289 2.0463243 -1.44478241
## C4 1.93553298 -2.86843955 4.8429547 1.6322762 -1.23128013
## C5 2.41170225 -2.05714667 4.7793417 2.0967727 -0.59099047
## D1 -3.06618045 -0.77675201 0.7938031 1.5116715 0.94196412
## D2 -3.36196611 -0.80817106 0.7938168 0.9432299 1.13672644
## E1 0.59350436 1.38165922 -0.3316782 1.9536058 -2.77096246
## E2 0.27218954 0.56687414 0.3876198 1.8826100 -1.69035681
## E3 -0.52332659 0.59848482 0.2504620 2.5939968 -1.65883260
## F1 0.22273476 -2.45198070 2.7746535 0.3861395 1.05274671
## G1 0.29489467 -2.73762657 2.6821865 -0.8493884 1.69457218
## G2 -0.65288642 -3.54374917 2.1397764 -1.0217453 0.94721876
## G3 -0.97408415 -3.32511230 2.7265237 -1.1669954 1.91906697
## G4 -0.02536859 -2.11640698 2.8554528 -1.2779744 1.76777639
## H1 -4.01209732 -0.23426385 -2.0376439 -0.3891521 1.00892916
## H2 -3.15960646 0.08868799 -2.0180995 -0.7439570 0.79702482
## H3 -3.70773942 0.14920207 -2.3082605 -1.0694013 -0.10941989
## H4 -3.46345662 0.22065346 -1.5988715 -1.0741026 0.11413459
## H5 -3.04286821 -0.24156083 -1.2535244 -1.1774794 1.38878052
## I1 0.99431344 -0.39255396 2.7078109 3.0684490 -2.19073296
## J1 1.48012052 -3.89868107 5.3744037 0.4072786 -1.20200149
## J2 0.48200711 -3.35309585 5.2460605 1.0674284 -1.45483287
## J3 1.83140764 -2.74096320 5.3074037 0.3876488 -1.27536042
## J4 1.11622382 -4.10303090 5.5272531 1.2245890 -1.04267349
## K1 -0.59879568 -2.64773568 3.8750676 0.8384714 0.27605077
## K2 0.21024899 -2.11642761 3.6165256 2.7382009 0.69900866
## K3 -1.09273604 -2.74959089 4.3711782 0.8696402 0.64300187
## K4 -0.04152818 -2.61839939 3.7803252 1.6808116 0.29258721
## K5 -0.72210266 -2.33429710 4.1279836 2.4090130 0.78003183
## L1 0.03509973 -2.05695368 1.9254794 -0.2567052 0.68410460
## L2 -1.05658416 -0.88037564 2.6933563 0.5511140 0.86919585
## L3 -0.51878917 -1.93693089 1.6358342 0.3728830 2.03744417
## L4 -0.67541145 -1.70725799 1.2117356 0.2004287 0.01558763
## L5 -1.02782819 -1.63031446 1.0975700 -0.4001964 1.63772015
## M1 2.20607975 -1.60221108 1.4350309 -1.0468959 0.24758727
## M2 2.62837719 -1.04389879 1.3517971 0.5392242 1.07560283
## N1 -2.93057050 -2.17509303 2.6040320 4.7347668 0.15589588
## N2 -3.15345861 -2.59937773 3.3869134 4.2462922 1.39983864
## N3 -3.60232496 -2.71874020 2.8319589 4.4270392 0.24840489
## N4 -4.25563290 -2.64746706 2.8833009 3.6301153 -0.27392402
## O1 -5.28683088 -3.70228410 2.3844570 1.3629796 2.39517794
## O2 -5.38104846 -4.87728898 2.8560162 2.0700145 1.50657055
## O3 -5.48787982 -4.25493642 2.8122042 1.3717118 2.28957850
## P1 3.09231576 -5.14787727 5.6295546 -0.8423818 1.52657614
## Q1 -3.82118671 -3.11733681 4.6756523 5.0033811 0.87620972
## Q2 -4.23022707 -4.21706371 3.7660017 5.0281017 0.93682444
## R1 -2.36400839 -2.27834062 1.9794830 1.4314134 0.60331716
## R2 -2.46014028 -1.08918230 1.6520861 0.9182345 -0.96725226
## R3 -2.73368752 -1.43893280 2.0945741 2.2804404 -0.69332985
## R4 -1.08796774 -2.09122908 1.1532129 1.8427559 0.13805477
## R5 -1.89238447 -2.26101902 1.9490731 1.3904670 -0.78965765
## S1 -3.12227407 -1.54462933 1.4331181 -0.3039058 1.86762871
## S2 -3.54825239 -1.57142641 0.2298638 1.6937560 1.14127352
## S3 -2.68354538 -1.00599847 0.5909795 0.9248922 1.83179640
## S4 -2.87604229 -1.72251468 0.3321336 1.0708315 1.71451830
## T1 -2.57118182 -1.90281553 1.4933224 3.0562073 -2.09442201
## U1 0.66466684 1.65714974 1.0003654 3.9438689 -0.86037147
## U2 0.45740087 2.04635918 1.3612159 4.2285696 -0.80167998
## U3 1.14681687 1.59757388 0.8111320 3.4097557 -1.34731730
## U4 0.52640220 1.64986943 1.7871073 4.3412065 -2.19397881
## V1 0.75658437 -1.87706996 4.4759921 1.5414931 -1.62431717
## W1 2.68267252 -3.19689644 7.4425853 3.7556724 -1.68887877
## W2 2.76196919 -3.74195183 6.1999013 2.9743391 -1.78427181
## X1 0.91234712 -0.72483560 3.6801507 4.4381717 -2.46478786
## X2 -0.05238740 -0.80465206 3.2210595 2.7922503 -3.38940095
## X3 0.45835800 -0.53257849 3.5476527 3.3640370 -2.73478436
## X4 1.49751687 -0.73809794 3.9058654 4.1334758 -3.66075267
## X5 0.28602440 -0.49840370 2.7803620 4.1733508 -2.10614716
## Y1 1.25534713 -5.10787364 4.3453981 -2.7789615 2.23696428
## Y2 3.18711033 -4.93714442 4.0112931 -2.4233436 2.50171151
## Y3 1.57399176 -4.46222484 4.3241397 -2.0585551 2.36389952
## Z1 2.57286210 -2.39561802 5.4534820 2.4485488 -1.16858447
## Z2 3.40395933 -2.58146747 5.9747592 0.9585838 -1.53438227
## Z3 2.55456063 -1.98275084 6.0612324 1.5865925 -2.48843615
## Z4 2.87982996 -0.97677109 5.7028841 1.2531958 -1.32016016
## Z5 2.60403104 -2.01564983 5.3533706 1.4698986 -1.86993584
## species means
Xm<-matrix(NA,Ntip(tree),ncol(Xi),
dimnames=list(tree$tip.label,
colnames(Xi)))
for(i in 1:Ntip(tree)){
ii<-grep(tree$tip.label[i],rownames(Xi))
Xm[i,]<-colMeans(Xi[ii,,drop=FALSE])
}
Xm
## [,1] [,2] [,3] [,4] [,5]
## Y 2.0054831 -4.835747633 4.2269436 -2.42028676 2.3675251
## K -0.4489827 -2.493290133 3.9542160 1.70722742 0.5381361
## R -2.1076377 -1.831740764 1.7656858 1.57266223 -0.3417736
## O -5.3852531 -4.278169835 2.6842258 1.60156865 2.0637757
## J 1.2274398 -3.523942753 5.3637802 0.77173620 -1.2437171
## W 2.7223209 -3.469424136 6.8212433 3.36500574 -1.7365753
## P 3.0923158 -5.147877265 5.6295546 -0.84238176 1.5265761
## N -3.4854967 -2.535169504 2.9265513 4.25955340 0.3825538
## Q -4.0257069 -3.667200260 4.2208270 5.01574141 0.9065171
## G -0.3393611 -2.930723754 2.6009849 -1.07902587 1.5821586
## L -0.6487026 -1.642366534 1.7127951 0.09350482 1.0488105
## X 0.6203718 -0.659713558 3.4270181 3.78025711 -2.8711746
## A 1.6655382 0.476585745 0.9086899 1.47340267 -1.6733055
## S -3.0575285 -1.461142224 0.6465237 0.84639345 1.6388042
## D -3.2140733 -0.792461533 0.7938099 1.22745074 1.0393453
## H -3.4771536 -0.003456233 -1.8432800 -0.89081847 0.6398898
## E 0.1141224 0.849006057 0.1021345 2.14340420 -2.0400506
## U 0.6988217 1.737738058 1.2399551 3.98085017 -1.3008369
## I 0.9943134 -0.392553964 2.7078109 3.06844904 -2.1907330
## B -3.5711206 3.567150066 -1.5020266 5.18877905 -3.9658338
## Z 2.8030486 -1.990451450 5.7091456 1.54336389 -1.6762998
## F 0.2227348 -2.451980697 2.7746535 0.38613948 1.0527467
## M 2.4172285 -1.323054933 1.3934140 -0.25383584 0.6615951
## T -2.5711818 -1.902815534 1.4933224 3.05620725 -2.0944220
## V 0.7565844 -1.877069964 4.4759921 1.54149311 -1.6243172
## C 2.4963601 -2.202611261 4.7859610 1.65747860 -1.1335972
## phylogenetic PCA with species means
pca<-phyl.pca(tree,Xm)
## scores for species means
scores(pca)
## PC1 PC2 PC3 PC4 PC5
## Y 4.8355633 -0.3294605 4.9204846 -1.214253071 0.48417330
## K 1.1271032 -1.8905319 1.1887193 0.146415414 0.52679413
## R -1.3776735 -0.5685501 1.0264633 -0.240509269 -0.23031699
## O -3.0957196 -3.2959470 4.5902049 0.079157313 0.17733263
## J 3.4199678 -2.3877114 0.5020696 -1.548803070 -0.50276096
## W 4.4388588 -4.0399357 -1.8032739 -1.402434823 0.98914395
## P 5.8491031 -1.8008209 3.2645258 -1.562917383 1.10626368
## N -2.6273140 -3.3123513 0.6706146 0.405660969 1.13010196
## Q -2.4989968 -5.2056862 1.2006453 0.378034447 1.59048227
## G 1.5317465 0.1939072 3.6138101 -0.224632080 -0.07356055
## L 0.3058007 0.7250770 2.1747536 0.317350313 0.15884578
## X 0.6098092 -1.3901896 -3.2212485 -0.563561025 0.04254944
## A 0.8665509 2.2910213 -1.9286367 -0.130411521 0.12191960
## S -2.3410551 0.4294088 2.8283221 0.915182583 0.25544998
## D -2.5983693 0.3485202 1.9828177 1.179837056 -0.11900397
## H -3.5720692 3.4015096 2.7538485 0.649728549 -1.19057929
## E -1.0349547 2.1898149 -2.1836836 0.006607482 -0.14215189
## U -0.5987227 1.1359716 -3.3898069 1.467907338 1.06531799
## I 0.7875657 -0.3033870 -2.6109169 -0.309751715 0.26272155
## B -6.0691018 1.7853416 -5.1025947 1.320105136 -1.02195377
## Z 4.4303034 -1.8337284 -1.3996878 -0.735809021 -0.14132283
## F 1.4969275 -0.1922748 2.0672492 -0.093322338 0.59191066
## M 2.5702312 2.2961278 1.0451950 -0.281481952 0.96142293
## T -2.5011852 -1.2398761 -0.6973528 -1.171601072 -0.35424555
## V 2.1926843 -1.5966486 -0.7306491 -0.654874153 -0.59622447
## C 3.7181929 -1.3800392 -0.8978589 -0.808252464 0.43215101
## scores for individuals
scores(pca,newdata=Xi)
## PC1 PC2 PC3 PC4 PC5
## A1 0.58137458 2.26833022 -1.50992755 0.05254930 -0.05285582
## A2 1.15172725 2.31371230 -2.34734580 -0.31337234 0.29669502
## B1 -5.69194457 1.27868306 -5.63929102 1.53355448 -1.13271412
## B2 -6.21729020 1.56489018 -4.80403093 1.57346089 -1.01085811
## B3 -6.09486112 2.37746282 -4.62585228 0.88981672 -1.46686067
## B4 -6.52217871 1.42718979 -5.04909125 1.66474417 -1.17136154
## B5 -5.81923440 2.27848201 -5.39470803 0.93894942 -0.32797441
## C1 4.04819458 -1.22316584 -0.98201155 -0.31625536 0.20825572
## C2 3.81380896 -1.50051504 -0.44363895 -1.31721174 -0.04694296
## C3 3.83684934 -0.78430503 -1.74656305 -0.80025429 0.76853764
## C4 3.35984901 -1.86986239 -0.46865739 -1.25673076 0.29950120
## C5 3.53226258 -1.52234773 -0.84842353 -0.35081017 0.93140344
## D1 -2.58766894 0.28198061 1.71060049 1.13387065 0.05883451
## D2 -2.60906967 0.41505973 2.25503481 1.22580347 -0.29684245
## E1 -0.91950745 2.94693733 -2.90222653 -0.22396817 -0.48265722
## E2 -0.63006722 2.04596141 -1.75443025 0.03041785 -0.05889452
## E3 -1.55528957 1.57654590 -1.89439396 0.21337277 0.11509608
## F1 1.49692754 -0.19227481 2.06724918 -0.09332234 0.59191066
## G1 1.98319119 0.33592552 3.26622081 -0.07905380 0.30794692
## G2 1.01866478 0.14291709 3.58606835 -1.12319680 -0.21212503
## G3 1.18807446 -0.22799820 4.24271622 -0.19775773 -0.12889419
## G4 1.93705568 0.52478448 3.36023520 0.50148001 -0.26116988
## H1 -4.21463202 3.06152517 2.98057354 0.74927758 -0.77191141
## H2 -3.46355581 3.61819193 2.63532881 0.71375796 -0.86001124
## H3 -4.00177798 3.77949585 2.41281338 0.23074439 -1.70697275
## H4 -3.42305167 3.37912219 2.41521653 0.58097673 -1.71788761
## H5 -2.75732844 3.16921261 3.32531008 0.97388609 -0.89611342
## I1 0.78756568 -0.30338695 -2.61091689 -0.30975171 0.26272155
## J1 3.77009981 -2.29606193 0.82355733 -1.87217506 -0.50930671
## J2 2.65643636 -2.62752255 0.36504020 -1.43004519 -0.74223991
## J3 3.93220043 -1.66531301 0.15517195 -1.12167970 -0.70897156
## J4 3.32113450 -2.96194826 0.66450880 -1.77131234 -0.05052566
## K1 1.25159512 -1.57925153 1.63321774 -0.17957517 -0.15125090
## K2 1.09822059 -1.70247634 0.34879191 0.38728142 1.51661199
## K3 1.15682815 -2.16043574 2.02139355 0.18036386 -0.24669436
## K4 1.34821280 -1.66402160 0.98964668 -0.20557961 0.62257941
## K5 0.78065922 -2.34647404 0.95054657 0.54958657 0.89272452
## L1 1.06414744 0.78680684 2.11315269 -0.27273027 0.06153136
## L2 0.29615528 -0.04188865 1.54620281 1.18255466 -0.30825301
## L3 0.36272587 0.61785755 2.72579932 0.62588492 1.02569981
## L4 -0.08817270 0.98616467 1.53609670 -0.44725043 -0.18624600
## L5 -0.10585249 1.27644467 2.95251641 0.49829269 0.20149675
## M1 2.68309527 2.41206799 1.42442617 -0.73468540 0.22034228
## M2 2.45736721 2.18018762 0.66596382 0.17172150 1.70250357
## N1 -2.56948782 -2.95017073 -0.05817857 0.39753741 1.52137482
## N2 -2.04353648 -3.52769425 1.21271646 1.00911592 1.69237269
## N3 -2.82187513 -3.43066649 0.61806938 0.19573987 1.19061494
## N4 -3.07435651 -3.34087382 0.90985106 0.02025067 0.11604538
## O1 -3.11192122 -2.70534277 4.64768619 0.57463847 0.18529130
## O2 -3.16038176 -3.87005819 4.24883516 -0.59826113 0.27313704
## O3 -3.01485587 -3.31243993 4.87409349 0.26109460 0.07356954
## P1 5.84910311 -1.80082091 3.26452583 -1.56291738 1.10626368
## Q1 -2.13304010 -5.24928996 0.85558856 0.87580561 1.36810917
## Q2 -2.86495354 -5.16208245 1.54570204 -0.11973672 1.81285536
## R1 -1.31674630 -0.90476921 1.96308744 0.04823315 0.15262143
## R2 -1.58422603 -0.05391324 0.78230593 -0.05693333 -1.28256572
## R3 -1.97457519 -1.18864009 0.41523442 0.11054128 -0.43355407
## R4 -0.95479129 0.03372347 0.99863380 -0.49640797 0.85483205
## R5 -1.05802876 -0.72915132 0.97305468 -0.80797947 -0.44291866
## S1 -1.56525324 0.29512563 3.63825338 1.14104633 -0.56828429
## S2 -3.25461074 0.14469106 2.26743678 0.57200356 0.46551160
## S3 -2.12348163 0.75013909 2.58217027 1.29353192 0.45436537
## S4 -2.42087490 0.52767940 2.82542790 0.65414852 0.67020723
## T1 -2.50118518 -1.23987611 -0.69735278 -1.17160107 -0.35424555
## U1 -0.70258767 1.29596288 -3.03986714 1.58536349 1.36163333
## U2 -0.78716461 0.98955203 -3.29104221 2.07335134 1.27860604
## U3 -0.28356663 1.79696944 -3.15636432 1.09327885 1.02709589
## U4 -0.62157182 0.46140212 -4.07195378 1.11963568 0.59393669
## V1 2.19268431 -1.59664865 -0.73064913 -0.65487415 -0.59622447
## W1 4.59323403 -4.57051267 -2.13032675 -0.95022639 0.98692254
## W2 4.28448359 -3.50935871 -1.47622110 -1.85464325 0.99136536
## X1 0.78289693 -1.77122312 -3.41011353 -0.30871842 0.73233760
## X2 0.28395177 -1.12040455 -2.70465518 -0.98695764 -1.02311779
## X3 0.69020005 -1.30387401 -2.91689514 -0.35978315 -0.27484074
## X4 1.37014904 -1.63790746 -4.16308027 -1.03919153 0.05889125
## X5 -0.07815197 -1.11753896 -2.91149857 -0.12315438 0.71947689
## Y1 4.44779172 -0.63479302 5.39272325 -1.36041582 -0.07126708
## Y2 5.65042965 0.20766608 4.69310049 -1.47307135 1.11650405
## Y3 4.40846844 -0.56125449 4.67562993 -0.80927204 0.40728292
## Z1 3.87645446 -2.25958453 -1.32575109 -0.72808203 0.78927693
## Z2 5.27985339 -1.79508555 -0.91340898 -1.16277358 -0.10879099
## Z3 4.35122084 -2.21607210 -1.86618529 -1.03597318 -0.76345277
## Z4 4.54419334 -1.27276565 -1.50322641 0.18997248 -0.36284858
## Z5 4.09979495 -1.62513421 -1.38986713 -0.94218880 -0.26079876
Cool. An interesting property of these data is that when we use the rotation from the species means the mean PC scores from the individual data are the same as the original PC scores from the means:
pca<-phyl.pca(tree,Xm)
S<-scores(pca,newdata=Xi)
Sa<-matrix(NA,26,5)
rownames(Sa)<-LETTERS
for(i in 1:26){
ii<-grep(LETTERS[i],rownames(S))
Sa[i,]<-if(length(ii)>1) colMeans(S[ii,]) else S[ii,]
}
plot(Sa,scores(pca)[rownames(Sa),],cex=1.5,pch=19,
col=make.transparent("darkblue",0.5),
xlab="mean scores",ylab="scores from means")
The data for this exercise were simulated using the packages clusterGeneration & phytools as follows:
library(clusterGeneration)
tree<-pbtree(n=26,scale=1,tip.label=sample(LETTERS))
k<-5
V<-genPositiveDefMat(k,covMethod="unifcorrmat")$Sigma
X<-sim.corrs(phytools:::lambdaTree(tree,0.7),V)[LETTERS,]
N<-sample(1:5,nrow(X),replace=TRUE)
Xi<-apply(X,2,sampleFrom,xvar=0.2,n=N)
spp<-as.factor(rownames(Xi))
rownames(Xi)<-paste(rownames(Xi),
unlist(sapply(N,function(n) 1:n)),
sep="")
That's it.
