Wednesday, August 21, 2019

Computing principal components scores for new data with phyl.pca

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")

plot of chunk unnamed-chunk-2

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.

Wednesday, July 31, 2019

Stochastic character mapping with polymorphic characters using make.simmap and fitpolyMk in phytools

I recently received an email with the following inquiry:

“I wondered if you had any insight or advice on how the make.simmap function in phytools works when the input data have multiple rows per species? It certainly does run but I am unclear on how example it's dealing with this. Specifically, I'm working with a dataset of a 3-state categorical trait but for which I have between 1 and 42 (mean 3) datapoints per species …. I initially fed in the trait data directly to make.simmap, and it did run fine and give reasonable looking estimates for internal nodes. However, it also allocated a single state (without exception) to the tips, so I'm not sure if this is choosing the modal value, the first value, or something else here?

As an alternative, I converted the raw data to a matrix of proportions of each state per species (in a similar manner as incorporating uncertainty in tip state as you described here) based on the observed proportions for each species. This gives very similar estimates for ancestral states, perhaps within the bounds of the inherent stochasticity of the method, though different enough in places that I was suspicious that it's doing something different from my first approach. It also (of course) gives tip states that are as variable as the data it derived from, as you'd expect. However I wonder whether this approach, which essentially equates uncertainty in tip state with intraspecific variability, is a valid way of treating this data.

In summary, do you have anything you can add concerning the appropriateness of either of these approaches (or any other suggestions) for dealing with N-state categorical traits which are variable within species?”

First - in the former case make.simmap is most definitely not taking intraspecific polymorphism into account. Rather, as the user suspected, it is most-likely just using the first observation for each species.

The second alternative is more interesting. Although this may seem reasonable it is also not correct. Saying that a particular tip taxa has a prior probability of being in state A that is P(A)=0.50 and P(B)=0.50 implies that the taxon could be either state A or B with equal probability - not that (say) 50% of the population is in state A & the other 50% in state B.

Though this may be somewhat unsatisfying, the simplest way to deal with intraspecific polymorphism is to re-code our character as an additional state. That is, instead of having character states A and B, we instead have A, A+B, and B.

Unfortunately, phytools::make.simmap presently doesn't automatically know how to handle data in this form. The function fitpolyMk does, however, and it is relatively straightforward to integrate it into make.simmap. I'll next demonstrate how.

Note that I just pushed a very small update to fitpolyMk to make this easier - thus to follow this demo exactly it is necessary to update phytools from GitHub.

Let's get some data & a tree for our demo. Here are mine:

library(phytools)
packageVersion("phytools")
## [1] '0.7.2'
poly.data<-read.csv(file="http://www.phytools.org/TS2019/ex/4c/polymorphic-data.csv",
    row.names=1)
