In response to some optimization issues a phytools user reported encountered for the
function evolvcv.lite, I just
pushed
a handful of changes, which should both improve optimization, and help users to identify
certain types of optimization failures.
library(phytools)
packageVersion("phytools")
## [1] '1.1.12'
As a reminder, evolvcv.lite fits a set of nested heterogeneous correlated Brownian
evolution models for two quantitative traits, as described in
this recent paper by Ken Toyama, Luke Mahler, and myself.
Here's a quick illustrative demo of how the function works.
In this case, I'm going to load a tree of sunfishes, with feeding mode (piscivory or non-piscivory) mapped onto the edges of the phylogeny. Then I'm going to fit a set of models in which the rate and evolutionary correlation (and both) between two aspects of the buccal morphology are free to vary across the tree according to feeding mode.
Our hypothesis is that piscivory constrains buccal evolution, so there should be a higher evolutionary correlation in piscivorous lineages.
## load tree
data(sunfish.tree)
data(sunfish.data)
## create plot showing data
par(mfrow=c(1,2))
plot(sunfish.tree,ftype="i",fsize=0.8)
## no colors provided. using the following legend:
## non pisc
## "black" "#DF536B"
par(mar=c(5.1,4.1,1.1,1.1),cex.axis=0.8,
las=1)
phylomorphospace(sunfish.tree,
sunfish.data[,2:3],ftype="off",
bty="n",xlab="relative gape width",
ylab="relative buccal length",
node.size=c(0,1.3),node.by.map=TRUE)
legend("topleft",c("non-piscivorous","piscivorous"),
pt.bg=palette()[1:2],pch=22,pt.cex=1.5,
cex=0.8,bty="n")
Next, let's run and print our set of models.
## fit models
sunfish.fit<-evolvcv.lite(sunfish.tree,sunfish.data[,2:3],models="all models")
## Fitting model 1: common rates, common correlation...
## Best log(L) from model 1: 72.1893.
## Fitting model 2: different rates, common correlation...
## Best log(L) from model 2: 77.9869.
## Fitting model 2b: different rates (trait 1), common correlation...
## Best log(L) from model 2b: 75.9661.
## Fitting model 2c: different rates (trait 2), common correlation...
## Best log(L) from model 2c: 75.3309.
## Fitting model 3: common rates, different correlation...
## Best log(L) from model 3: 73.6144.
## Fitting model 3b: different rates (trait 1), different correlation...
## Best log(L) from model 3b: 76.4875.
## Fitting model 3c: different rates (trait 2), different correlation...
## Best log(L) from model 3c: 80.7484.
## Fitting model 4: no common structure...
## Best log(L) from model 4: 81.2393.
print(sunfish.fit)
## Model 1: common rates, common correlation
## R[1,1] R[1,2] R[2,2] k log(L) AIC
## fitted 0.114 0.033 0.0556 5 72.1893 -134.3786
##
## (R thinks it has found the ML solution for model 1.)
##
## Model 2: different rates, common correlation
## R[1,1] R[1,2] R[2,2] k log(L) AIC
## non 0.183 0.0272 0.0199 7 77.9869 -141.9738
## pisc 0.0489 0.0296 0.0889
##
## (R thinks it has found the ML solution for model 2.)
##
## Model 2b: different rates (trait 1), common correlation
## R[1,1] R[1,2] R[2,2] k log(L) AIC
## non 0.2006 0.0578 0.0553 6 75.9661 -139.9321
## pisc 0.044 0.0271 0.0553
##
## (R thinks it has found the ML solution for model 2b.)
##
## Model 2c: different rates (trait 2), common correlation
## R[1,1] R[1,2] R[2,2] k log(L) AIC
## non 0.1142 0.0144 0.0165 6 75.3309 -138.6619
## pisc 0.1142 0.0357 0.1011
##
## (R thinks it has found the ML solution for model 2c.)
##
## Model 3: common rates, different correlation
## R[1,1] R[1,2] R[2,2] k log(L) AIC
## non 0.0991 0.0129 0.0634 6 73.6144 -135.2287
## pisc 0.0991 0.0538 0.0634
##
## (R thinks it has found the ML solution for model 3.)
##
## Model 3b: different rates (trait 1), different correlation
## R[1,1] R[1,2] R[2,2] k log(L) AIC
## non 0.1677 0.0343 0.0556 7 76.4875 -138.975
## pisc 0.0454 0.0324 0.0556
##
## (R thinks it has found the ML solution for model 3b.)
##
## Model 3c: different rates (trait 2), different correlation
## R[1,1] R[1,2] R[2,2] k log(L) AIC
## non 0.1128 -0.003 0.0115 7 80.7484 -147.4968
## pisc 0.1128 0.1138 0.1577
##
## (R thinks it has found the ML solution for model 3c.)
##
## Model 4: no common structure
## R[1,1] R[1,2] R[2,2] k log(L) AIC
## non 0.1386 -0.0022 0.012 8 81.2393 -146.4785
## pisc 0.0758 0.0795 0.129
##
## (R thinks it has found the ML solution for model 4.)
anova(sunfish.fit)
## log(L) d.f. AIC weight
## model 1 72.18932 5 -134.3786 0.0008255187
## model 2 77.98690 7 -141.9738 0.0368124897
## model 2b 75.96605 6 -139.9321 0.0132631624
## model 2c 75.33095 6 -138.6619 0.0070278692
## model 3 73.61437 6 -135.2287 0.0012627698
## model 3b 76.48752 7 -138.9750 0.0082190369
## model 3c 80.74838 7 -147.4968 0.5824932088
## model 4 81.23927 8 -146.4785 0.3500959444
Our best-supported model is model 3c: different rates (trait 2), different correlation.
Let's see the evolutionary correlations by feeding mode.
cov2cor(sunfish.fit$model3c$R$pisc)
## [,1] [,2]
## [1,] 1.0000000 0.8530147
## [2,] 0.8530147 1.0000000
cov2cor(sunfish.fit$model3c$R$non)
## [,1] [,2]
## [1,] 1.00000000 -0.08183551
## [2,] -0.08183551 1.00000000
Indeed, the evolutionary correlation is substantially higher in piscivorous lineages in our best-supported evolutionary scenario for the traits.
One of most significant updates in this function version is an attempt to specify reasonable bounds on optimization, and to indicate if optimization reaches those bounds.
Bounds are set using the arguments lower and upper: and the first two values for each bound
vector are for the rates of the two traits (on a log scale – i.e., -1 means exp(-1)); while
the third elements indicate the lower and upper bounds of the correlation.
Let's see what happens when we set lower to 0 for the correlation. (Here I'll just optimize
the four basic models.)
sunfish.bounded<-evolvcv.lite(sunfish.tree,sunfish.data[,2:3],
lower=c(-4,-4,0))
## Fitting model 1: common rates, common correlation...
## Best log(L) from model 1: 72.1893.
## Fitting model 2: different rates, common correlation...
## Best log(L) from model 2: 77.9869.
## Fitting model 3: common rates, different correlation...
## Best log(L) from model 3: 73.6144.
## Fitting model 4: no common structure...
## Best log(L) from model 4: 80.815.
sunfish.bounded
## Model 1: common rates, common correlation
## R[1,1] R[1,2] R[2,2] k log(L) AIC
## fitted 0.114 0.033 0.0556 5 72.1893 -134.3786
##
## (R thinks it has found the ML solution for model 1.)
##
## Model 2: different rates, common correlation
## R[1,1] R[1,2] R[2,2] k log(L) AIC
## non 0.183 0.0271 0.0199 7 77.9869 -141.9738
## pisc 0.0489 0.0296 0.0889
##
## (R thinks it has found the ML solution for model 2.)
##
## Model 3: common rates, different correlation
## R[1,1] R[1,2] R[2,2] k log(L) AIC
## non 0.0991 0.0129 0.0634 6 73.6144 -135.2287
## pisc 0.0991 0.0538 0.0634
##
## (R thinks it has found the ML solution for model 3.)
##
## Model 4: no common structure
## R[1,1] R[1,2] R[2,2] k log(L) AIC
## non 0.1413 7e-04 0.0183 8 80.815 -145.6299
## pisc 0.0714 0.0729 0.1193
##
## (Model 4 optimization may be at bounds.)
We can see that the first three models are OK, but the fourth one hit the boundary – which is what we'd expect based on our previous analysis.
Now, let's do the same but re-analyzing the dataset from our recent study.
In this analysis, we're going to compare the evolutionary correlation between body size and relative body depth in tropidurid lizards that are rock-dwelling and non-rock-dwelling. We might expect this relationship to become decoupled (or even turn negative) in rock-dwelling forms.
trop.tree<-read.simmap(file=
"https://raw.githubusercontent.com/liamrevell/evolvcv.lite.figures/main/tropidurid-tree.tre",
version=1.5)
trop.data<-read.csv(file=
"https://raw.githubusercontent.com/liamrevell/evolvcv.lite.figures/main/tropidurid-data.csv",
row.names=1)
cols<-setNames(c("white","black"),c("n_rock","rock"))
plot(trop.tree,cols,ftype="i",fsize=0.6,lwd=3,outline=TRUE,
direction="upwards")
legend("bottomright",c("non-rock-dwelling","rock-dwelling"),
pt.bg=cols,pch=22,pt.cex=1.5,cex=0.8,bty="n")
tropidurid.fits<-evolvcv.lite(trop.tree,trop.data,models="all models")
## Fitting model 1: common rates, common correlation...
## Best log(L) from model 1: 52.3056.
## Fitting model 2: different rates, common correlation...
## Best log(L) from model 2: 54.3968.
## Fitting model 2b: different rates (trait 1), common correlation...
## Best log(L) from model 2b: 52.7752.
## Fitting model 2c: different rates (trait 2), common correlation...
## Best log(L) from model 2c: 53.8969.
## Fitting model 3: common rates, different correlation...
## Best log(L) from model 3: 55.1105.
## Fitting model 3b: different rates (trait 1), different correlation...
## Best log(L) from model 3b: 55.3174.
## Fitting model 3c: different rates (trait 2), different correlation...
## Best log(L) from model 3c: 56.0383.
## Fitting model 4: no common structure...
## Best log(L) from model 4: 56.2877.
tropidurid.fits
## Model 1: common rates, common correlation
## R[1,1] R[1,2] R[2,2] k log(L) AIC
## fitted 0.2224 0.0154 0.0589 5 52.3056 -94.6112
##
## (R thinks it has found the ML solution for model 1.)
##
## Model 2: different rates, common correlation
## R[1,1] R[1,2] R[2,2] k log(L) AIC
## n_rock 0.2025 0.0187 0.0456 7 54.3968 -94.7936
## rock 0.3042 0.0382 0.1263
##
## (R thinks it has found the ML solution for model 2.)
##
## Model 2b: different rates (trait 1), common correlation
## R[1,1] R[1,2] R[2,2] k log(L) AIC
## n_rock 0.2029 0.0167 0.0589 6 52.7752 -93.5505
## rock 0.3019 0.0203 0.0589
##
## (R thinks it has found the ML solution for model 2b.)
##
## Model 2c: different rates (trait 2), common correlation
## R[1,1] R[1,2] R[2,2] k log(L) AIC
## n_rock 0.2225 0.0181 0.0457 6 53.8969 -95.7938
## rock 0.2225 0.03 0.1256
##
## (R thinks it has found the ML solution for model 2c.)
##
## Model 3: common rates, different correlation
## R[1,1] R[1,2] R[2,2] k log(L) AIC
## n_rock 0.2256 0.0394 0.0588 6 55.1105 -98.221
## rock 0.2256 -0.0354 0.0588
##
## (R thinks it has found the ML solution for model 3.)
##
## Model 3b: different rates (trait 1), different correlation
## R[1,1] R[1,2] R[2,2] k log(L) AIC
## n_rock 0.2122 0.0372 0.0587 7 55.3174 -96.6348
## rock 0.2738 -0.04 0.0587
##
## (R thinks it has found the ML solution for model 3b.)
##
## Model 3c: different rates (trait 2), different correlation
## R[1,1] R[1,2] R[2,2] k log(L) AIC
## n_rock 0.2251 0.0345 0.0488 7 56.0383 -98.0766
## rock 0.2251 -0.0467 0.1
##
## (R thinks it has found the ML solution for model 3c.)
##
## Model 4: no common structure
## R[1,1] R[1,2] R[2,2] k log(L) AIC
## n_rock 0.2108 0.0325 0.0485 8 56.2877 -96.5753
## rock 0.2794 -0.0564 0.1009
##
## (R thinks it has found the ML solution for model 4.)
anova(tropidurid.fits)
## log(L) d.f. AIC weight
## model 1 52.30560 5 -94.61119 0.04619846
## model 2 54.39681 7 -94.79362 0.05061060
## model 2b 52.77525 6 -93.55050 0.02718320
## model 2c 53.89688 6 -95.79376 0.08344880
## model 3 55.11048 6 -98.22096 0.28085276
## model 3b 55.31738 7 -96.63476 0.12706885
## model 3c 56.03828 7 -98.07656 0.26129020
## model 4 56.28765 8 -96.57530 0.12334712
We see that all the model support is divided among the four models (model 3, 3b, 3c, and 4)
that allow the evolutionary correlation to differ between lineages.
That's it.
No comments:
Post a Comment
Note: due to the very large amount of spam, all comments are now automatically submitted for moderation.