Analysis of locomotion and forelimb morphology in carnivorans using the mvSLOUCH R-package

Jesualdo Fuentes Gonzalez, Jason Pienaar and Krzysztof Bartoszek

2023-11-20

Introduction

We will explore several attributes of the R (R Core Team 2020) package mvSLOUCH (Bartoszek et al. 2012), that allows fitting multivariate models for phylogenetic comparative data with emphasis on those based on an Ornstein-Uhlenbeck process. Versions of mvSLOUCH from 2.0.0 run models at considerably faster speeds through using the computational engine provided by PCMBase (Mitov et al. 2020), so let us start by attaching mvSLOUCH as well as
ggplot2 (Wickham 2016), which have some useful functions for data processing and plotting, respectively (PCMBase with its suggested PCMBaseCpp, if installed, C++ backend will be loaded by mvSLOUCH):

library(ggplot2)
library(ape)
library(mvSLOUCH)
## Loading required package: abind

As a phylogenetic comparative analysis with mvSLOUCH can run a long time, we first load required fragments of the precalculated objects. The unreduced objects are the results of the estimation procedures run here with the set random seed. We preload them so that building the vignette does not take over a day and we store within mvSLOUCH only the necessary parts as the full objects would take up too much space. The complete objects can be found in the KJVJMRKK_mvSLOUCH GitHub repository. However, the readers are encouraged to run the presented code by themselves.

load("./Carnivora_mvSLOUCH_objects.RData")

Using an example on carnivoran locomotion and forelimb morphology, we will explore and compare the basic models of mvSLOUCH, going through some key tasks associated with its inputs (e.g. setting up the data and the phylogeny, specifying selective regimes for adaptive hypotheses) and outputs (e.g. identifying key statistics, optimizing parameter estimates of a preferred model, computing confidence intervals under parametric bootstrap).

The function set.seed() allows specifying a starting point in a sequence of randomly generated numbers so that a user can obtain the same outputs under a given process. For the purposes of the current vignette, if you want to replicate the outputs below (for mvSLOUCH 2.6.2), set up the following seed:

RNGversion(min(as.character(getRversion()),"3.6.1"))
set.seed(5, kind = "Mersenne-Twister", normal.kind = "Inversion")

Data: Locomotion and forelimb morphology in carnivorans

The order Carnivora has colonized a wide variety of habitats. The specific challenges of moving through these habitats are reflected in the diversity of locomotor strategies they exhibit, such as running fast or for long distances (i.e. cursorial locomotion, such as hyenas and wolves), climbing (e.g. from the scansorial raccoon to the fully arboreal kinkajou) swimming (e.g. from the semiaquatic otter to the fully aquatic seal) and digging (i.e. semifossorial locomotion, such as some skunks and badgers). Although some morphological attributes are useful across different locomotor types (e.g. swimmers and diggers are similarly benefited by an increased area of the paws and high force outputs of the limbs), others can be at odds with each other. A good example of the latter involves cursorial and semifossorial carnivorans. At the core of the contrast between these two locomotor types lies a trade-off in limb mechanics, as many adaptations that maximize velocity transmissions are at odds with those that maximize force outputs. Limb bones selected for strength exhibit pronounced crests that increase the area of insertion for locomotor muscles, and limb segments tend to be short given that smaller output levers result in higher force outputs. On the other hand, elongated limb segments result in longer strides and larger output levers that favor runners (because larger output levers increase relative velocity transmissions). Runners also benefit from lighter limbs that maximize the distance gained per force input of each stride i.e. big muscles and conspicuous crests are less advantageous for runners as they are for diggers. The mvSLOUCH package offers analytical tools for evaluating evolutionary hypotheses that both:

  1. relate limb morphologies to locomotor types (e.g. runners vs diggers); and
  2. test for trade-offs between force outputs and velocity transmissions in limb mechanics (i.e. strength vs speed).

We will explore these ideas using a subset of the dataset collected by Samuels, Meachen, and Sakai (2013), available at the Dryad Data Repository dx.doi.org/10.5061/dryad.77tm4. We first download the data, remove fossil samples (lacking locomotor ecology data; Urocyon cinereoargenteus is removed separately as the current and fossil samples have the same entry as species name in the data file), species with missing measurements (for radius length, see below) and Lycalopex sp. (as species identification is required for branch length assignation), rename species according to Johnson et al. (2006) and Wilson and Reeder (2005), and match the tip labels in the phylogeny (spaces replaced by underscores), and finally keep only the columns that we need (locomotor habits, humerus length, deltopectoral crest length and radius length of the forelimb; see below).

b_correct_dryad_download<-FALSE
temp <- tempfile()
tryCatch({
download.file("datadryad.org/api/v2/datasets/doi%253A10.5061%252Fdryad.77tm4/download",temp)
b_correct_dryad_download<-TRUE
},error=function(e){cat("Could not download data file from Dryad! No analysis can be done! Vignette not built!")},
warning=function(w){cat("Problem with downloading data file from Dryad! No analysis can be done! Vignette not built!")})