y<-setNames(poly.data[,1],rownames(poly.data))
y
##  t70  t71  t47 t169 t170  t67  t68  t62  t59  t60 t173 t174  t19  t20  t21 
##  C+D  C+D  C+D  C+D  C+D  C+D  B+C  C+D    C    C    B    B    B  B+C    B 
##  t15  t50 t109 t123 t124 t104 t105  t52  t53 t167 t168 t141  t74  t66  t48 
##    B  B+C  B+C  B+C  B+C  B+C  B+C  B+C  B+C    C    C    C    C  C+D    C 
## t181 t182   t2  t39  t40   t4  t27 t191 t192  t28 t145 t146  t77   t8 t193 
##  C+D  C+D    C    C    C    D    C  B+C    C    D  C+D  C+D  C+D    D    D 
## t194  t96  t45  t86  t87  t13   t6 t125 t183 t184  t16  t17 t135 t136 t147 
##    D    D    D    D    D  C+D  C+D  C+D    C    C    C    C    C    C    C 
## t148  t83  t88  t89  t78  t79  t24 t187 t188  t36  t29  t99 t189 t190 t115 
##    C    C  B+C    C    C    C    C    C    C    C    D  C+D    C    C    C 
## t119 t120 t121 t122  t46  t49 t165 t166 t112   t5 t106 t126 t142 t143  t98 
##    C    C    C    C    C    C  C+D  C+D  C+D  B+C  B+C  B+C  B+C  B+C    C 
## t179 t180 t144 t113 t114  t18  t30 t155 t156  t11 t128 t129 t127  t64  t65 
##  B+C  B+C  B+C  B+C  B+C  B+C    C    C    C  B+C    C  C+D    C  C+D  C+D 
## t153 t154 t185 t186  t80  t81  t33  t43 t151 t152  t82  t97 t132 t133  t61 
##    C    C    C    C    B    B    C    D  C+D  C+D  C+D    D    D    D    D 
##   t3 t102 t103  t35  t69 t149 t150   t7  t92  t93 t110 t111 t197 t198 t161 
##  A+B  B+C  B+C  B+C  A+B    B    B  C+D  C+D  C+D  C+D  C+D    C    C    C 
## t157 t158  t12   t9  t41  t42  t25  t37  t38  t26  t10  t90  t94  t95  t54 
##    C    C  B+C    C    C    C    C    C    C    C  B+C  C+D  C+D  C+D    C 
##  t55  t91 t171 t172 t162 t163  t51  t32 t177 t178  t63  t84  t85  t44  t56 
##    C    D    C    C  C+D  C+D  C+D    C    C    C  B+C    C    C    C    C 
## t199 t200 t107 t108 t130 t131  t31 t137 t138  t22 t139 t140 t118  t72  t73 
##    C    C    C    C    C    C    C  C+D    C    C  C+D  C+D  C+D  C+D    C 
##  t76 t159 t160  t14  t75 t134 t195 t196 t164  t34 t100 t101 t116 t117  t23 
##  C+D  C+D  C+D    D  C+D  C+D  C+D  C+D  C+D    D    D    D  C+D    C  C+D 
## t175 t176  t57  t58   t1 
##  C+D  C+D  B+C  B+C  C+D 
## Levels: A+B B B+C C C+D D
poly.tree<-read.tree(file="http://www.phytools.org/TS2019/ex/4c/polymorphic-tree.phy")

Just for fun - and to show ourselves that fitpolyMk can handle it - let's flip the order of some of our polymorphic conditions. In other words, we will make state A + B a mixture of taxa coded as A+B and B+A and so on.

y<-setNames(as.character(y),names(y))
ii<-which(y=="A+B")
y[sample(ii,round(0.5*length(ii)))]<-"B+A"
ii<-which(y=="C+D")
y[sample(ii,round(0.5*length(ii)))]<-"D+C"
y<-as.factor(y)
y
##  t70  t71  t47 t169 t170  t67  t68  t62  t59  t60 t173 t174  t19  t20  t21 
##  D+C  C+D  C+D  C+D  D+C  C+D  B+C  C+D    C    C    B    B    B  B+C    B 
##  t15  t50 t109 t123 t124 t104 t105  t52  t53 t167 t168 t141  t74  t66  t48 
##    B  B+C  B+C  B+C  B+C  B+C  B+C  B+C  B+C    C    C    C    C  C+D    C 
## t181 t182   t2  t39  t40   t4  t27 t191 t192  t28 t145 t146  t77   t8 t193 
##  C+D  D+C    C    C    C    D    C  B+C    C    D  D+C  C+D  D+C    D    D 
## t194  t96  t45  t86  t87  t13   t6 t125 t183 t184  t16  t17 t135 t136 t147 
##    D    D    D    D    D  C+D  C+D  C+D    C    C    C    C    C    C    C 
## t148  t83  t88  t89  t78  t79  t24 t187 t188  t36  t29  t99 t189 t190 t115 
##    C    C  B+C    C    C    C    C    C    C    C    D  C+D    C    C    C 
## t119 t120 t121 t122  t46  t49 t165 t166 t112   t5 t106 t126 t142 t143  t98 
##    C    C    C    C    C    C  D+C  C+D  C+D  B+C  B+C  B+C  B+C  B+C    C 
## t179 t180 t144 t113 t114  t18  t30 t155 t156  t11 t128 t129 t127  t64  t65 
##  B+C  B+C  B+C  B+C  B+C  B+C    C    C    C  B+C    C  D+C    C  C+D  C+D 
## t153 t154 t185 t186  t80  t81  t33  t43 t151 t152  t82  t97 t132 t133  t61 
##    C    C    C    C    B    B    C    D  C+D  D+C  D+C    D    D    D    D 
##   t3 t102 t103  t35  t69 t149 t150   t7  t92  t93 t110 t111 t197 t198 t161 
##  A+B  B+C  B+C  B+C  B+A    B    B  D+C  D+C  D+C  D+C  D+C    C    C    C 
## t157 t158  t12   t9  t41  t42  t25  t37  t38  t26  t10  t90  t94  t95  t54 
##    C    C  B+C    C    C    C    C    C    C    C  B+C  D+C  C+D  C+D    C 
##  t55  t91 t171 t172 t162 t163  t51  t32 t177 t178  t63  t84  t85  t44  t56 
##    C    D    C    C  D+C  C+D  D+C    C    C    C  B+C    C    C    C    C 
## t199 t200 t107 t108 t130 t131  t31 t137 t138  t22 t139 t140 t118  t72  t73 
##    C    C    C    C    C    C    C  D+C    C    C  D+C  C+D  C+D  D+C    C 
##  t76 t159 t160  t14  t75 t134 t195 t196 t164  t34 t100 t101 t116 t117  t23 
##  D+C  C+D  D+C    D  C+D  C+D  C+D  C+D  D+C    D    D    D  D+C    C  D+C 
## t175 t176  t57  t58   t1 
##  D+C  D+C  B+C  B+C  D+C 
## Levels: A+B B B+A B+C C C+D D D+C

