Fitting and Comparing Evolutionary Models: Continuous Data
Dean Adams (dcadams@iastate.edu)
2022-04-14
Overview
Today we explore methods for comparing continuous models of trait evolution. Many ‘classic’ macroevolutionary hypotheses are based on a particular model of trait change: does the trait evolve rapidly or slowly? Does a trait evolve faster in one sub-lineage as compared to another? Does trait variation evolve quickly early in the history of a lineage? Does a trait evolve to one or more phenotypic optima?
All of these questions can be addressed by fitting different evolutionary models to the data, and determining which displays the highest support. The approach falls squarely in statistical Model Comparison methods, where the fit of different models is obtained, and evaluated using various approaches: likelihood ratio tests (LRT) and AIC comparisons being the most common.
IMPORTANT STATISTICAL NOTE: Statistical model comparison is of the form: \(Y=X\beta + E\), where the data are fit to different \(X\) variables (i.e., differing explanatory hypotheses). With evolutionary model comparison, all models are the algebraic form: \(Y=1\beta + E\), but where \(E\) is modeled using different expected covariance structures (BM, OU1, etc.). Thus the latter fits different error structures to an intercept model, rather than comparing the fit of the data to differing explanatory variables (see Adams and Collyer 2018; 2019).
Learning Objectives
- Appreciate the diversity of trait evolution models and understand the basic concepts of BM, OU, and EB models.
- Be able to fit BM, OU, and EB models to continuous data.
- Compare which models fit best using AIC or log likelihood
- Understand how to set up models where evolution dynamics differ based on regime shifts
Download data files
Today we will use data from the following files:
Additionally, The R-script for running the code found in this tutorial is found here:
As before, on your local computer, place them in a folder named ‘TutorialData’, which is found in the same directory as this tutorial.
Fitting Evolutionary Models
Here we explore methods for fitting evolutionary models to continuous characters. First, we bring our data into R:
library(phytools)## Loading required package: ape## Loading required package: mapslibrary(geiger)
library(OUwie)## Loading required package: corpcor## Loading required package: nloptr## Loading required package: RColorBrewertree<-read.tree("../data/anole.gp.tre",tree.names=T)
group<-read.csv('../data/anole.gp.csv', row.names=1, header=TRUE,colClasses=c('factor'))
gp<-as.factor(t(group)); names(gp)<-row.names(group)
svl<-read.csv('../data/anole.svl2.csv', row.names=1, header=TRUE)
svl<-as.matrix(treedata(phy = tree,data = svl, warnings=FALSE)$data)[,1] #match data to tree##BM plot data
tree.col<-contMap(tree,svl,plot=FALSE) #runs Anc. St. Est. on branches of tree
plot(tree.col,legend=0.7*max(nodeHeights(tree)),
fsize=c(0.7,0.9))
#PLOT data
tree.col<-contMap(tree,svl,plot=FALSE) #runs Anc. St. Est. on branches of tree
plot(tree.col,type="fan",legend=0.7*max(nodeHeights(tree)),
fsize=c(0.7,0.9))
cols<-setNames(palette()[1:length(unique(gp))],sort(unique(gp)))
tiplabels(pie=model.matrix(~gp-1),piecol=cols,cex=0.3)
add.simmap.legend(colors=cols,prompt=FALSE,x=0.8*par()$usr[1],
y=-max(nodeHeights(tree)-.6),fsize=0.8)
Here we see our phylogeny with the continuous body size character (SVL) mapped to the branches, and displayed as a heat map. We also have several groups that define habitat use (ecomorphs). These are displayed at the tips of the phylogeny.
BM1: Brownian motion (one group)
A useful null model of evolutionary change is Brownian motion. For continuous characters, this is considered a neutral model, in the sense that there is no selection included: it simply represents random perturbations of the trait over time. This is akin to what might be expected under random genetic drift. The model assumes that changes are independent from time step to time step, so variance increases over time, but the mean does not change. An example is below:
nsim <- 100
t <- 0:100 # time
sig2 <- 0.005
X <- matrix(rnorm(n = nsim * (length(t) - 1), sd = sqrt(sig2)), nsim, length(t) - 1)
X <- cbind(rep(0, nsim), t(apply(X, 1, cumsum)))
plot(t, X[1, ], xlab = "time", ylab = "phenotype", ylim = c(-2, 2), type = "l")
apply(X[2:nsim, ], 1, function(x, t) lines(t, x), t = t)
## NULLMany phylogenetic comparative methods (PCMs) for continuous trait data have Brownian motion as their underlying model. For instance, methods for phylogenetic regression (phylogenetically independent contrasts, phylogenetic generalized least squares, phylogenetic transform), and methods for ancestral state estimation of continuous data (maximum likelihood, squared change parsimony) are all based on the Brownian model.
To fit this model and many others in R we can use the fitContinuous() function.
fitContinuous()
This function fits various models of continuous character evolution to trees. This function has many parameters but we’ll be using the following ones:
- \(phy\) a phylo object
- \(dat\) data vector for a single trait, with names matching tips in phy
- \(model\) A string with the type of model we will attempt to fit. We will be using the following models:
"BM"Brownian motion"OU"Ornstein-Uhlenbeck"EB"Early burst"lambda"Lambda"kappa"Early Burst
- \(bounds\) These place bounds on relevant parameters. This is particularly useful on likelihood surfaces that are long flat ridges, making it difficult for optimization. See examples in the tutorial for usage of \(bounds\)
This function returns the fit model with parameter estimates and summary statistics. Some elements of interest are:
- \(\$lik\) A likelihood function used to compute the likelihood of parameter values for the model
- \(\$opt\) This contains the results of optimization, namely: the resulting parameter estimates, maximum likelihood estimate, AIC, and sample-size corrected AIC (AICC)
fit.BM1<-fitContinuous(tree, svl, model="BM") #Brownian motion model
fit.BM1## GEIGER-fitted comparative model of continuous data
## fitted 'BM' model parameters:
## sigsq = 0.018223
## z0 = 4.053507
##
## model summary:
## log-likelihood = 5.256010
## AIC = -6.512019
## AICc = -6.360121
## free parameters = 2
##
## Convergence diagnostics:
## optimization iterations = 100
## failed iterations = 0
## number of iterations with same best fit = 100
## frequency of best fit = 1.00
##
## object summary:
## 'lik' -- likelihood function
## 'bnd' -- bounds for likelihood search
## 'res' -- optimization iteration summary
## 'opt' -- maximum likelihood parameter estimatesNotice that fitting the data to the phylogeny under the BM1 model provides the likelihood of that fit, its corresponding AIC (and AICc), as well as the parameters of the Brownian motion model: the rate parameter (sigma2), which describes the pace of evolutionary change (how fast the trait evolves along the phylogeny), and the phylogenetic mean (the estimate at the root of the phylogeny).
OU1: Ornstein-Uhlenbeck Model with a single optimum
Many times, a more realistic model is one containing both selection and drift. This may be found in an OU model. Here, the drift parameter remains, but there is also a selective ‘pull’ towards an optimum. The further the trait value is from the optimum, the stronger the pull. One can view the consequences of this model with the following simulations: first where the starting value and the optimum are the same (demonstrates the constraining force of the OU), and second when the starting value and optimum are not the same:
#from Lars Shmitz tutorial: http://schmitzlab.info/BMandOU.html
OU.sim <- function(n, theta, alpha, sigma,x0){
dw <- rnorm(n, 0)
x <- c(x0)
for (i in 2:(n+1)) {
x[i] <- x[i-1] + alpha*(theta-x[i-1]) + sigma*dw[i-1]
}
return(x);
}
OU1.sim1 <- replicate(100, OU.sim(n=100, theta=0.5, alpha=0.5, sigma=0.03, x0=0.5), simplify=FALSE) #No change
OU1.sim2 <- replicate(100, OU.sim(n=100, theta=0.75, alpha=0.5, sigma=0.03, x0=0.25), simplify=FALSE) #start and theta differ
par(mfcol = c(1, 2))
plot(OU1.sim1[[1]],xlab = "time", ylab = "phenotype", ylim = c(0,1), type = "l")
for(i in 2:100){lines(OU1.sim1[[i]])}
plot(OU1.sim2[[1]],xlab = "time", ylab = "phenotype", ylim = c(0,1), type = "l")
for(i in 2:100){lines(OU1.sim2[[i]])}
par(mfcol = c(1,1))Notice the impotant components of the OU model as illustrated here. First, unlike BM where variation increases proportional to time, the OU model displays ‘constrained’ variation over time. This is because of the ‘pull’ of selection to the optimum. Second, if the starting value differs from the optimum, there is a shift of the mean over time. This further illustrates the power of selection in this model. To fit an OU1 model for our data, we do the following:
options(warn=-1)
fit.OU1<-fitContinuous(tree, svl,model="OU") #OU1 model
fit.OU1## GEIGER-fitted comparative model of continuous data
## fitted 'OU' model parameters:
## alpha = 0.000000
## sigsq = 0.018227
## z0 = 4.053503
##
## model summary:
## log-likelihood = 5.258507
## AIC = -4.517014
## AICc = -4.209321
## free parameters = 3
##
## Convergence diagnostics:
## optimization iterations = 100
## failed iterations = 44
## number of iterations with same best fit = NA
## frequency of best fit = NA
##
## object summary:
## 'lik' -- likelihood function
## 'bnd' -- bounds for likelihood search
## 'res' -- optimization iteration summary
## 'opt' -- maximum likelihood parameter estimatesAs before, the logL, AIC, and parameter estimates are returned. Armed with these, one may compare the OU1 fit to the null-model BM1. This is accomplished using either using a likelihood ratio test (LRT), or by comparing AIC values (recall a difference > 4 is usually treated as supporting the alternative model).
LRT<- -(2*(fit.BM1$opt$lnL-fit.OU1$opt$lnL))
prob<-pchisq(LRT, 1, lower.tail=FALSE)
LRT## [1] 0.004994324prob## [1] 0.94366fit.BM1$opt$aic## [1] -6.512019fit.OU1$opt$aic## [1] -4.517014Here we find that under the OU1 model, the fit of the data to the phylogeny is not significantly improved. Thus we prefer the simpler model.
Other models
There are many alternative models that can be envisioned; many of which are easily implemented. In the category of ‘single group’ models, several common ones include:
- Brownian motion with a directional trend
- Early-burst models (rapid trait diversification early in clade history)
- Lambda models (covaration of phylogeny and trait change)
- Kappa models (trait change is puncutational and tends to occur at speciation events)
Models 2,3, and 4 are accomplished by adjusting the branch lengths of the phylogeny in some manner. As before, models are fit and then compared using AIC.
fit.BMtrend<-fitContinuous(tree, svl, model="trend") #Brownian motion with a trend
fit.EB<-fitContinuous(tree, svl, model="EB") #Early-burst model
fit.lambda<-fitContinuous(tree, svl,
bounds = list(lambda = c(min = exp(-5), max = 2)), model="lambda") #Lambda model
options(warn=-1)
fit.K<-fitContinuous(tree, svl, model="kappa") #Early-burst model
#Examine AIC
c(fit.BM1$opt$aic,fit.BMtrend$opt$aic,fit.EB$opt$aic,fit.lambda$opt$aic,fit.OU1$opt$aic)## [1] -6.512019 -6.981431 -7.235181 -5.517514 -4.517014 #NOTE: none of these generate dAIC > 4. So go with simplest model (BM1)OUM: Multi-group OU model
The above evolutionary models all assume that the taxa belong to a single ‘group.’ However, there may be more than one evolutionary optimum, and species may be evolving towards these optima. This requires and OUM approach (multi-group OU model). By convention, groups are referred to as regimes. The first step is to assign branches to regimes; commonly known as ‘painting’ the regimes on the tree. This may be done using some explicit biological hypothesis (e.g., ancestral branches are assumed to be one phenotype, while extant taxa belong to one or more groups). Alternatively, one may use ancestral state estimation to assign ancestral regimes, and then assign branches based on these and some rule of state change. Several examples of both are in Butler & King (2004).
Another approach is to use stochastic mapping, which is used here. We can fit multi-group OU models using the OUwie() function.
OUwie()
This function fits various OU models to trees. We’ll primarily be using this function to model regime shifts, each with their own OU parameters and dynamics, as such we won’t go over all the parameters but only the ones relevant for our goals
- \(phy\) This should be a
simmapobject, the output frommake.simmap(). See the ancestral state estimation tutorial for guidance on how to make a simmap. Essentially, it simulates some changing regime on the phylogeny. - \(data\) A data frame with trait and regime information. The trait data.frame must have column entries in the following order: [,1] species names, [,2] current selective regime, and [,3] the continuous trait of interest.
- \(model\) A string with the type of model we will attempt to fit. We will be using the following models:
"BM1"Brownian motion"BMS"Brownian motion with each regime having a different rate parameter"OU1"Ornstein-Uhlenbeck"OUM"Ornstein-Uhlenbeck with each regime having a different optimum"OUMV"Ornstein-Uhlenbeck with each regime having a different optimum and different rate
This function returns the fit model with parameter estimates and summary statistics. Some elements of interest are:
- \(\$loglik\) A likelihood function used to compute the likelihood of parameter values for the model
- \(\$AIC\) Akaike information criterion.
- \(\$AICc\) sample-size corrected AIC
- \(\$solution\) Maximum likelihood estimates for the parameters
data<-data.frame(Genus_species=names(svl),Reg=gp,X=svl) #input data.frame for OUwie
tree.simmap<-make.simmap(tree,gp) # perform & plot stochastic maps (we would normally do this x100)## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## CG GB TC TG Tr Tw
## CG -0.05727128 0.02402041 0.03325087 0.00000000 0.00000000 0.00000000
## GB 0.02402041 -0.20033484 0.04238120 0.06009264 0.01341265 0.06042794
## TC 0.03325087 0.04238120 -0.13585957 0.00000000 0.02427593 0.03595157
## TG 0.00000000 0.06009264 0.00000000 -0.06009264 0.00000000 0.00000000
## Tr 0.00000000 0.01341265 0.02427593 0.00000000 -0.03768858 0.00000000
## Tw 0.00000000 0.06042794 0.03595157 0.00000000 0.00000000 -0.09637951
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## CG GB TC TG Tr Tw
## 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667## Done.plot(tree.simmap,type='fan')## no colors provided. using the following legend:
## CG GB TC TG Tr Tw
## "black" "#DF536B" "#61D04F" "#2297E6" "#28E2E5" "#CD0BBC"
fitBM1<-OUwie(tree.simmap,data,model="BM1",simmap.tree=TRUE)## Initializing...
## Finished. Begin thorough search...
## Finished. Summarizing results.fitOU1<-OUwie(tree.simmap,data,model="OU1",simmap.tree=TRUE)## Initializing...
## Finished. Begin thorough search...
## Finished. Summarizing results.fitOUM<-OUwie(tree.simmap,data,model="OUM",simmap.tree=TRUE)## Initializing...
## Finished. Begin thorough search...
## Finished. Summarizing results.fitBM1##
## Fit
## lnL AIC AICc BIC model ntax
## 5.256014 -6.512028 -6.36013 -1.69859 BM1 82
##
## Rates
## alpha sigma.sq
## NA 0.01822168
##
## Optima
## 1
## estimate 4.0535069
## se 0.1009994
##
## Arrived at a reliable solutionfitOU1##
## Fit
## lnL AIC AICc BIC model ntax
## 5.25601 -4.512021 -4.204328 2.708137 OU1 82
##
##
## Rates
## alpha sigma.sq
## 0.000000001 0.018223440
##
## Optima
## 1
## estimate 4.0535069
## se 0.1010043
##
##
## Half life (another way of reporting alpha)
## alpha
## 693147181
##
## Arrived at a reliable solutionfitOUM #OUM is strongly preferred (examine AIC)##
## Fit
## lnL AIC AICc BIC model ntax
## 39.36652 -62.73304 -60.76044 -43.47929 OUM 82
##
##
## Rates
## CG GB TC TG Tr Tw
## alpha 0.27521974 0.27521974 0.27521974 0.27521974 0.27521974 0.27521974
## sigma.sq 0.01623802 0.01623802 0.01623802 0.01623802 0.01623802 0.01623802
##
## Optima
## CG GB TC TG Tr Tw
## estimate 5.8659001 3.73723424 4.2389072 4.25088507 3.7734415 4.0750891
## se 0.1511051 0.05473244 0.1025593 0.08680096 0.1952303 0.1157977
##
##
## Half life (another way of reporting alpha)
## CG GB TC TG Tr Tw
## 2.518523 2.518523 2.518523 2.518523 2.518523 2.518523
##
## Arrived at a reliable solution#How it SHOULD be run
#trees.simmap<-make.simmap(tree = tree,x = gp,nsim = 100) # 100 simmaps
#fitOUM.100<-lapply(1:100, function(j) OUwie(trees.simmap[[j]],data,model="OUM",simmap.tree=TRUE))
#OUM.AIC<-unlist(lapply(1:100, function(j) fitOUM.100[[j]]$AIC))
#hist(OUM.AIC)
#abline(v=fitBM1$AIC, lwd=2) #add value for BM1Try it out! OUwie over many ancestral state histories
Remember that the ancestral state estimation we maded with the make.simmap() function only recreates one possible state history given our discrete trait evolution model when in reality there are infinitely many histories that could generate our data! Here you will try to fit the various models in OUwie with different simmaps to see whether our results are robust to different ancestral state estimates.
- Create 100 simmaps with ancestral state histories using the
make.simmap()function. NOTE you can do this either in aforloop or by modifying thensimargument in themake.simmap()function itself. - For each simmap:
- Use the
OUwie()function to fit the ‘BM1’, ‘OU1’, and ’OUM` models to the simmap - Record the AIC for each fit model
- Make histograms of the AIC for each model. Which model on average fits the best over all the simmaps? Which model is the most robust to different evolutionary histories (where do we see the smallest variation in AIC)?
BMM: Multi-group BM model
Alternatively, there may be multiple rates across the phylogeny. For example, one may wish to test the hypothesis that trait evolution is faster in one group as compared to another (e.g., Do traits evolve faster in island taxa than mainland taxa?). Again groups are assigned, BM1 model and BMM models are fit to the data, and compared using LRT and AIC.
As was done for OU models, we can fit the multi-BM model with the OUwie() function by using the argument model="BMS" (in this case, BMS refers to Brownian Motion for each State).
# 3A: BM1 vs BMM: Comparing evolutionary rates
tree.simmap<-make.simmap(tree,gp) # perform & plot stochastic maps (we would normally do this x100)## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## CG GB TC TG Tr Tw
## CG -0.05727128 0.02402041 0.03325087 0.00000000 0.00000000 0.00000000
## GB 0.02402041 -0.20033484 0.04238120 0.06009264 0.01341265 0.06042794
## TC 0.03325087 0.04238120 -0.13585957 0.00000000 0.02427593 0.03595157
## TG 0.00000000 0.06009264 0.00000000 -0.06009264 0.00000000 0.00000000
## Tr 0.00000000 0.01341265 0.02427593 0.00000000 -0.03768858 0.00000000
## Tw 0.00000000 0.06042794 0.03595157 0.00000000 0.00000000 -0.09637951
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## CG GB TC TG Tr Tw
## 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667## Done.BMM.res<-OUwie(tree.simmap,data,model="BMS",simmap.tree = TRUE)## Initializing...
## Finished. Begin thorough search...
## Finished. Summarizing results.BMM.res##
## Fit
## lnL AIC AICc BIC model ntax
## 27.10382 -40.20764 -38.69413 -23.36061 BMS 82
##
## Rates
## CG GB TC TG Tr Tw
## alpha NA NA NA NA NA NA
## sigma.sq 0.1266219 0.01296043 0.02175058 0.005944775 0.001073071 0.001149187
##
## Optima
## CG GB TC TG Tr Tw
## estimate 3.92405727 3.92405727 3.92405727 3.92405727 3.92405727 3.92405727
## se 0.07896232 0.07896232 0.07896232 0.07896232 0.07896232 0.07896232
##
## Arrived at a reliable solutionExploratory Models: Identifying Rate Shifts
The above multi-group methods (OUM and BMM) are * a priori* methods, in that groups are pre-defined by the evolutionary biologist, and hypotheses tested based on those pre-specified groups. An alternative is to mine the data to identify regions on the phylogeny where putative rate shifts are most evident. Two MCMC Bayesian approaches have been proposed for this (Revell et al. 2012: Eastman et al. 2011). These methods are highly exploratory, and are found in the functions ‘evol.rate.mcmc’ in phytools, and the function ‘rjmcmc.bm’ in geiger .
# 3B: Identifying rate shifts on phylogeny
#Bayesian MCMC: single rate shift (Revell et al. 2012: Evol.)
BM.MCMC<-evol.rate.mcmc(tree=tree,x=svl, quiet=TRUE)## Starting MCMC run....## Done MCMC run.post.splits<-minSplit(tree,BM.MCMC$mcmc) #summarize
MCMC.post<-posterior.evolrate(tree=tree,mcmc=BM.MCMC$mcmc,tips = BM.MCMC$tips,ave.shift = post.splits)
#Plot rescaled to rates
ave.rates(tree,post.splits,extract.clade(tree, node=post.splits$node)$tip.label,
colMeans(MCMC.post)["sig1"], colMeans(MCMC.post)["sig2"],post.splits)
## $sig1
## [1] 0.01755507
##
## $sig2
## [1] 0.04166577 #Reversible-jump MCMC: multiple rate shifts (Eastman et al. 2011: Evol.)
BM.RJMC<-rjmcmc.bm(phy = tree,dat = svl)## | GENERATIONS COMPLETE |
## . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
##
## SAMPLING SUMMARY
## proposed accepted adoptrate
## mergesplit 9765 2 0.0002
## rootstate 1226 199 0.1623
## ratetune 95 18 0.1895
## moveshift 87 33 0.3793
## ratescale 8950 4529 0.5060
## movejump 9691 2469 0.2548
## incjump 4899 1641 0.3350
## decjump 4777 1636 0.3425
## jumpvar 8697 6150 0.7071
## SE 1813 639 0.3525#Run again for plotting
r <- paste(sample(letters,9,replace=TRUE),collapse="")
rjmcmc.bm(phy=tree, dat=svl, prop.width=1.5, ngen=20000, samp=500, filebase=r, simple.start=TRUE, type="rbm")## | GENERATIONS COMPLETE |
## . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
##
## SAMPLING SUMMARY
## proposed accepted adoptrate
## mergesplit 5572 15 0.0027
## rootstate 701 107 0.1526
## ratetune 5655 807 0.1427
## moveshift 4683 1491 0.3184
## ratescale 2351 784 0.3335
## movejump 0 0 0.0000
## incjump 0 0 0.0000
## decjump 0 0 0.0000
## jumpvar 0 0 0.0000
## SE 1038 347 0.3343## $method
## [1] "direct"
##
## $rate.lim
## $rate.lim$min
## [1] 0
##
## $rate.lim$max
## [1] Inf
##
##
## $root.lim
## $root.lim$min
## [1] -1000
##
## $root.lim$max
## [1] 1000
##
##
## $se.lim
## $se.lim$min
## [1] 0
##
## $se.lim$max
## [1] Inf
##
##
## $constrainJUMP
## [1] 0
##
## $dlnSHIFT
## function (x)
## sapply(x, function(xi) if (xi %in% yc) y[which(yc == xi)] else NA)
## <bytecode: 0x00000000326a9fb8>
## <environment: 0x000000002cfa2f40>
##
## $dlnJUMP
## function (x)
## sapply(x, function(xi) if (xi %in% yc) y[which(yc == xi)] else NA)
## <bytecode: 0x00000000326a9fb8>
## <environment: 0x000000002cf9dc98>
##
## $dlnROOT
## function (x)
## dunif(x, min = con$root.lim$min, max = con$root.lim$max, log = TRUE)
## <bytecode: 0x0000000033b7e2e8>
## <environment: 0x000000002d56a718>
##
## $dlnRATE
## function (x)
## dexp(x, rate = 1/(10^3), log = TRUE)
## <bytecode: 0x0000000033b79748>
## <environment: 0x000000002d56a718>
##
## $dlnPULS
## function (x)
## dexp(x, rate = 1/(10^3), log = TRUE)
## <bytecode: 0x0000000033b79748>
## <environment: 0x000000002d56a718>
##
## $dlnSE
## function (x)
## dexp(x, rate = 1/(10^3), log = TRUE)
## <bytecode: 0x0000000033b79748>
## <environment: 0x000000002d56a718>
##
## $jump.lim
## [1] 1
##
## $excludeSHIFT
## NULL
##
## $excludeJUMP
## NULL
##
## $bm.jump
## [1] 0.5
##
## $mergesplit.shift
## [1] 0.5
##
## $tune.scale
## [1] 0.65
##
## $slide.mult
## [1] 0.25
##
## $prob.dimension
## [1] 0.65
##
## $prob.effect
## [1] 0.3
##
## $prob.root
## [1] 0.02
##
## $prob.SE
## [1] 0.03
##
## $prop.width
## [1] 1.5
##
## $sample.priors
## [1] FALSE
##
## $simple.start
## [1] TRUE
##
## $summary
## [1] TRUE
##
## $beta
## [1] 1
##
## $filebase
## [1] "smcuqybol"
##
## $primary.parameter
## [1] "shifts"
##
## $lik
## likelihood function for relaxed-rates univariate continuous trait evolution
## argument names: rates, root, SE
##
## definition:
## function (rates, root, SE)
## {
## check.argn(rates, root)
## vv[mm] = rates
## datc_se = datc_init
## var[which(adjvar == 1)] = SE^2
## datc_se$var = var
## datc_se$y[rootidx] = root
## ll = ll.bm.direct(pars = vv, datc_se)
## return(ll)
## }
## <bytecode: 0x00000000301e0718>
## <environment: 0x000000002d31c5b8>
##
##
## NOTE: 'rates' vector should be given in order given by the function 'argn()'
## $model
## [1] "relaxedBM"
##
## $algo
## [1] "rjmcmc"
##
## $prop.cs
## dim effect root se
## 0.65 0.95 0.97 1.00
##
## $n.subprop
## mergesplit rootstate ratetune moveshift ratescale movejump incjump
## 5572 701 5655 4683 2351 0 0
## decjump jumpvar SE
## 0 0 1038
##
## $n.subaccept
## mergesplit rootstate ratetune moveshift ratescale movejump incjump
## 15 107 807 1491 784 0 0
## decjump jumpvar SE
## 0 0 347
##
## $parmbase
## [1] "relaxedBM.smcuqybol"
##
## $parlogger
## function (delta, values = NULL)
## {
## xx = which(delta != 0)
## nd = c(nm[rr], nm[xx])
## if (length(delta) != length(nm))
## stop(paste("Expecting 'delta' of length ", length(nm),
## ".", sep = ""))
## if (is.null(values)) {
## values = delta
## }
## else {
## if (!length(values) == length(delta))
## stop("'delta' and 'values' are mismatched.")
## }
## vv = values[c(rr, xx)]
## ests = paste(nd, vv, collapse = "\t", sep = "\t")
## ests
## }
## <bytecode: 0x000000002e276748>
## <environment: 0x000000002ccb48a8>
##
## $thin
## [1] 500
##
## $ngen
## [1] 20000outdir <- paste("relaxedBM", r, sep=".")
ps <- load.rjmcmc(outdir)
plot(x=ps, par="shifts", burnin=0.25, legend=TRUE, show.tip=FALSE, edge.width=2,type='fan')
## branch shift.direction shift.probability
## 9e9e6f56c19c8c9f957e6077d58fd6cf 84 0 0.00000000
## ba7ea74079079929e17d82684c2b62bd 85 0 0.00000000
## 2ce4bf143fed53bf2b155b0eef7fc0e9 86 0 0.00000000
## 004ac217dd86ef893f9382d3a787eaef 87 0 0.00000000
## 97a9968a581368a986b7c7bad9078fdb 88 0 0.00000000
## f4fef45845bc1d781dcfbfa6b69831ad 89 0 0.00000000
## 603900d7d6c439214c5d6acd83da3d5e 90 0 0.00000000
## 95c2de48b98c17822b5841935c358b76 1 0 0.00000000
## bf8b62ed9e8b10245aa2c355ae08acbf 2 0 0.00000000
## 8ea1316df44bea07ccd8aaf33affeab1 3 0 0.00000000
## 0973b4d527844bb03d7b022ddb210644 4 0 0.00000000
## 91d6c3151ae3eba740750b0975615625 91 0 0.00000000
## 569bc41ab2a37b45e036eb32f8cc726c 92 0 0.00000000
## 0cb7a72d69899366edd15a3d99305360 93 0 0.00000000
## 60f198711d22960e9b986eea14e05ba8 94 0 0.00000000
## e6f6e7b9725804e4b133bc439748a41c 5 0 0.00000000
## 8bf06f2925d0ab96fec9e5f76b623742 95 0 0.00000000
## dd4e455b4916fdf887994e6460209684 6 0 0.00000000
## 9a4c8e7a276ed5b93d359c9f05d6214d 7 0 0.00000000
## 2b8349af4fc2e0355c9aea7aec4ccbdb 8 0 0.00000000
## 656880f1afda9aa525e0676a01c931c5 9 0 0.00000000
## 6b0ef0ee0bdbb5a29cbf4e4ce1de8f05 96 0 0.00000000
## 264d974f4a6a232465fab628db06d2e3 97 0 0.00000000
## 431583761f26dc63ad8a65c3e1014066 98 0 0.00000000
## f4f409705cfbeebed10b25094c7a8bdd 10 0 0.00000000
## 6565d1578c75c75f961edceac12b120b 11 0 0.00000000
## af7240ca03d1fb0a8dd90279d3ff9fb1 12 0 0.00000000
## 576baf5405148ec072fb5cdf8380cea8 13 0 0.00000000
## 3da2aeeedcceade46ad08988e48ed2c2 99 0 0.00000000
## 6a63add632896bf7d3d528c14ba4d14c 100 0 0.00000000
## 25ca9224382b47dd42763ceb57254425 101 0 0.00000000
## 40d51a39dbc43b3ba80e91b14cf36d95 102 1 0.16666667
## b42807f5ee6d6b0b6a20f88fea16feab 14 1 0.83333333
## 36f45dc485575dd799011bbd3b32e63e 15 0 0.00000000
## b4965f960ee721b68ea56a26e9237d8d 16 0 0.00000000
## 047aa8e63e98f543272170eb0e036fb3 17 0 0.00000000
## 92ac903ba89ce37a6d0329d3399c7c36 18 0 0.00000000
## ff7637498521a016194681c647dbe9e9 103 0 0.00000000
## 9dfb9442f171550122a94f704e73db0f 104 0 0.00000000
## a07f786ba39612c433f50d15bfe20304 105 0 0.00000000
## 77852d59a89ff60fa6f97c90a9018ffe 19 0 0.00000000
## 87ae66a8bc2655f8f411eafe8e6920cb 20 0 0.00000000
## 02cabc58f2db0710ee1fd3f8773939c9 106 0 0.00000000
## 3fed3a12910ad5828665fcd22bcc4848 107 0 0.00000000
## 14ef3a67e511d3bbd57ee52a1ea04eac 108 0 0.00000000
## 2fe6778ac111070ee4bdc7d139295a23 21 0 0.00000000
## ba4832911e37f1cf41430a03c6ceef5f 22 0 0.00000000
## 57d80a343acec3e69bcfbc46dfe9341a 109 0 0.00000000
## 8d8042fc58f4d5c157df04d6c10e4e1c 23 0 0.00000000
## d3ecc9d765519ea7da9ba6813012d70e 24 0 0.00000000
## c04c191dbc1c8b9419997785a3bcce09 110 0 0.00000000
## 0a6a7dce1a1191a928e6f9bd089ea331 25 0 0.00000000
## bfd581f502d01b7dfd7410275b088617 26 0 0.00000000
## 67354635ec5ec66f778820b953cd1540 111 0 0.00000000
## 7054ebd2f005e3326fa63bb61e57b884 112 0 0.00000000
## 92a54765bdb3903efa5c5f464aa648ec 113 0 0.00000000
## 960fb6f7fbf37676c4def46f92e92561 27 0 0.00000000
## 42fa18e6a0afc3de663daee528721ed0 114 0 0.00000000
## 83b2a2a5a73e11f84f013b6316526c6c 28 0 0.00000000
## 87c6d1721e26e0eaf9c8ef863f3b529b 29 0 0.00000000
## 6d0ea960023b658217fc3ccff0919755 30 0 0.00000000
## fb006c4bc7a286192ca4d43a2bece288 31 0 0.00000000
## 52894f3b31688bc8eb9ab5d63a02604e 115 0 0.00000000
## 71945eb2c8d49781ea563f9a1ba914ab 116 0 0.00000000
## ed693e97b591ab5f6b25ebf9c01a18db 117 0 0.00000000
## c3ab933e0e1eee5d950706cbe9b3f8d6 118 0 0.00000000
## d6813367d75d6e444779a28d466b9a87 119 0 0.00000000
## f7c616f76557208f623500980d91d108 32 0 0.00000000
## 9916c107f2cd71d8463ac2d2f840c05b 33 0 0.00000000
## 441eab8d21a722fc1cc74f98631a24b4 34 0 0.00000000
## 5dd9c024ff99d34c5277810fd103d82e 35 0 0.00000000
## ff325966dd9d0f3a1c3821e05917a8d9 120 0 0.00000000
## 1ae0f500f6e841545b524cbc8d496ed5 121 0 0.00000000
## 6439edc5bd92a66335e99bafcf07ad8d 122 0 0.00000000
## 918187102d3365572c59d51db99a464f 123 0 0.00000000
## 755cc2aa986d15f548643ba66632f4db 36 0 0.00000000
## 02a0b4ad260b7c933c431e6b5b4e8d1d 124 0 0.00000000
## 5663a710ebdd340a53b3fe6a35d5e971 125 0 0.00000000
## 62f09823ef21ed57d125d676e105cb44 37 0 0.00000000
## 12c8a839ca57290722c587672467dfdb 38 0 0.00000000
## 499d1845410997d6c732766cb4c786ea 126 0 0.00000000
## e10c35a4e314b8b27811ff9b25f3a51a 127 0 0.00000000
## 022548499ae362530bf1c4b1a803daaf 39 0 0.00000000
## 78e954d36db5f5b6f9e2e93ff93ecec4 40 0 0.00000000
## 0ffa3d6a4883f1c6c91e714545aa1ea1 41 0 0.00000000
## 3e0d8cbab26f2ee684135baf25a27982 128 0 0.00000000
## be56f5bf258a5931dc7421de7737cf5a 42 0 0.00000000
## 9a6586721c1d94c0339c49b76bc40120 43 0 0.00000000
## 0c6fa5bd7e77ad5fa9d51925e2f4a70f 44 0 0.00000000
## 9902dbf078b25aefa62485146ab9a93d 129 1 0.03333333
## 89c2c8db45812f43b4d41ef5040a32ed 130 0 0.00000000
## 8e77892a4494a82cb88ef298cfffc43e 131 0 0.00000000
## 0b6a64daa4e227f5302536091da1b515 45 0 0.00000000
## 3e89682c9067b34c74ff50c4b88699f8 46 0 0.00000000
## e4b296be23ddbb232c40a7003faa0098 47 0 0.00000000
## 300e149ddd44d1673f5f93377b7ae5dc 48 0 0.00000000
## 8204cae4dcd001c2f77dcec55f4dd8f9 132 0 0.00000000
## b8ee22e4656bec7421347b732a9f29cf 133 0 0.00000000
## ec4e1a14d20b7306662a42c2daea951d 134 0 0.00000000
## a856a8dbb9f80621d0a42c31a3aabce4 135 0 0.00000000
## 883989cb504be968432fc341db70ddc0 136 0 0.00000000
## 69693ad89a162354c8a6a0cbe08de7e3 137 0 0.00000000
## 9181fd77007c6869eded8548480fa3a7 138 0 0.00000000
## c94fe834ab9033aec4f82a85d4940c5b 49 0 0.00000000
## f0bd6ff3ed42e4ebf2ee7c60ed9fad8d 50 0 0.00000000
## a7d0b17a669a714d98d4a5ffebda167b 51 0 0.00000000
## 5251efc6e47c651bcce94c9ea5fb0145 139 0 0.00000000
## c5cdb7b7591179ab2ef87f58ca03f02b 52 0 0.00000000
## 9e3f3fd9ccc0d37c151343dfb8c1798b 53 0 0.00000000
## ee28b018af326ec8f948f38867d379e9 54 0 0.00000000
## 3a54efeafddd88bf9a54e12d56cfb3d3 55 0 0.00000000
## c09ed9692bff14853c7301c113d71a4e 140 0 0.00000000
## b4150c746d3fe184e8b93940c2f2dede 141 0 0.00000000
## f988bb267cab0bc73e329c543e24eb38 56 0 0.00000000
## 6990f145a8d187917b15fd059a499555 57 0 0.00000000
## cadd0753cc08ded56c06480638173cb0 142 0 0.00000000
## bce1b5aa8956d47e911daecbcfe6fc55 58 0 0.00000000
## c1878160be7a5f962469d7878d8fe7ed 143 0 0.00000000
## d00bca2c5e7c6ff7e5c8e3c94b02bdd5 59 0 0.00000000
## dc874bb6a84c70fbd13729af46eccb2a 60 0 0.00000000
## d9a70b0089a2bc87e32f3b28b32806ab 144 0 0.00000000
## 929f5540b4c3a3c94d97b3efe3dfdf05 145 0 0.00000000
## f0616b5889ad6a3f494990fb5bb91bb6 61 0 0.00000000
## 5daf4de974c85a1c1978169726468e4a 62 0 0.00000000
## ec59fe99c41364fd3fa5f3137b5ad381 146 0 0.00000000
## 4fffc1c4e83bc3a691098f88cd616ba6 147 0 0.00000000
## ac57cf8d2c470c6c19087db803eae9c7 148 0 0.00000000
## 0a734a506eafab0ff4f8301f1f7ca37f 63 0 0.00000000
## 05f4d067dec4e85cd3324def026005ba 149 0 0.00000000
## 5c700795497368998fe7ed6afc744ed6 64 0 0.00000000
## ebb2580313c6233d02822a69e564306a 65 0 0.00000000
## 1e039da3165ec5f7d83699c11c0f83a5 150 0 0.00000000
## 41f18d58cb8448f005f578c8bb0c2a12 66 0 0.00000000
## 1baf6058dccd1dc09abd2a8b0ea6f17e 67 0 0.00000000
## 7a8cbde95131ebc383d2546d622901d3 68 0 0.00000000
## 29ac11334461bc059d5b0e177f5d5129 151 0 0.00000000
## cbe83cce090d5cec35cb4652229a5481 152 0 0.00000000
## 824893c7e70d576e32428790740336b2 153 0 0.00000000
## 17c2a21ccc426e7c2c9f90f7ae036f6f 154 0 0.00000000
## c6e8bf77689b2a281544a7de2b68b2b7 155 0 0.00000000
## a6efbc1bb4581c63e074939131a0dd57 69 0 0.00000000
## 58be1f0b2a7997ce5d1742574dc47793 156 0 0.00000000
## 42619a887689308e0f3b5a49a043dffe 70 0 0.00000000
## a46c03417feae16030e0fb82f82ca69a 71 0 0.00000000
## 2ae34efc4141fdacf26d52895aed2069 72 0 0.00000000
## a1595d98874906937cffac11a8b7366e 73 0 0.00000000
## b5c6696860f4debb2f2caef531c7770c 157 0 0.00000000
## 541c6d1b4906b2abce12fc3d65176783 158 0 0.00000000
## e67909d62f94990fcf77142dd561fbcb 159 0 0.00000000
## 91f3d57e75e7dce641d65b47ca76ea76 74 0 0.00000000
## 564b5b1be4dcf93ce611d8824d1ecb40 160 0 0.00000000
## 6f0955bd9c4b2f4924dbf077396b9e6b 75 0 0.00000000
## 68c8b67f4d8111eb1af13a412c9d7819 76 0 0.00000000
## 4017179115849ff0f3aaa4b44d26eb27 77 0 0.00000000
## 3953908561d77fbdba3d5d11b7d6bd41 161 0 0.00000000
## be2a2b5fb1a367178f785ce6d62cc639 162 0 0.00000000
## 18cc77c7e638fecfe5eb7c7023fe25c6 163 0 0.00000000
## 3ce7e59f3d75fbeddcf39402c44549b6 78 0 0.00000000
## 500c38508cb0d393abd472c506cbe2f9 79 0 0.00000000
## 23b505862258a1437c0e059290ca2c6d 80 0 0.00000000
## e4e19e8484f31d2449b8cc65e1e47627 81 0 0.00000000
## 59537e87050d7ef8f71d5a5fb742a63a 82 0 0.00000000
## median.rates
## 9e9e6f56c19c8c9f957e6077d58fd6cf 0.01485279
## ba7ea74079079929e17d82684c2b62bd 0.01485279
## 2ce4bf143fed53bf2b155b0eef7fc0e9 0.01485279
## 004ac217dd86ef893f9382d3a787eaef 0.01485279
## 97a9968a581368a986b7c7bad9078fdb 0.01485279
## f4fef45845bc1d781dcfbfa6b69831ad 0.01485279
## 603900d7d6c439214c5d6acd83da3d5e 0.01485279
## 95c2de48b98c17822b5841935c358b76 0.01485279
## bf8b62ed9e8b10245aa2c355ae08acbf 0.01485279
## 8ea1316df44bea07ccd8aaf33affeab1 0.01485279
## 0973b4d527844bb03d7b022ddb210644 0.01485279
## 91d6c3151ae3eba740750b0975615625 0.01485279
## 569bc41ab2a37b45e036eb32f8cc726c 0.01485279
## 0cb7a72d69899366edd15a3d99305360 0.01485279
## 60f198711d22960e9b986eea14e05ba8 0.01485279
## e6f6e7b9725804e4b133bc439748a41c 0.01485279
## 8bf06f2925d0ab96fec9e5f76b623742 0.01485279
## dd4e455b4916fdf887994e6460209684 0.01485279
## 9a4c8e7a276ed5b93d359c9f05d6214d 0.01485279
## 2b8349af4fc2e0355c9aea7aec4ccbdb 0.01485279
## 656880f1afda9aa525e0676a01c931c5 0.01485279
## 6b0ef0ee0bdbb5a29cbf4e4ce1de8f05 0.01485279
## 264d974f4a6a232465fab628db06d2e3 0.01485279
## 431583761f26dc63ad8a65c3e1014066 0.01485279
## f4f409705cfbeebed10b25094c7a8bdd 0.01485279
## 6565d1578c75c75f961edceac12b120b 0.01485279
## af7240ca03d1fb0a8dd90279d3ff9fb1 0.01485279
## 576baf5405148ec072fb5cdf8380cea8 0.01485279
## 3da2aeeedcceade46ad08988e48ed2c2 0.01485279
## 6a63add632896bf7d3d528c14ba4d14c 0.01485279
## 25ca9224382b47dd42763ceb57254425 0.01485279
## 40d51a39dbc43b3ba80e91b14cf36d95 0.01514985
## b42807f5ee6d6b0b6a20f88fea16feab 12.58675528
## 36f45dc485575dd799011bbd3b32e63e 0.01514985
## b4965f960ee721b68ea56a26e9237d8d 0.01485279
## 047aa8e63e98f543272170eb0e036fb3 0.01485279
## 92ac903ba89ce37a6d0329d3399c7c36 0.01485279
## ff7637498521a016194681c647dbe9e9 0.01485279
## 9dfb9442f171550122a94f704e73db0f 0.01485279
## a07f786ba39612c433f50d15bfe20304 0.01485279
## 77852d59a89ff60fa6f97c90a9018ffe 0.01485279
## 87ae66a8bc2655f8f411eafe8e6920cb 0.01485279
## 02cabc58f2db0710ee1fd3f8773939c9 0.01485279
## 3fed3a12910ad5828665fcd22bcc4848 0.01485279
## 14ef3a67e511d3bbd57ee52a1ea04eac 0.01485279
## 2fe6778ac111070ee4bdc7d139295a23 0.01485279
## ba4832911e37f1cf41430a03c6ceef5f 0.01485279
## 57d80a343acec3e69bcfbc46dfe9341a 0.01485279
## 8d8042fc58f4d5c157df04d6c10e4e1c 0.01485279
## d3ecc9d765519ea7da9ba6813012d70e 0.01485279
## c04c191dbc1c8b9419997785a3bcce09 0.01485279
## 0a6a7dce1a1191a928e6f9bd089ea331 0.01485279
## bfd581f502d01b7dfd7410275b088617 0.01485279
## 67354635ec5ec66f778820b953cd1540 0.01485279
## 7054ebd2f005e3326fa63bb61e57b884 0.01485279
## 92a54765bdb3903efa5c5f464aa648ec 0.01485279
## 960fb6f7fbf37676c4def46f92e92561 0.01485279
## 42fa18e6a0afc3de663daee528721ed0 0.01485279
## 83b2a2a5a73e11f84f013b6316526c6c 0.01485279
## 87c6d1721e26e0eaf9c8ef863f3b529b 0.01485279
## 6d0ea960023b658217fc3ccff0919755 0.01485279
## fb006c4bc7a286192ca4d43a2bece288 0.01485279
## 52894f3b31688bc8eb9ab5d63a02604e 0.01485279
## 71945eb2c8d49781ea563f9a1ba914ab 0.01485279
## ed693e97b591ab5f6b25ebf9c01a18db 0.01485279
## c3ab933e0e1eee5d950706cbe9b3f8d6 0.01485279
## d6813367d75d6e444779a28d466b9a87 0.01485279
## f7c616f76557208f623500980d91d108 0.01485279
## 9916c107f2cd71d8463ac2d2f840c05b 0.01485279
## 441eab8d21a722fc1cc74f98631a24b4 0.01485279
## 5dd9c024ff99d34c5277810fd103d82e 0.01485279
## ff325966dd9d0f3a1c3821e05917a8d9 0.01485279
## 1ae0f500f6e841545b524cbc8d496ed5 0.01485279
## 6439edc5bd92a66335e99bafcf07ad8d 0.01485279
## 918187102d3365572c59d51db99a464f 0.01485279
## 755cc2aa986d15f548643ba66632f4db 0.01485279
## 02a0b4ad260b7c933c431e6b5b4e8d1d 0.01485279
## 5663a710ebdd340a53b3fe6a35d5e971 0.01485279
## 62f09823ef21ed57d125d676e105cb44 0.01485279
## 12c8a839ca57290722c587672467dfdb 0.01485279
## 499d1845410997d6c732766cb4c786ea 0.01485279
## e10c35a4e314b8b27811ff9b25f3a51a 0.01485279
## 022548499ae362530bf1c4b1a803daaf 0.01485279
## 78e954d36db5f5b6f9e2e93ff93ecec4 0.01485279
## 0ffa3d6a4883f1c6c91e714545aa1ea1 0.01485279
## 3e0d8cbab26f2ee684135baf25a27982 0.01485279
## be56f5bf258a5931dc7421de7737cf5a 0.01485279
## 9a6586721c1d94c0339c49b76bc40120 0.01485279
## 0c6fa5bd7e77ad5fa9d51925e2f4a70f 0.01485279
## 9902dbf078b25aefa62485146ab9a93d 0.01496415
## 89c2c8db45812f43b4d41ef5040a32ed 0.01496415
## 8e77892a4494a82cb88ef298cfffc43e 0.01496415
## 0b6a64daa4e227f5302536091da1b515 0.01496415
## 3e89682c9067b34c74ff50c4b88699f8 0.01496415
## e4b296be23ddbb232c40a7003faa0098 0.01496415
## 300e149ddd44d1673f5f93377b7ae5dc 0.01496415
## 8204cae4dcd001c2f77dcec55f4dd8f9 0.01485279
## b8ee22e4656bec7421347b732a9f29cf 0.01485279
## ec4e1a14d20b7306662a42c2daea951d 0.01485279
## a856a8dbb9f80621d0a42c31a3aabce4 0.01485279
## 883989cb504be968432fc341db70ddc0 0.01485279
## 69693ad89a162354c8a6a0cbe08de7e3 0.01485279
## 9181fd77007c6869eded8548480fa3a7 0.01485279
## c94fe834ab9033aec4f82a85d4940c5b 0.01485279
## f0bd6ff3ed42e4ebf2ee7c60ed9fad8d 0.01485279
## a7d0b17a669a714d98d4a5ffebda167b 0.01485279
## 5251efc6e47c651bcce94c9ea5fb0145 0.01485279
## c5cdb7b7591179ab2ef87f58ca03f02b 0.01485279
## 9e3f3fd9ccc0d37c151343dfb8c1798b 0.01485279
## ee28b018af326ec8f948f38867d379e9 0.01485279
## 3a54efeafddd88bf9a54e12d56cfb3d3 0.01485279
## c09ed9692bff14853c7301c113d71a4e 0.01485279
## b4150c746d3fe184e8b93940c2f2dede 0.01485279
## f988bb267cab0bc73e329c543e24eb38 0.01485279
## 6990f145a8d187917b15fd059a499555 0.01485279
## cadd0753cc08ded56c06480638173cb0 0.01485279
## bce1b5aa8956d47e911daecbcfe6fc55 0.01485279
## c1878160be7a5f962469d7878d8fe7ed 0.01485279
## d00bca2c5e7c6ff7e5c8e3c94b02bdd5 0.01485279
## dc874bb6a84c70fbd13729af46eccb2a 0.01485279
## d9a70b0089a2bc87e32f3b28b32806ab 0.01485279
## 929f5540b4c3a3c94d97b3efe3dfdf05 0.01485279
## f0616b5889ad6a3f494990fb5bb91bb6 0.01485279
## 5daf4de974c85a1c1978169726468e4a 0.01485279
## ec59fe99c41364fd3fa5f3137b5ad381 0.01485279
## 4fffc1c4e83bc3a691098f88cd616ba6 0.01485279
## ac57cf8d2c470c6c19087db803eae9c7 0.01485279
## 0a734a506eafab0ff4f8301f1f7ca37f 0.01485279
## 05f4d067dec4e85cd3324def026005ba 0.01485279
## 5c700795497368998fe7ed6afc744ed6 0.01485279
## ebb2580313c6233d02822a69e564306a 0.01485279
## 1e039da3165ec5f7d83699c11c0f83a5 0.01485279
## 41f18d58cb8448f005f578c8bb0c2a12 0.01485279
## 1baf6058dccd1dc09abd2a8b0ea6f17e 0.01485279
## 7a8cbde95131ebc383d2546d622901d3 0.01485279
## 29ac11334461bc059d5b0e177f5d5129 0.01485279
## cbe83cce090d5cec35cb4652229a5481 0.01485279
## 824893c7e70d576e32428790740336b2 0.01485279
## 17c2a21ccc426e7c2c9f90f7ae036f6f 0.01485279
## c6e8bf77689b2a281544a7de2b68b2b7 0.01485279
## a6efbc1bb4581c63e074939131a0dd57 0.01485279
## 58be1f0b2a7997ce5d1742574dc47793 0.01485279
## 42619a887689308e0f3b5a49a043dffe 0.01485279
## a46c03417feae16030e0fb82f82ca69a 0.01485279
## 2ae34efc4141fdacf26d52895aed2069 0.01485279
## a1595d98874906937cffac11a8b7366e 0.01485279
## b5c6696860f4debb2f2caef531c7770c 0.01485279
## 541c6d1b4906b2abce12fc3d65176783 0.01485279
## e67909d62f94990fcf77142dd561fbcb 0.01485279
## 91f3d57e75e7dce641d65b47ca76ea76 0.01485279
## 564b5b1be4dcf93ce611d8824d1ecb40 0.01485279
## 6f0955bd9c4b2f4924dbf077396b9e6b 0.01485279
## 68c8b67f4d8111eb1af13a412c9d7819 0.01485279
## 4017179115849ff0f3aaa4b44d26eb27 0.01485279
## 3953908561d77fbdba3d5d11b7d6bd41 0.01485279
## be2a2b5fb1a367178f785ce6d62cc639 0.01485279
## 18cc77c7e638fecfe5eb7c7023fe25c6 0.01485279
## 3ce7e59f3d75fbeddcf39402c44549b6 0.01485279
## 500c38508cb0d393abd472c506cbe2f9 0.01485279
## 23b505862258a1437c0e059290ca2c6d 0.01485279
## e4e19e8484f31d2449b8cc65e1e47627 0.01485279
## 59537e87050d7ef8f71d5a5fb742a63a 0.01485279Skull Thickness of Fizzgigs
Fizzgigs are little balls of fur that live in the Endless Forest of the Dark Crystal world. They have four stubby paws that they use to travel short distances but will roll up into a ball for longer journeys They have sharp teeth used for catching and gripping prey and even have a second set of teeth near their uvula. It is said that they have the largest bite radii of the animals in Thra. There is much diversity among the various species of Fizzgigs when considering their bite radii but not much is known how this diversity came to be.
Try it out! Evolution of biting Fizzgigs
We’re going to try and fit various evolutionary models to a phylogenetic tree of Fizzgigs to see how their bite radii might have evolved
- Use
read.tree()to read in the tree titledfizzgig.tre. This file can be found in the data folder of this tutorial: - Use the
read.csv()function to read in the data set calledfizzgig_data.csv. NOTE there are two columns in this data set, one labeledbitethat contains the bit radii and another labeledtype. For now we’ll only need thebiteinformation - use
fitContinuous()to fit various evolutionary models and compare their AICs. Which model fit best? What are the maximum likelihood estimates for that model?
Now, let us suppose some of these Fizzgigs had different diets and that some are insectivores and others are avivores (eats birds). If we believed that these diets affected their bite radii evolution, how would we model it?
- Use the
typeinformation of thefizzgig_data.csvto make a simmap of how their diets may have evolved over time. - Use
OUwie()along with the simmap to test out various multiple-state models. Which fits the best and what might this say about their diets?