Possible models:
model=BM1 # a Brownian motion model applied over the entire tree
model=BMS # a different Brownian rate parameter (sigma2) for each regime
model=OU1 # a single peak Ornstein-Uhlenbeck model across the entire tree
model=OUM # a multi-peak Ornstein-Uhlenbeck model with different optima (theta) for each regime
model=OUMV # a multi-peak Ornstein-Uhlenbeck model with different optima (theta) and different Brownian rate parameter (sigma2) for each regime
model=OUMA # a multi-peak Ornstein-Uhlenbeck model with different optima (theta) and different strength of selection parameter (alpha) for each regime
model=OUMVA # a multi-peak Ornstein-Uhlenbeck model with different optima (theta), different Brownian rate parameter (sigma2), and different strength of selection parameter (alpha) for each regime
A great step forward for all of us using comparative methods. To start playing, all we need to specify are a tree, a data frame, and a model. Katie and I had already been playing around with make.simmap from the package phytools, so we had trees with stochastically mapped discrete traits. The data frame needs to be in the format where the first column is the species names, the second column is the discrete character, and the third column is the continuous character. Specifying that can be as easy as this:
> df1 = data.frame(speciesNames, discrete, PC1)
We had multiple mapped trees (as we should), so we took one to go for a spin, starting with the simplest model first:
> bm1Output = OUwie(mouthPsimmap[[1]], df1, model = "BM1", simmap.tree = TRUE)
Begin subplex optimization routine -- Starting value: 1
Finished. Summarizing results.
Finished. Performing diagnostic tests.
> bmsOutput = OUwie(mouthPsimmap[[1]], df1, model = "BMS", simmap.tree = TRUE)
Begin subplex optimization routine -- Starting value: 1
Finished. Summarizing results.
Finished. Performing diagnostic tests.
For this model, though, we need to be sure to iterate over the different mappings as this depends on the selective regime. We should get different results from different trees as we can see here:
> bmsOutput2 = OUwie(mouthPsimmap[[2]], df1, model = "BMS", simmap.tree = TRUE)
Begin subplex optimization routine -- Starting value: 1
Finished. Summarizing results.
Finished. Performing diagnostic tests.
> bmsOutput$loglik
[1] -30.24609
> bmsOutput2$loglik
[1] -29.44005
We have a multiphylo of 500 mapped trees, so let's just sample 20 of them to see what's going on.
> sampledTrees = sample(mouthPsimmap, 20)
Now, we can just throw this into an apply to get a list of output for each of these twenty trees.
> bmsSampledOutput = lapply(sampledTrees, function(x) {
+ OUwie(x, df1, model = "BMS", simmap.tree = TRUE)
+ })
This might take a while, but once it's done, everything seems great! We can use some more apply functions to get out particular bits of the data we want.
Unfortunately, as soon as I add varying alpha into my model, I run into some problems:
> oumaOutput1 = OUwie(mouthPsimmap[[1]], df1, model = "OUMA", simmap.tree = TRUE)
Begin subplex optimization routine -- Starting value: 1
Error in solve.default(V, res) :
system is computationally singular: reciprocal condition number = 1.27328e-103
I've tried a variety of debugging strategies, including changing the optimization algorithm in the OUwie function, but I can't seem to get anything that varies alpha to work, unfortunately. This is probably a function of having very little data, and thus no power to estimate parameters. Unfortunate, but I would love to try this again on a larger dataset!
Hi Tomomi~
ReplyDeleteThis error has been fixed in a newer version of OUwie. Feel free to send me an email if you are interested.
All the best,
Jeremy