Plot our data:

plotTree(poly.tree,ftype="off",lwd=1,type="fan")
X<-strsplit(setNames(as.character(y),names(y)),"+",fixed=TRUE)
pies<-matrix(0,Ntip(poly.tree),4,dimnames=list(poly.tree$tip.label,
    c("A","B","C","D")))
for(i in 1:Ntip(poly.tree)) 
    pies[poly.tree$tip.label[i],X[[poly.tree$tip.label[i]]]]<-
        rep(1/length(X[[poly.tree$tip.label[i]]]),
        length(X[[poly.tree$tip.label[i]]]))
tiplabels(pie=pies,piecol=c("black","yellow","red","blue"),
    cex=0.35)
legend(x="topleft",legend=c("A","B","C","D"),pt.cex=2,pch=21,
    pt.bg=c("black","yellow","red","blue"))

plot of chunk unnamed-chunk-3

Now we'll go ahead & fit our polymorphic evolution model. Here, I'm going to use an symmetric-rates ordered model - just because I happen to know that this model has good fit to the data:

fit<-fitpolyMk(poly.tree,y,model="SYM",ordered=TRUE)
## 
## This is the design matrix of the fitted model. Does it make sense?
## 
##     A A+B B B+C C C+D D
## A   0   1 0   0 0   0 0
## A+B 1   0 2   0 0   0 0
## B   0   2 0   3 0   0 0
## B+C 0   0 3   0 4   0 0
## C   0   0 0   4 0   5 0
## C+D 0   0 0   0 5   0 6
## D   0   0 0   0 0   6 0
plot(fit,lwd=2)

plot of chunk unnamed-chunk-4

Now, how can we apply this to stochastic character mapping? Well, this is easier than you'd think simply because the function fitpolyMk (now) returns both the design matrix of the model, and a data matrix of the input data to the user as components of the fitted model object: $index.matrix and $data, respectively. Let's examine these & then use these as follows:

fit$index.matrix
##     A A+B B B+C C C+D D
## A   0   1 0   0 0   0 0
## A+B 1   0 2   0 0   0 0
## B   0   2 0   3 0   0 0
## B+C 0   0 3   0 4   0 0
## C   0   0 0   4 0   5 0
## C+D 0   0 0   0 5   0 6
## D   0   0 0   0 0   6 0
fit$data
##      A A+B B B+C C C+D D
## t70  0   0 0   0 0   1 0
## t71  0   0 0   0 0   1 0
## t47  0   0 0   0 0   1 0
## t169 0   0 0   0 0   1 0
## t170 0   0 0   0 0   1 0
## t67  0   0 0   0 0   1 0
## t68  0   0 0   1 0   0 0
## t62  0   0 0   0 0   1 0
## t59  0   0 0   0 1   0 0
## t60  0   0 0   0 1   0 0
## t173 0   0 1   0 0   0 0
## t174 0   0 1   0 0   0 0
## t19  0   0 1   0 0   0 0
## t20  0   0 0   1 0   0 0
## t21  0   0 1   0 0   0 0
## t15  0   0 1   0 0   0 0
## t50  0   0 0   1 0   0 0
## t109 0   0 0   1 0   0 0
## t123 0   0 0   1 0   0 0
## t124 0   0 0   1 0   0 0
## t104 0   0 0   1 0   0 0
## t105 0   0 0   1 0   0 0
## t52  0   0 0   1 0   0 0
## t53  0   0 0   1 0   0 0
## t167 0   0 0   0 1   0 0
## t168 0   0 0   0 1   0 0
## t141 0   0 0   0 1   0 0
## t74  0   0 0   0 1   0 0
## t66  0   0 0   0 0   1 0
## t48  0   0 0   0 1   0 0
## t181 0   0 0   0 0   1 0
## t182 0   0 0   0 0   1 0
## t2   0   0 0   0 1   0 0
## t39  0   0 0   0 1   0 0
## t40  0   0 0   0 1   0 0
## t4   0   0 0   0 0   0 1
## t27  0   0 0   0 1   0 0
## t191 0   0 0   1 0   0 0
## t192 0   0 0   0 1   0 0
## t28  0   0 0   0 0   0 1
## t145 0   0 0   0 0   1 0
## t146 0   0 0   0 0   1 0
## t77  0   0 0   0 0   1 0
## t8   0   0 0   0 0   0 1
## t193 0   0 0   0 0   0 1
## t194 0   0 0   0 0   0 1
## t96  0   0 0   0 0   0 1
## t45  0   0 0   0 0   0 1
## t86  0   0 0   0 0   0 1
## t87  0   0 0   0 0   0 1
## t13  0   0 0   0 0   1 0
## t6   0   0 0   0 0   1 0
## t125 0   0 0   0 0   1 0
## t183 0   0 0   0 1   0 0
## t184 0   0 0   0 1   0 0
## t16  0   0 0   0 1   0 0
## t17  0   0 0   0 1   0 0
## t135 0   0 0   0 1   0 0
## t136 0   0 0   0 1   0 0
## t147 0   0 0   0 1   0 0
## t148 0   0 0   0 1   0 0
## t83  0   0 0   0 1   0 0
## t88  0   0 0   1 0   0 0
## t89  0   0 0   0 1   0 0
## t78  0   0 0   0 1   0 0
## t79  0   0 0   0 1   0 0
## t24  0   0 0   0 1   0 0
## t187 0   0 0   0 1   0 0
## t188 0   0 0   0 1   0 0
## t36  0   0 0   0 1   0 0
## t29  0   0 0   0 0   0 1
## t99  0   0 0   0 0   1 0
## t189 0   0 0   0 1   0 0
## t190 0   0 0   0 1   0 0
## t115 0   0 0   0 1   0 0
## t119 0   0 0   0 1   0 0
## t120 0   0 0   0 1   0 0
## t121 0   0 0   0 1   0 0
## t122 0   0 0   0 1   0 0
## t46  0   0 0   0 1   0 0
## t49  0   0 0   0 1   0 0
## t165 0   0 0   0 0   1 0
## t166 0   0 0   0 0   1 0
## t112 0   0 0   0 0   1 0
## t5   0   0 0   1 0   0 0
## t106 0   0 0   1 0   0 0
## t126 0   0 0   1 0   0 0
## t142 0   0 0   1 0   0 0
## t143 0   0 0   1 0   0 0
## t98  0   0 0   0 1   0 0
## t179 0   0 0   1 0   0 0
## t180 0   0 0   1 0   0 0
## t144 0   0 0   1 0   0 0
## t113 0   0 0   1 0   0 0
## t114 0   0 0   1 0   0 0
## t18  0   0 0   1 0   0 0
## t30  0   0 0   0 1   0 0
## t155 0   0 0   0 1   0 0
## t156 0   0 0   0 1   0 0
## t11  0   0 0   1 0   0 0
## t128 0   0 0   0 1   0 0
## t129 0   0 0   0 0   1 0
## t127 0   0 0   0 1   0 0
## t64  0   0 0   0 0   1 0
## t65  0   0 0   0 0   1 0
## t153 0   0 0   0 1   0 0
## t154 0   0 0   0 1   0 0
## t185 0   0 0   0 1   0 0
## t186 0   0 0   0 1   0 0
## t80  0   0 1   0 0   0 0
## t81  0   0 1   0 0   0 0
## t33  0   0 0   0 1   0 0
## t43  0   0 0   0 0   0 1
## t151 0   0 0   0 0   1 0
## t152 0   0 0   0 0   1 0
## t82  0   0 0   0 0   1 0
## t97  0   0 0   0 0   0 1
## t132 0   0 0   0 0   0 1
## t133 0   0 0   0 0   0 1
## t61  0   0 0   0 0   0 1
## t3   0   1 0   0 0   0 0
## t102 0   0 0   1 0   0 0
## t103 0   0 0   1 0   0 0
## t35  0   0 0   1 0   0 0
## t69  0   1 0   0 0   0 0
## t149 0   0 1   0 0   0 0
## t150 0   0 1   0 0   0 0
## t7   0   0 0   0 0   1 0
## t92  0   0 0   0 0   1 0
## t93  0   0 0   0 0   1 0
## t110 0   0 0   0 0   1 0
## t111 0   0 0   0 0   1 0
## t197 0   0 0   0 1   0 0
## t198 0   0 0   0 1   0 0
## t161 0   0 0   0 1   0 0
## t157 0   0 0   0 1   0 0
## t158 0   0 0   0 1   0 0
## t12  0   0 0   1 0   0 0
## t9   0   0 0   0 1   0 0
## t41  0   0 0   0 1   0 0
## t42  0   0 0   0 1   0 0
## t25  0   0 0   0 1   0 0
## t37  0   0 0   0 1   0 0
## t38  0   0 0   0 1   0 0
## t26  0   0 0   0 1   0 0
## t10  0   0 0   1 0   0 0
## t90  0   0 0   0 0   1 0
## t94  0   0 0   0 0   1 0
## t95  0   0 0   0 0   1 0
## t54  0   0 0   0 1   0 0
## t55  0   0 0   0 1   0 0
## t91  0   0 0   0 0   0 1
## t171 0   0 0   0 1   0 0
## t172 0   0 0   0 1   0 0
## t162 0   0 0   0 0   1 0
## t163 0   0 0   0 0   1 0
## t51  0   0 0   0 0   1 0
## t32  0   0 0   0 1   0 0
## t177 0   0 0   0 1   0 0
## t178 0   0 0   0 1   0 0
## t63  0   0 0   1 0   0 0
## t84  0   0 0   0 1   0 0
## t85  0   0 0   0 1   0 0
## t44  0   0 0   0 1   0 0
## t56  0   0 0   0 1   0 0
## t199 0   0 0   0 1   0 0
## t200 0   0 0   0 1   0 0
## t107 0   0 0   0 1   0 0
## t108 0   0 0   0 1   0 0
## t130 0   0 0   0 1   0 0
## t131 0   0 0   0 1   0 0
## t31  0   0 0   0 1   0 0
## t137 0   0 0   0 0   1 0
## t138 0   0 0   0 1   0 0
## t22  0   0 0   0 1   0 0
## t139 0   0 0   0 0   1 0
## t140 0   0 0   0 0   1 0
## t118 0   0 0   0 0   1 0
## t72  0   0 0   0 0   1 0
## t73  0   0 0   0 1   0 0
## t76  0   0 0   0 0   1 0
## t159 0   0 0   0 0   1 0
## t160 0   0 0   0 0   1 0
## t14  0   0 0   0 0   0 1
## t75  0   0 0   0 0   1 0
## t134 0   0 0   0 0   1 0
## t195 0   0 0   0 0   1 0
## t196 0   0 0   0 0   1 0
## t164 0   0 0   0 0   1 0
## t34  0   0 0   0 0   0 1
## t100 0   0 0   0 0   0 1
## t101 0   0 0   0 0   0 1
## t116 0   0 0   0 0   1 0
## t117 0   0 0   0 1   0 0
## t23  0   0 0   0 0   1 0
## t175 0   0 0   0 0   1 0
## t176 0   0 0   0 0   1 0
## t57  0   0 0   1 0   0 0
## t58  0   0 0   1 0   0 0
## t1   0   0 0   0 0   1 0
simmap.trees<-make.simmap(poly.tree,fit$data,model=fit$index.matrix,nsim=100)
## make.simmap is sampling character histories conditioned on the transition matrix
## 
## Q =
##     A       A+B         B       B+C         C       C+D        D
## A   0  0.000000  0.000000  0.000000  0.000000  0.000000  0.00000
## A+B 0 -1.271756  1.271756  0.000000  0.000000  0.000000  0.00000
## B   0  1.271756 -2.333135  1.061379  0.000000  0.000000  0.00000
## B+C 0  0.000000  1.061379 -1.768122  0.706743  0.000000  0.00000
## C   0  0.000000  0.000000  0.706743 -1.913028  1.206285  0.00000
## C+D 0  0.000000  0.000000  0.000000  1.206285 -2.252835  1.04655
## D   0  0.000000  0.000000  0.000000  0.000000  1.046550 -1.04655
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##         A       A+B         B       B+C         C       C+D         D 
## 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571
## Done.
simmap.trees
## 100 phylogenetic trees with mapped discrete characters