if (b_correct_dryad_download){
    b_correct_dryad_download<-FALSE
    tryCatch({
    dfcarnivores_postcranial <- read.table(unz(temp, "Carnivore Postcranial Morphology Data Samuels et al. 2013.txt"),header=TRUE,sep="\t",stringsAsFactors =FALSE)
    b_correct_dryad_download<-TRUE
    },error=function(e){cat("Corrupt data file from Dryad! No analysis can be done! Vignette not built!")},
    warning=function(w){cat("Problem with accessing data file from Dryad! No analysis can be done! Vignette not built!")})
}
unlink(temp)
Urocyon_cinereoargenteus_duplicated<-which(dfcarnivores_postcranial[,2]=="Urocyon cinereoargenteus")[2]
dfcarnivores_postcranial<-dfcarnivores_postcranial[-Urocyon_cinereoargenteus_duplicated,]
v_species_to_remove<-c("Bdeogale crassicauda","Lycalopex sp.","Daphoenus vetus","Barbourofelis loveorum","Archaeocyon leptodus","Canis armbrusteri","Canis dirus","Canis latrans orcutti","Canis lupus (Pleistocene)","Desmocyon thomsoni","Hesperocyon gregarius","Mesocyon coryphaeus","Paraenhydrocyon josephi","Phlaocyon leucosteus","Vulpes macrotis (Pleistocene)","Vulpes vulpes (Pleistocene)","Homotherium ischyrus","Homotherium serum","Leopardus wiedii (Pleistocene)","Lynx rufus (Pleistocene)","Miracinonyx inexpectatus","Panthera atrox","Puma concolor (Pleistocene)","Puma lacustris","Smilodon fatalis","Miacis gracilis","Mephitis mephitis (Pleistocene)","Spilogale gracilis (Pleistocene)","Spilogale putorius (Pleistocene)","Gulo gulo (Pleistocene)","Martes americana (Pleistocene)","Martes nobilis (Pleistocene)","Mustela nigripes (Pleistocene)","Neovison frenata (Pleistocene)","Neovison vison (Pleistocene)","Satherium piscinarium","Taxidea taxus (Pleistocene)","Dinictis felina","Dinictis major","Hoplophoneus primaevus","Nimravus brachyops","Agriotherium africanum","Arctodus simus","Ursus arctos (Pleistocene)")
v_indices_to_remove<-which(sapply(dfcarnivores_postcranial[,2],function(x,v_species_to_remove){is.element(x,v_species_to_remove)},v_species_to_remove=v_species_to_remove,simplify=TRUE))
dfcarnivores_postcranial<-dfcarnivores_postcranial[-v_indices_to_remove,]
m_names_change<-rbind(c("Alopex lagopus","Vulpes lagopus"),c("Lycalopex gymnocerus","Lycalopex gymnocercus"),c("Caracal serval","Leptailurus serval"),c("Felis silvestris libyca","Felis silvestris"),c("Atilax palundinosus","Atilax paludinosus"),c("Gallerella pulverulenta","Galerella pulverulenta"),c("Gallerella sanguinea","Galerella sanguinea"),c("Conepatus mesoleucus","Conepatus leuconotus"),c("Mephitis macroura vittata","Mephitis macroura"),c("Aonyx cinereus","Aonyx cinerea"),c("Paradoxurus hemaphroditus","Paradoxurus hermaphroditus"))
for (i in 1:nrow(m_names_change)){
    dfcarnivores_postcranial[which(dfcarnivores_postcranial[,2]==m_names_change[i,1]),2]<-m_names_change[i,2]
}
dfcarnivores_postcranial[,2]<-gsub(" ", "_", dfcarnivores_postcranial[,2])
row.names(dfcarnivores_postcranial)<-dfcarnivores_postcranial[,2]
dat<-dfcarnivores_postcranial[,c("Ecology","HuL","HuPCL","RaL")]

The subset consists of a categorization of locomotor habits (Ecology) and three variables associated with forelimb morphology, measured in mm (HuL = humerus length; HuPCL = deltopectoral crest length; RaL = radius length):

head(dat)
##                        Ecology    HuL HuPCL    RaL
## Ailurus_fulgens       arboreal 112.82 50.35  87.84
## Vulpes_lagopus      generalist 106.47 44.08 101.89
## Atelocynus_microtis generalist 116.01 52.55 107.78
## Canis_adustus        cursorial 127.84 45.03 135.86
## Canis_latrans        cursorial 160.06 67.03 168.32
## Canis_lupus          cursorial 212.12 92.11 210.94

The categorical variable includes six different locomotor preferences (generalist, scansorial, arboreal, semiaquatic, semifossorial, cursorial). The three morphological variables are functionally interesting because a larger dectopectoral crest relative to humerus length is indicative of a large shoulder moment, which facilitates mechanical advantage of the deltoid and pectoral muscles acting across the shoulder joint. A larger radius relative to humerus length is indicative of distal segments that are longer than proximal segments reflecting a large output lever (linked to high relative velocity transmissions) and overall limb elongation. The relative lengths of the radius and the dectopectoral crest are thus informative of speed and strength maximization where humerus length operates effectively as a scaling factor for both attributes. This setup can be explored in mvSLOUCH through multivariate models in which the lengths of the dectopectoral crest and the radius are specified as responses, and the length of the humerus is specified as explanatory variable in those models that allow continuous predictors (fundamentally OUBM models; see OUBM Section).