We can plot a summary of this analysis as follows:

cols<-setNames(colorRampPalette(c("black","yellow","red","blue"))(7),
    colnames(fit$data))
plot(summary(simmap.trees),colors=cols,type="fan",ftype="off")
legend(x="topleft",legend=colnames(fit$data),pt.cex=2.4,pch=21,
    pt.bg=cols)

plot of chunk unnamed-chunk-6

However, it might make more sense to do things like combine states before computing posterior probabilities at node. The following is just an example that combines A+B, B, and B+C and then plots them compared to a new state not-B:

b.trees<-mergeMappedStates(simmap.trees,c("A+B","B","B+C"),"B")
b.trees<-mergeMappedStates(b.trees,c("A","C","C+D","D"),"not-B")
cols<-setNames(c("black","yellow"),c("not-B","B"))
plot(summary(b.trees),colors=cols,type="fan",ftype="off")
legend(x="topleft",legend=names(cols),pt.cex=2.4,pch=21,
    pt.bg=cols)

plot of chunk unnamed-chunk-7

Or, similarly using densityMap:

object<-densityMap(b.trees,states=names(cols),plot=FALSE)
## sorry - this might take a while; please be patient
object<-setMap(object,c("black","yellow"))
plot(object,lwd=c(2,6),ftype="off")

plot of chunk unnamed-chunk-8

That's kind of the idea. We can similarly analyze our posterior sample of trees in all of the usual ways, e.g.:

object<-density(simmap.trees)
object
## 
## Distribution of changes from stochastic mapping:
## 
##  A+B->B      A+B->B+C
##  Min.   :0   Min.   :0
##  Median :0   Median :0
##  Mean   :0.34    Mean   :0
##  Max.   :2   Max.   :0
## 
##  A+B->C      A+B->C+D
##  Min.   :0   Min.   :0
##  Median :0   Median :0
##  Mean   :0   Mean   :0
##  Max.   :0   Max.   :0
## 
##  A+B->D      B->A+B
##  Min.   :0   Min.   :2
##  Median :0   Median :2
##  Mean   :0   Mean   :2.29
##  Max.   :0   Max.   :4
## 
##  B->B+C      B->C
##  Min.   :0   Min.   :0
##  Median :1.5 Median :0
##  Mean   :1.85    Mean   :0
##  Max.   :6   Max.   :0
## 
##  B->C+D      B->D
##  Min.   :0   Min.   :0
##  Median :0   Median :0
##  Mean   :0   Mean   :0
##  Max.   :0   Max.   :0
## 
##  B+C->A+B        B+C->B
##  Min.   :0   Min.   :4
##  Median :0   Median :7
##  Mean   :0   Mean   :6.67
##  Max.   :0   Max.   :10
## 
##  B+C->C      B+C->C+D
##  Min.   :1   Min.   :0
##  Median :3   Median :0
##  Mean   :3.54    Mean   :0
##  Max.   :9   Max.   :0
## 
##  B+C->D      C->A+B
##  Min.   :0   Min.   :0
##  Median :0   Median :0
##  Mean   :0   Mean   :0
##  Max.   :0   Max.   :0
## 
##  C->B        C->B+C
##  Min.   :0   Min.   :13
##  Median :0   Median :15
##  Mean   :0   Mean   :15.49
##  Max.   :0   Max.   :21
## 
##  C->C+D      C->D
##  Min.   :18  Min.   :0
##  Median :24  Median :0
##  Mean   :23.9    Mean   :0
##  Max.   :30  Max.   :0
## 
##  C+D->A+B        C+D->B
##  Min.   :0   Min.   :0
##  Median :0   Median :0
##  Mean   :0   Mean   :0
##  Max.   :0   Max.   :0
## 
##  C+D->B+C        C+D->C
##  Min.   :0   Min.   :8
##  Median :0   Median :14
##  Mean   :0   Mean   :13.68
##  Max.   :0   Max.   :26
## 
##  C+D->D      D->A+B
##  Min.   :7   Min.   :0
##  Median :11  Median :0
##  Mean   :10.64   Mean   :0
##  Max.   :15  Max.   :0
## 
##  D->B        D->B+C
##  Min.   :0   Min.   :0
##  Median :0   Median :0
##  Mean   :0   Mean   :0
##  Max.   :0   Max.   :0
## 
##  D->C        D->C+D
##  Min.   :0   Min.   :0
##  Median :0   Median :3
##  Mean   :0   Mean   :3.51
##  Max.   :0   Max.   :9
## 
## 95% HPD interval(A+B->B): [0, 2]
## 95% HPD interval(A+B->B+C): [0, 0]
## 95% HPD interval(A+B->C): [0, 0]
## 95% HPD interval(A+B->C+D): [0, 0]
## 95% HPD interval(A+B->D): [0, 0]
## 95% HPD interval(B->A+B): [2, 4]
## 95% HPD interval(B->B+C): [0, 5]
## 95% HPD interval(B->C): [0, 0]
## 95% HPD interval(B->C+D): [0, 0]
## 95% HPD interval(B->D): [0, 0]
## 95% HPD interval(B+C->A+B): [0, 0]
## 95% HPD interval(B+C->B): [5, 9]
## 95% HPD interval(B+C->C): [1, 7]
## 95% HPD interval(B+C->C+D): [0, 0]
## 95% HPD interval(B+C->D): [0, 0]
## 95% HPD interval(C->A+B): [0, 0]
## 95% HPD interval(C->B): [0, 0]
## 95% HPD interval(C->B+C): [13, 18]
## 95% HPD interval(C->C+D): [20, 29]
## 95% HPD interval(C->D): [0, 0]
## 95% HPD interval(C+D->A+B): [0, 0]
## 95% HPD interval(C+D->B): [0, 0]
## 95% HPD interval(C+D->B+C): [0, 0]
## 95% HPD interval(C+D->C): [8, 19]
## 95% HPD interval(C+D->D): [8, 14]
## 95% HPD interval(D->A+B): [0, 0]
## 95% HPD interval(D->B): [0, 0]
## 95% HPD interval(D->B+C): [0, 0]
## 95% HPD interval(D->C): [0, 0]
## 95% HPD interval(D->C+D): [0, 7]
plot(object)

plot of chunk unnamed-chunk-9

Or from our B or not-B analysis:

object<-density(b.trees)
object
## 
## Distribution of changes from stochastic mapping:
##  B->not-B        not-B->B
##  Min.   :1   Min.   :13
##  Median :3   Median :15
##  Mean   :3.54    Mean   :15.49
##  Max.   :9   Max.   :21
## 
## 95% HPD interval(B->not-B): [1, 7]
## 95% HPD interval(not-B->B): [13, 18]
par(family="mono")
plot(object)

plot of chunk unnamed-chunk-10

You get the idea.

The very heterogeneous sampling (“between 1 and 42 (mean 3) datapoints per species”) raises another interesting point. What if we only have one sample for a species? How can we genuinely say that it is monomorphic? The truth is - we can't; however, if we have a model for our uncertainty we can easily incorporate it in the same way that we specify uncertainty in the values of monomorphic tip states (as previously discussed).

That's all I have to say about this - for now.