Samuels, Meachen, and Sakai (2013) included \(107\) extant carnivoran taxa in their dataset, but here we use a reduced subset of species after removing those which phenotypic data were not complete, defined at the species level, or available in the phylogenetic tree. The phylogeny is pruned from the mammalian calibrated tree presented by Hedges et al. (2015), which we download from the supplementary material made available at the Center for Biodiversity ( http://www.biodiversitycenter.org/ttol ). Then we rename species according to Wilson and Reeder (2005) and Johnson et al. (2006) for taxonomic matching, remove the tips not present in our phenotypic data, and finally arrange the rows in our data matrix according to the order of the tips in the phylogeny. In this case the last step is not actually needed but in general such a reordering is recommended for users to make sure the data correspond to the appropriate tips of the phylogeny.

b_correct_tree_download<-FALSE
tryCatch({
phyltree_mammals<-ape::read.tree("http://www.biodiversitycenter.org/files/ttol/8.TTOL_mammals_unsmoothed.nwk")
b_correct_tree_download<-TRUE
},error=function(e){cat("Could not download tree file! No analysis can be done! Vignette not built!")},
warning=function(w){cat("Problem with downloading tree file! No analysis can be done! Vignette not built!")}
)
phyltree_mammals$tip.label[which(phyltree_mammals$tip.label=="Uncia_uncia")]<-"Panthera_uncia"
phyltree_mammals$tip.label[which(phyltree_mammals$tip.label=="Parahyaena_brunnea")]<-"Hyaena_brunnea"
phyltree_mammals$tip.label[which(phyltree_mammals$tip.label=="Bdeogale_crassicauda")]<-"Bdeogale_jacksoni"
tips_todrop<-setdiff(phyltree_mammals$tip.label,rownames(dat))
PrunedTree<-ape::drop.tip(phyltree_mammals,tips_todrop)
Tree<-ape::di2multi(PrunedTree)
dat<-dat[Tree$tip.label,]

You will notice that the number of internal nodes (90) is unexpectedly low given the number of tips (105). This is because the tree includes polytomies. Polytomies are not a problem for mvSLOUCH as its design does not depend on the number of descendants. Including polytomies will not affect the values of the likelihood during optimization, and in fact can result in more stable estimations than when using phylogenies with very short branches (close to zero). Another practice that can lead to more stable estimations is using trees scaled to maximum tree height, as the smaller numerical branch length values allow the estimator to find the maximum likelihood peak more easily. As currently loaded, this tree is not scaled:

mvSLOUCH::phyltree_paths(Tree)$tree_height
## [1] 54.96021

This is the total tree height, i.e. the amount of time (in millions of years) from the root to the tips (the function mvSLOUCH::phyltree_paths() will be addressed in more detail in Models section). We can scale branch lengths by total tree height to make them a proportion of one:

tree_height<-mvSLOUCH::phyltree_paths(Tree)$tree_height
ScaledTree<-Tree
ScaledTree$edge.length<-ScaledTree$edge.length/tree_height
mvSLOUCH::phyltree_paths(ScaledTree)$tree_height
## [1] 1

Now the branches are proportionately scaled to the tree height (rather than absolute time), easing the estimation procedure. It is important to ensure that the taxa names in the data and in the tree match, and that the order of names correspond to each other:

isTRUE(all.equal(ScaledTree$tip.label,rownames(dat)))
## [1] TRUE

When this is not the case (i.e. when there is a list of mismatches instead of the “OK” sign), the data and/or the tree have to be pruned (see the ape package ape::drop.tip(), ape::keep.tip(), for example) or renamed accordingly.

Regime specification

We will use the ecological habitat preferences as hypothesized selective regime for the limb traits. mvSLOUCH implements an unordered parsimony algorithm for this purpose of reconstructing the regime states on the phylogeny. First, we need to ensure that the locomotor categories match the correct species names in the tree. This is typically not the case as the data frame and the phylogenetic tree are usually from independent sources and have species names listed in different orders:

row.names(dat)==ScaledTree$tip.label
##   [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [16] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [31] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [46] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [61] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [76] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [91] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE

So, we will store the locomotor categories in a new object where the order is specified according to tree tip names:

regimes<-dat$Ecology[order(match(row.names(dat), ScaledTree$tip.label))]

With this new object we can run the parsimony reconstruction:

regimesFitch<-mvSLOUCH::fitch.mvsl(ScaledTree,regimes)
## The number of ambiguous nodes are:
## [1] 8
##  and they are at nodes:
## [1]  4  5  6 24 25 26 87 88
## Try deltran OR acctran = TRUE in the function call in order to implement a delayed or accelerated character transformation

When the reconstruction involves ambiguous nodes (as in this case), mvSLOUCH offers two options for resolving the character optimization: ACCTRAN (accelerated transformations) and DELTRAN (delayed transformations). The former assigns changes as close to the root of the phylogenetic tree as possible (thus favoring reversals), and the latter as close to the tips as possible (thus favoring convergences). Here we will use DELTRAN for the purposes of illustration:

regimesFitch<-mvSLOUCH::fitch.mvsl(ScaledTree,regimes,deltran=TRUE)

mvSLOUCH also offers the possibility of fixing the root to a given character state, which is particularly useful if the root node is ambiguous. We can visualize the reconstruction by painting the branches of the phylogenetic tree with different colors:

  1. generalist as purple,
  2. cursorial as red,
  3. arboreal as green,
  4. scansorial as orange,
  5. semiaquatic as blue and
  6. semifossorial as brown

according to the reconstruction:

reg.col<-regimesFitch$branch_regimes
reg.col[reg.col=="generalist"]<-"purple"
reg.col[reg.col=="arboreal"]<-"green"
reg.col[reg.col=="cursorial"]<-"red"
reg.col[reg.col=="scansorial"]<-"orange"
reg.col[reg.col=="semiaquatic"]<-"blue"
reg.col[reg.col=="semifossorial"]<-"brown"
plot(ScaledTree, cex = 1,  edge.color = reg.col, edge.width=3.5, type="fan", font=4)

Main models

Before running analyses, we will log transform the morphological variables so that they are less susceptible to scaling effects:

mvData<-data.matrix(dat[,c("HuPCL","RaL","HuL")])
mvData<-log(mvData)

mvSLOUCH works faster when the phylogeny object includes information on the paths and distances for each node. This information can be obtained with the function mvSLOUCH::phyltree_paths():

mvStree<-mvSLOUCH::phyltree_paths(ScaledTree)

Brownian motion (BM)

Now we are ready to explore some models. Let us start with a multivariate Brownian motion model, under which the morphological variables accumulate variation over time in the absence of systematic selection towards deterministic optima. The main inputs are the modified data and tree that we just created:

BMestim<-mvSLOUCH::BrownianMotionModel(mvStree,mvData)

Key parameter estimates for this model are the diffusion component of the stochastic differential equation (important for obtaining the infinitesimal covariance matrix) and the ancestral trait values:

BMestim$ParamsInModel$Sxx
##           HuPCL       RaL        HuL
## HuPCL 0.7960162 0.0000000 0.00000000
## RaL   0.6373211 0.3075142 0.00000000
## HuL   0.5995487 0.2079952 0.09433783
BMestim$ParamsInModel$vX0
##           [,1]
## HuPCL 3.900998
## RaL   4.455671
## HuL   4.616998

Ornstein-Uhlenbeck Brownian-motion (OUBM)

We can also specify a model with a multivariate regression setup in which humerus length is used as a continuous random explanatory variable using the argument predictors in the mvSLOUCH::mvslouchModel() function. Here the predictor variable is modeled as a Brownian motion process on the phylogeny and the response variables are modeled as a multivariate Ornstein-Uhlenbeck process. A given response trait’s optimum is affected both by the other response variable trait’s optimum as well as the state of the predictor variable. Humerus length is the \(3^{\mathrm{rd}}\) variable in the data matrix, where the first two are the response variables (indicated with the argument kY), and the model is set up as follows:

OU1BM<-mvSLOUCH::mvslouchModel(mvStree, mvData, kY = 2, predictors = c(3))

The output for this model has more components than the simpler BM only model as there are many more parameters estimated. We will cover several of them later on (Numerical optimization section), but describe some key ones here starting with the rate of adaptation matrix \((\mathbf{A})\):

OU1BM$FinalFound$ParamsInModel$A 
##           HuPCL         RaL
## HuPCL 0.1859072 -1.68908585
## RaL   0.2906626  0.09983092

This matrix contains information about phylogenetic inertia, which is easier to interpret with the half-lives:

OU1BM$FinalFound$ParamSummary$phyl.halflife$halflives
##                               [,1]                   [,2]
## eigenvalues   0.1428691+0.6993581i   0.1428691-0.6993581i
## halflife      4.8516250+0.0000000i   4.8516250+0.0000000i
## %treeheight 485.1625000+0.0000000i 485.1625000+0.0000000i

#Under this model, it takes close to five times the tree height for the response variables to lose half of their ancestral influence. #So, if we assume that humerus length evolves randomly, it takes a very long time for the response variables to track it. #This can be confirmed by another set of key parameters in this model, the optimal and evolutionary regressions: Under this model, it takes close to five times the tree height for the response variables to lose half of their ancestral influence. Note that, unlike the \(\mathbf{A}\) matrix, the half-life entries are numbered ([,1] and [,2]) rather than tied to specific variables (HuPCL and RaL). This is because half-lives are reported in the eigenvector directions rather than in trait space. So, if we assume that humerus length evolves randomly, it takes a very long time for the response variables to track it. This can be confirmed by another set of key parameters in this model, the optimal and evolutionary regressions:

OU1BM$FinalFound$ParamSummary$optimal.regression
##             HuL
## HuPCL 2.6686996
## RaL   0.5232399
OU1BM$FinalFound$ParamSummary$evolutionary.regression
##              HuL
## HuPCL 0.03673609
## RaL   0.40216998

The optimal regression describes the predicted association if the responses (HuPCL and RaL) could adapt instantaneously to changes in the explanatory variable free of ancestral trait influences (HuL). The evolutionary regression shows the observed relationship, after accounting for general phylogenetic effects. The two sets of coefficients are very different, indicating that the observed association is far shallower than the theoretical expectation of instantaneous adaptation. Consistent with the half-lives, this suggests that changes in the response variables towards a randomly evolving humerus length are very slow. However, note that this model also assumes that the carnivorans under consideration evolve in the same type of environment. This can be easily recognized by checking the deterministic part of the primary optimum of the response variables:

OU1BM$FinalFound$ParamsInModel$mPsi
##           reg.1
## HuPCL -8.461001
## RaL    2.021930

There is a single regime for the primary optimum, but if the habitats these carnivorans occupy have had any impact on the evolution of their limbs, several niches should be considered (one for each habitat type, analogous to MANCOVA). We can account for this habitat contribution by using the locomotor preferences as a selective regime (with the argument regimes below; note that we are calling the regimesFitch object created earlier in the Regime specification section). Note also that we are specifying the niche at the root of the tree as “generalist” (with the argument root.regime below) as indicated by the parsimony reconstruction (the base of the tree is purple in the figure in the Regime specification section, corresponding to the generalist niche):

OUBMestim <- mvSLOUCH::mvslouchModel(mvStree, mvData, kY = 2, predictors = c(3), regimes = regimesFitch$branch_regimes, root.regime = "generalist")

Before going further, note there is a big difference in the output of this model compared with the Brownian motion one concerning the likelihood calculations. Under the OUBM model, two sets of outputs are reported:

OUBMestim$FinalFound
OUBMestim$MaxLikFound

The former (OUBMestim$FinalFound) stores the estimates corresponding the point where the likelihood optimization stopped. The latter (OUBMestim$MaxLikFound) stores the estimates under the maximum likelihood point found during the optimization. When these points are the same (as in OU1BM), mvSLOUCH will report it explicitly ("Same as final found"). When they are not, the outputs for each point will be stored separately. The discrepancy between the two points is indicative of likelihood convergence issues. We will discuss the issue of convergence in the Model comparison section. However, this does not seem to be the case for the current model (OUBMestim) and we may start comparing this model (OUBMestim) with the single regime specification (OU1BM):

OUBMestim$FinalFound$ParamsInModel$mPsi
##        arboreal cursorial generalist scansorial semiaquatic semifossorial
## HuPCL -7.368512 52.353475   5.391112  -2.850995   52.129339     -7.326189
## RaL    4.353962  4.984678   4.354435   4.663430    4.200482      3.932279

The model estimated a very different deterministic part of the primary optimum for each variable under particular locomotor types. This habitat contribution will also affect all other parameter estimates. The half-lives are now:

OUBMestim$FinalFound$ParamSummary$phyl.halflife$halflives
##                     [,1]         [,2]
## eigenvalues 8.224048e+03 2.587842e-02
## halflife    8.428297e-05 2.678476e+01
## %treeheight 8.428297e-03 2.678476e+03

When accounting for the locomotor types, the response variables lose their ancestral effects faster. The lower phylogenetic inertia is also observed in terms of the optimal regression relationship with the explanatory variable:

OUBMestim$FinalFound$ParamSummary$optimal.regression
##                HuL
## HuPCL -0.308507324
## RaL    0.001366637
OUBMestim$FinalFound$ParamSummary$evolutionary.regression
##                HuL
## HuPCL -0.003959721
## RaL    0.001272818

The difference between the two regressions is less extreme here than in the single regime specification (OU1BM), in particular regarding radius length. So, the response variables (especially radius length) are less influenced (having smaller, in magnitude, regression coefficients) by evolving humerus length when locomotor types are accounted for. But what if humerus length does not evolve independently of the other two variables? Then using a BM process to model the evolution of humerus length would not be appropriate (as in OU1BM and OUBMestim), requiring a different model that we describe next.

Ornstein-Uhlenbeck Ornstein-Uhlenbeck (OUOU)

Instead of assuming that humerus length evolves as BM, we can model all variables as an OU process:

OU1OU <- mvSLOUCH::ouchModel(mvStree, mvData, predictors = c(3))

We can easily verify that humerus length is now following an OU process by checking the rate of adaptation matrix (note that, unlike the OUBM models, HuL is now part of the matrix):

OU1OU$FinalFound$ParamsInModel$A
##           HuPCL        RaL        HuL
## HuPCL 0.4036434  0.8025262 -1.8698215
## RaL   0.6923298  0.1095877 -0.7467194
## HuL   0.3180505 -0.3959240  0.1555626

We can now look at the half-life estimates for the vector of traits:

OU1OU$FinalFound$ParamSummary$phyl.halflife$halflives
##                      [,1]                    [,2]                    [,3]
## eigenvalues  0.9930676+0i   -0.1621369+0.4040448i   -0.1621369-0.4040448i
## halflife     0.6979859+0i   -4.2750725+0.0000000i   -4.2750725+0.0000000i
## %treeheight 69.7985910+0i -427.5072490+0.0000000i -427.5072490+0.0000000i

Note that two half-lives (second and third) are negative (as opposed to the OUBM model explored previously), associated with the negative eigenvalues from the \(\mathbf{A}\) matrix. The interpretation of these negative eigenvalues is tricky, although one way of looking at them is in terms of character displacement (Bartoszek et al. 2012). Before trying any sort of interpretation, however, let us take a look at the conditional on predictors non-phylogenetic \(\mathrm{R}^{2}\) that is computed by mvSLOUCH:

OU1OU$FinalFound$ParamSummary$RSS$R2_non_phylogenetic_conditional_on_predictors
## [1] 0.9799846

It seems much better, compared to the OUBM models explored above:

OU1BM$FinalFound$ParamSummary$RSS$R2_non_phylogenetic_conditional_on_predictors
## [1] 0.5961821
OUBMestim$FinalFound$ParamSummary$RSS$R2_non_phylogenetic_conditional_on_predictors
## [1] NaN

In the OUBMestin object we have a NaN value as there are huge numerical issues in the calculation. This is due to the optimizer ending in a local maximum with the \(\mathbf{A}\) matrix

OUBMestim$FinalFound$ParamsInModel$A
##            HuPCL        RaL
## HuPCL 0.02166171  -13.71051
## RaL   2.52932467 8224.05186

having entry [2,2] orders of magnitude larger than the other entries and resulting in one extremely short half-life (essentially instantaneous adaptation) and one extremely long half-life

OUBMestim$FinalFound$ParamSummary$phyl.halflife$halflives
##                     [,1]         [,2]
## eigenvalues 8.224048e+03 2.587842e-02
## halflife    8.428297e-05 2.678476e+01
## %treeheight 8.428297e-03 2.678476e+03

Such situations cause numerical issues when summarizing parameters, calculating summary statistics or conditional distributions. In fact, the value of the R2_non_phylogenetic_conditional_on_predictors could actually be negative. This could be due to the fact that the non-phylogenetic conditional \(\mathrm{R}^{2}\) is calculated under a completely different, non-nested, model (compared to the model in which the parameters are estimated). As the name of the field suggests, this statistic ignores the phylogenetic correlations. Hence, since the parameters were estimated with the phylogenetic correlations, it could happen that a model with fewer degrees of freedom (only the grand mean) had a lower residual sum of squares (but it is non-nested as it ignores the phylogenetic correlations). The reason why we calculate the conditional \(\mathrm{R}^{2}\) without the phylogeny is that the linear likelihood evaluation algorithm cannot be carried over to conditional distributions. However, the hope is that if the model fit is good, then the non-phylogenetic conditional \(R^{2}\) will be high, and tell us something about the explained variance. We will not interpret this model right now, as other models could explain the data better (Model comparison section). So, let us rather focus on comparing some of its outputs with the OUBM counterpart instead. Recall that under the OUBM model, the association of the responses with the explanatory variable was established by contrasting two regressions (i.e. the optimal and evolutionary regressions). Recall also that the optimal regression described a theoretical association in which the responses adapted instantaneously to a randomly evolving humerus length. But humerus length does not evolve randomly under the OUOU model, so this theoretical contrast is no longer in place and mvSLOUCH reports only the observed relationship:

OU1OU$FinalFound$ParamSummary$evolutionary.regression
##            HuL
## HuPCL 1.358657
## RaL   1.071534

In fact, the primary optimum value can now be estimated for humerus length, given that it does not evolve under BM under this model:

OU1OU$FinalFound$ParamsInModel$mPsi
##          reg.1
## HuPCL 3.891970
## RaL   4.444546
## HuL   4.604391

There is a single primary optimum value for each morphological variable, reflecting the constant regime we specified in the OUOU model. Let us now a fit an OUOU model accounting for the locomotor preferences as selective regime:

OUOUestim <- mvSLOUCH::ouchModel(mvStree, mvData, regimes = regimesFitch$branch_regimes, root.regime = "generalist", predictors = c(3))

Now we have several regimes, one for each locomotor category:

OUOUestim$FinalFound$ParamsInModel$mPsi
##         arboreal cursorial generalist scansorial semiaquatic semifossorial
## HuPCL -89.124658 482.07835   3.957083 -68.227523   227.91093    -185.75094
## RaL   -35.894705 212.48162   4.542468 -26.790948   100.73129     -78.00468
## HuL    -9.651532  78.68383   4.653273  -6.599909    39.40867     -24.95484

Before trying to interpret the primary optima on different locomotor types, let us take a look at the eigenvalues from the \(\mathbf{A}\) matrix:

OUOUestim$FinalFound$ParamSummary$phyl.halflife$halflives
##                               [,1]                   [,2]            [,3]
## eigenvalues   0.3141781+0.3541126i   0.3141781-0.3541126i 9.373519e-03+0i
## halflife      2.2062240+0.0000000i   2.2062240+0.0000000i 7.394738e+01+0i
## %treeheight 220.6223997+0.0000000i 220.6223997+0.0000000i 7.394738e+03+0i

and at the the output of the non-phylogenetic \(\mathrm{R}^{2}\):

OUOUestim$FinalFound$ParamSummary$RSS$R2_non_phylogenetic_conditional_on_predictors
## [1] 0.5983495

The evolutionary regression has positive coefficients:

OUOUestim$FinalFound$ParamSummary$evolutionary.regression
##            HuL
## HuPCL 1.626443
## RaL   1.557123

We have up to now merely scratched the surface of the main models implemented by mvSLOUCH. Before saying anything about their performance, we need to conduct a more thorough model comparison, which is the topic of the next section.

Model comparison

In the previous section we explored a basic set of models applied on the carnivoran dataset. These models differ not only in their assumptions but also in their level of complexity, as revealed by their degrees of freedom:

cbind(BMestim$ParamSummary$dof,
      OU1BM$FinalFound$ParamSummary$dof,
      OU1OU$FinalFound$ParamSummary$dof,
      OUBMestim$FinalFound$ParamSummary$dof,
      OUOUestim$FinalFound$ParamSummary$dof)
##      [,1] [,2] [,3] [,4] [,5]
## [1,]    9   13   18   23   33

BM (1) is the simplest candidate, and the OU models accounting for the selective regime (4, 5) are the most complex. However, the OU models themselves can vary considerably in complexity depending on their parameter specifications. mvSLOUCH offers a variety of parameter specifications that allows adjusting the level of complexity of the models as well as contrasting different evolutionary scenarios on the data. The main two arguments involved in this parameter modification are Atype and Syytype (see below for an example), which allow setting the type of \(\mathbf{A}\) and \(\mathbf{\Sigma}_{yy}\) matrices, respectively. mvSLOUCH sets up \(\mathbf{\Sigma}_{yy}\) as "UpperTri" by default (i.e. upper triangular), which can be easily verified by checking this matrix in any of the OUOU models fitted above (this one in particular corresponds to the OUOU model with the locomotor regime):

OUOUestim$FinalFound$ParamsInModel$Syy
##           HuPCL       RaL      HuL
## HuPCL 0.4298954 0.1447712 1.247322
## RaL   0.0000000 0.2808094 1.194409
## HuL   0.0000000 0.0000000 1.134235

Other options are: "SingleValueDiagonal", "Diagonal", "LowerTri", "Symmetric", and "Any". In the case of the \(\mathbf{A}\) matrix, mvSLOUCH sets it up as "Invertible" by default, with other options being: "SingleValueDiagonal", "Diagonal", "UpperTri", "LowerTri", "Symmetric", "SymmetricPositiveDefinite", "DecomposablePositive", "DecomposableNegative", "DecomposableReal", "TwoByTwo", and "Any". The default setting (i.e. “Invertible”) is very general and makes the fewest biological assumptions on the traits, but can often lead the estimation procedure to getting stuck at a local likelihood peak. Let us try a model with the matrices \(\mathbf{\Sigma}_{yy}\) as lower triangular, and \(\mathbf{A}\) as diagonal:

OUBMestim.mod <- mvSLOUCH::mvslouchModel(mvStree, mvData, kY = 2, predictors = c(3), regimes = regimesFitch$branch_regimes, root.regime = "generalist", Syytype = "LowerTri", Atype = "Diagonal")

The maximum likelihood and final points are the same:

OUBMestim.mod$MaxLikFound
## [1] "Same as final found"

It is important to keep in mind, however, that no parameter specification guarantees attaining the maximum likelihood peak. Additional measures should be taken towards this end, such as increasing the number of iterations (which can be specified with the argument maxiter, see the manual for details) or conducting the search from different starting points (e.g. running the analysis several times, as a new starting point will be used for each run). It starts to become obvious that a thorough model comparison is challenging, not only because several parameter combinations are possible, but also because each of these combinations should be ran from different starting points. mvSLOUCH facilitates this process by offering a wrapper function that runs different types of models on the data from different starting points, which we will explore next.

Global regime

Let us apply the wrapper function to the constant regime for the whole tree:

OU1 <- mvSLOUCH::estimate.evolutionary.model(mvStree, mvData, repeats = 5, model.setups = "basic", predictors = c(3), kY = 2, doPrint = TRUE)

This command takes some time to run, as many models are fitted consequently. By setting the doPrint argument to TRUE, we can visualize what type of model is being fitted at each moment (this visualization can be omitted by leaving the default option: FALSE). The different OU model settings are specified through the argument model.setups. The "basic" option, despite being the simplest, is sufficient for most purposes and was selected for the current example. Other options increase progressively the model combinations to be tried in the following order: "fundamental", "extended", and finally "all" (the latter taking considerable time to run, as all possible model combinations are tried). The option "basic" consists of all possible combinations of "Diagonal" and "UpperTri" settings for Syytype, with "Diagonal", "UpperTri", "LowerTri", "DecomposablePositive", and "DecomposableReal" settings for Atype. Each of these settings (\(10\) combinations) is tried for OUBM and OUOU models (\(20\) combinations in total), and each combination is run from the number of starting points specified in repeats. Since we specified \(5\) starting points in this case, we will have \(100\) OU models (the \(20\) combinations described earlier running from \(5\) different starting points), plus BM (for a total of \(101\) models). Thus, the output is large and it might be a good idea to store it in a file:

capture.output(OU1,file = "OU1.txt")

The best candidate from the \(101\) possibilities is an OUOU model (ouch) with diagonal \(\mathbf{A}\) and upper triangular \(\mathbf{\Sigma}_{yy}\):

OU1$BestModel$model
## $evolmodel
## [1] "ouch"
## 
## $Atype
## [1] "Diagonal"
## 
## $Syytype
## [1] "UpperTri"
## 
## $diagA
## [1] "Positive"
## 
## $parameter_signs
## list()

Under the "basic" setting of model.setups, the diagonal elements of \(\mathbf{A}\) (diagA) are always positive ("Positive") and the signs of other parameters (parameter_signs) are not coerced in any particular way (and thus the list is empty). The argument model.setups also allows modifying these signs by providing customized lists (see manual for details), but this should be done with caution as some model specifications rely on particular sign settings. Therefore, the user should make sure that customized sign settings do not conflict other model specifications (e.g. Atype and Syytype settings). The model setting of this preferred candidate behaves better than the OUOU model under a global regime we fitted earlier (OU1OU; see section OUOU). The eigenvalues from the \(\mathbf{A}\) matrix are all positive for the current model:

OU1$BestModel$BestModel$ParamSummary$phyl.halflife$halflives
##                   [,1]       [,2]       [,3]
## eigenvalues  2.0552790  1.8363154  1.8175491
## halflife     0.3372521  0.3774663  0.3813637
## %treeheight 33.7252106 37.7466292 38.1363667

As well as the non-phylogenetic conditional on predictors \(\mathrm{R}^{2}\)

OU1$BestModel$key.properties$R2_non_phylogenetic_conditional_on_predictors
## [1] 0.9724327

The value of \(\mathrm{R}^{2}\) is similar to that of the earlier model (OU1OU) but the wrapper function has shown us that a different parameterization results in a simpler structure with a diagonal \(\mathbf{A}\)

OU1$BestModel$BestModel$ParamsInModel$A
##          HuPCL      RaL      HuL
## HuPCL 2.055279 0.000000 0.000000
## RaL   0.000000 1.836315 0.000000
## HuL   0.000000 0.000000 1.817549

unlike the OU1OU case. However, from our simulations studies there seems to be some bias towards models with diagonal \(\mathbf{A}\) so a careful consideration of competing models has to be done. A comprehensive output for this improved model can be found in the BestModel field of the output list:

OU1$BestModel

The list also includes the outputs of all the compared models. These models can be found in the testedModels element of the list. For example, the best candidate corresponds to model \(11\) of the compared models:

OU1$BestModel$i
## [1] 11

So, you can also visualize the outputs of the preferred candidate in the list of tested models:

OU1$testedModels[[11]]

BM corresponds to model 21:

OU1$testedModels[[21]]
## $result
## $result$ParamsInModel
## $result$ParamsInModel$Sxx
##           HuPCL       RaL        HuL
## HuPCL 0.7960162 0.0000000 0.00000000
## RaL   0.6373211 0.3075142 0.00000000
## HuL   0.5995487 0.2079952 0.09433783
## 
## $result$ParamsInModel$vX0
##           [,1]
## HuPCL 3.900998
## RaL   4.455671
## HuL   4.616998
## 
## 
## $result$ParamSummary
## $result$ParamSummary$StS
##           HuPCL       RaL       HuL
## HuPCL 0.6336418 0.5073179 0.4772505
## RaL   0.5073179 0.5007432 0.4460665
## HuL   0.4772505 0.4460665 0.4116203
## 
## $result$ParamSummary$LogLik
## [1] 145.4019
## 
## $result$ParamSummary$dof
## [1] 9
## 
## $result$ParamSummary$m2loglik
## [1] -290.8038
## 
## $result$ParamSummary$aic
## [1] -272.8038
## 
## $result$ParamSummary$aic.c
## [1] -272.2136
## 
## $result$ParamSummary$sic
## [1] -239.0306
## 
## $result$ParamSummary$bic
## [1] -239.0306
## 
## $result$ParamSummary$RSS
## $result$ParamSummary$RSS$RSS
## [1] 315
## 
## $result$ParamSummary$RSS$R2
## [1] 0.676023
## 
## 
## $result$ParamSummary$confidence.interval
## $result$ParamSummary$confidence.interval$regression.summary
## $result$ParamSummary$confidence.interval$regression.summary$X0.regression.confidence.interval
##      Lower.end Estimated.Point Upper.end
## [1,]  3.264607        3.900998  4.537390
## [2,]  3.889940        4.455671  5.021402
## [3,]  4.104076        4.616998  5.129919
## 
## $result$ParamSummary$confidence.interval$regression.summary$regression.covariance.matrix
##            X0_1       X0_2       X0_3
## X0_1 0.10542711 0.08440898 0.07940629
## X0_2 0.08440898 0.08331505 0.07421780
## X0_3 0.07940629 0.07421780 0.06848655
## 
## $result$ParamSummary$confidence.interval$regression.summary$regression.confidence.interval.comment
## [1] "These are confidence intervals for parameters estimated by a GLS conditional on the A and diffusion matrix parameters. In the full covariance matrix of the regression estimators the matrix parameters are entered column wise for the deterministic optimum and row wise for the B matrix. Be careful if some of the parameters were set at fixed values in the user's input, i.e. check if the variances correctly correspond to the presented 95% CIs. These are calculated  as estimate =/- (qnorm(0.975), i.e. 0.975 quantile of N(0,1))*(square root of appropriate diagonal element of the covariance matrix), i.e. standard deviation."
## 
## 
## 
## 
## 
## $aic.c
## [1] -272.2136
## 
## $bic
## [1] -239.0306
## 
## $model
## $model$evolmodel
## [1] "bm"

The settings of all the models tried can be found in the model.setups element of the output list (where: “bm” = BM; “mvslouch” = OUBM; “ouch” = OUOU):

OU1$model.setups

Now, how did mvSLOUCH select among all these \(101\) candidates? It used information criteria, in particular, the second-order bias correction of the Akaike information criterion (AICc):

OU1$BestModel$BestModel$ParamSummary$aic.c
## [1] -275.8235

With lower values representing better model fit. Although AICc (aic.c) is used for identifying the top candidate (OU1$BestModel), other criteria are reported as well in case the user is interested in conducting alternative comparisons (aic = Akaike information criterion; sic = Schwarz information criterion; bic = Bayesian information criterion):

OU1$BestModel$BestModel$ParamSummary$aic
## [1] -276.8566
OU1$BestModel$BestModel$ParamSummary$sic
## [1] -231.8257
OU1$BestModel$BestModel$ParamSummary$bic
## [1] -231.8257

Full regime

The previous comparison was conducted among a set of models that did not account for locomotor preferences. We can do so by running a new comparison, this time accounting for the selective regime (and then comparing the results with the above outputs under the global regime):

OUf <- mvSLOUCH::estimate.evolutionary.model(mvStree, mvData, regimes = regimesFitch$branch_regimes, root.regime = "generalist", repeats = 5, model.setups = "basic", predictors = c(3), kY = 2, doPrint = TRUE)

Note that, other than adding information on the selective regimes (through the arguments regimes and root.regime), the specifications are the same as for the previous comparison (OU1). Once again, storing the results in file might be a good way of readily accessing the outputs later on:

capture.output(OUf, file = "OUf.txt")

Note that the preferred candidate under the new comparison corresponds to the same model setting of the uniform regime case:

OUf$BestModel$model
## $evolmodel
## [1] "ouch"
## 
## $Atype
## [1] "Diagonal"
## 
## $Syytype
## [1] "UpperTri"
## 
## $diagA
## [1] "Positive"
## 
## $parameter_signs
## list()

We now look at the the eigenvalues from the \(\mathbf{A}\) matrix and the non-phylogenetic conditional on predictors \(\mathrm{R}^{2}\):

OUf$BestModel$BestModel$ParamSummary$phyl.halflife$halflives
##                  [,1]       [,2]       [,3]
## eigenvalues  1.579286  1.4193091  1.3738891
## halflife     0.438899  0.4883694  0.5045146
## %treeheight 43.889901 48.8369432 50.4514635
OUf$BestModel$key.properties$R2_non_phylogenetic_conditional_on_predictors
## [1] 0.9407261

The difference between the two models is clearly reinforced by AICc, with the new candidate showing a remarkably lower value than the old one:

OUOUestim$FinalFound$ParamSummary$aic.c
## [1] -41.918
OUf$BestModel$BestModel$ParamSummary$aic.c
## [1] -283.5499

This newer model does also better than the top candidate of the wrapper function under a global regime (AICc = -275.8234892), although the difference of values in these two comparisons is considerable:

OUOUestim$FinalFound$ParamSummary$aic.c-OUf$BestModel$BestModel$ParamSummary$aic.c
## [1] 241.6319
OU1$BestModel$BestModel$ParamSummary$aic.c-OUf$BestModel$BestModel$ParamSummary$aic.c
## [1] 7.726385

Burnham and Anderson (2002) provided some rules of thumb (but see also Burnham, Anderson, and Huyvaert (2011)) for ranking support for candidate models based on such differences (\(\Delta\)AICc):

  1. Support \(0 - 2\): Substantial similarity
  2. \(4 - 7\): Considerably less similarity
  3. \(>10\): Essentially no similarity

Based on these rules of thumb, there is little empirical support for the top model of the global regime comparison (\(\Delta\)AICc = 7.7263854), and essentially no support for the OUOU accounting for locomotor types under default settings (\(\Delta\)AICc = 241.6318762), when compared with the preferred candidate of the new comparison (OUf$BestModel). These results, besides highlighting the advantages of a thorough model comparison (facilitated by the wrapper function of mvSLOUCH), suggest that locomotor preferences work effectively as a selective regime for the limb attributes considered in this example.

Lumped regimes

The preferred model presented above uses the regime specification inspired by the categorization of Samuels, Meachen, and Sakai (2013), with six locomotor ecologies. As introduced earlier (in the Data section), however, some morphological attributes are advantageous for different locomotor types, so it is possible that sets of niches have experienced similar selective pressures. For example, we mentioned how both scansorial and arboreal forms climb, and how both swimmers and diggers benefit from high force outputs of the limbs. If these locomotor types have experienced similar selective pressures, the model specified above (with six different niches) is probably too complex. This is not a trivial issue, as these multivariate models are inherently complex and every effort to simplify them is worth a try (this is one of the advantages of conducting the model comparison under the wrapper function). This can be achieved by lumping some of the niches together, and comparing the results with the outputs of the full regime specification. Besides avoiding overparameterization, the simplified regimes can offer better insights on the adaptive significance of the traits. For example, if a model lumping the semiaquatic and semifossorial niches has better fit than the full regime specification, it would be indicative that the selective pressure is more associated with the type of motion (stroking) than the habitat per se (i.e. water or soil). For this demonstration, we will compare the full regime specification with three simpler alternatives:

  1. lumping the arboreal and scansorial niches (climbers),
  2. lumping the semiaquatic and semifossorial niches (strokers), and
  3. combining the two (climbers and strokers).

Let us start by lumping the arboreal and scansorial niches:

climb.reg <- regimesFitch$branch_regimes
climb.reg[climb.reg=="arboreal"] <- "climber"
climb.reg[climb.reg=="scansorial"] <- "climber"

Now we have five niches: generalist, cursorial, climber, semiaquatic, semifossorial.

climb.col <- climb.reg
climb.col[climb.col=="generalist"] <- "purple"
climb.col[climb.col=="climber"] <- "green"
climb.col[climb.col=="cursorial"] <- "red"
climb.col[climb.col=="semiaquatic"] <- "blue"
climb.col[climb.col=="semifossorial"] <- "brown"
plot(ScaledTree, cex = 1,  edge.color = climb.col, edge.width=3.5, type="fan", font=4)

Let us conduct a model comparison under this regime specification:

OUc <- mvSLOUCH::estimate.evolutionary.model(mvStree, mvData, regimes = climb.reg, root.regime = "generalist", repeats = 5, model.setups = "basic", predictors = c(3), kY = 2, doPrint = TRUE)

The top candidate has the same parameter structure as the preferred models described above (OU1 and OUf):

OUc$BestModel$model
## $evolmodel
## [1] "ouch"
## 
## $Atype
## [1] "Diagonal"
## 
## $Syytype
## [1] "UpperTri"
## 
## $diagA
## [1] "Positive"
## 
## $parameter_signs
## list()

This consolidates the properties of this parameterization for describing the data. The best model of this comparison does better than the global regime (AICc = -275.8234892), but not as well as the full regime specification (AICc = -283.5498746):

OUc$BestModel$BestModel$ParamSummary$aic.c
## [1] -278.0746

So, at least for now, the extra complexity of the full regime specification is justified. Let us see if the same applies when lumping the semiaquatic and semifossorial niches:

strok.reg<-regimesFitch$branch_regimes
strok.reg[strok.reg=="semiaquatic"]<-"stroker"
strok.reg[strok.reg=="semifossorial"]<-"stroker"

Once again, we have five niches: generalist, cursorial, arboreal, scansorial, stroker.

strok.col<-strok.reg
strok.col[strok.col=="generalist"]<-"purple"
strok.col[strok.col=="stroker"]<-"blue"
strok.col[strok.col=="cursorial"]<-"red"
strok.col[strok.col=="arboreal"]<-"green"
strok.col[strok.col=="scansorial"]<-"orange"
plot(ScaledTree, cex = 1,  edge.color = strok.col, edge.width=3.5, type="fan", font=4)