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)

The model comparison under this new lumped regime specification:

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

The preferred candidate of this comparison confirms that the model type found in previous analyses has good properties for describing the data:

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

This model does better than the alternative lumped strategy (AICc = -278.0746499), but only extremely marginally better than the full regime specification (AICc = -283.5498746):

OUs$BestModel$BestModel$ParamSummary$aic.c
## [1] -283.5876

These two regime specifications (OUf and OUs) are supported then, although the lumped strategy provides, besides a slightly better fit, a simpler model:

OUf$BestModel$BestModel$ParamSummary$dof
## [1] 27
OUs$BestModel$BestModel$ParamSummary$dof
## [1] 24

Before saying anything decisive on this issue, however, let us fit an even simpler model by using a reduced regime specification combining the two lumping strategies described above:

red.reg<-regimesFitch$branch_regimes
red.reg[red.reg=="arboreal"]<-"climber"
red.reg[red.reg=="scansorial"]<-"climber"
red.reg[red.reg=="semiaquatic"]<-"stroker"
red.reg[red.reg=="semifossorial"]<-"stroker"

In this reduced regime we have four niches: generalist, cursorial, climber, stroker.

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

Model comparison for this reduced regime specification:

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

The data under the reduced regime is also well explained by the model type selected under other specifications:

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

But it is not as well supported as the previous model (AICc = -278.0746499) or as the full regime specification (AICc = -283.5498746):

OUr$BestModel$BestModel$ParamSummary$aic.c
## [1] -273.0162

So, the specifications that are better supported as adaptive regimes are

  1. the one with strokers as a lumped niche (OUs), and
  2. the one without lumped niches (OUf).

Numerical optimization

Given that the more complex model (OUf) is not helping us to explain the data better than the simpler model (OUs), we could narrow our attention on the latter. But as mentioned earlier (Model comparison section), there is no guarantee that the maximum likelihood peak has been reached, not even after using the wrapper function. It is possible that one of the two models (or both) are stuck at a local likelihood peak. A way to confirm this is using the previous outputs to conduct a more focused search in which the starting points are based on the optimized estimates of progressively complex models. We can start this search from a simple model (without regimes) that can then be used as starting point for both regime specifications:

OUOUstart<-mvSLOUCH::ouchModel(mvStree, mvData, predictors = c(3), Atype = "Diagonal", diagA = NULL)

We are using the ouchModel function because both preferred regime specifications (OUf and OUs) are based on an OUOU model of evolution. We are using the defaults for Syytype and a simple specification for Atype, allowing the signs of the \(\mathbf{A}\) matrix to vary (diagA = NULL). We will use the resulting \(\mathbf{A}\) and \(\mathbf{\Sigma}_{yy}\) matrices from this model:

OUOUstart$FinalFound$ParamsInModel$A
##          HuPCL     RaL      HuL
## HuPCL 1.512175 0.00000 0.000000
## RaL   0.000000 1.33721 0.000000
## HuL   0.000000 0.00000 1.278613
OUOUstart$FinalFound$ParamsInModel$Syy
##           HuPCL         RaL       HuL
## HuPCL 0.3173801 -0.08977163 1.1058162
## RaL   0.0000000  0.15781145 0.9962170
## HuL   0.0000000  0.00000000 0.9158166

As starting points for models specifying the selective regimes under the specifications identified by the wrapper function:

OptOUs1<-mvSLOUCH::ouchModel(mvStree, mvData, regimes = strok.reg, root.regime = "generalist", predictors = c(3), Atype = OUs$BestModel$model$Atype, Syytype = OUs$BestModel$model$Syytype, diagA = OUs$BestModel$model$diagA, start_point_for_optim=list(A = OUOUstart$FinalFound$ParamsInModel$A, Syy = OUOUstart$FinalFound$ParamsInModel$Syy))
OptOUf1<-mvSLOUCH::ouchModel(mvStree, mvData, regimes = regimesFitch$branch_regimes, root.regime = "generalist", predictors = c(3), Atype = OUf$BestModel$model$Atype, Syytype = OUf$BestModel$model$Syytype, diagA = OUf$BestModel$model$diagA, start_point_for_optim=list(A = OUOUstart$FinalFound$ParamsInModel$A, Syy = OUOUstart$FinalFound$ParamsInModel$Syy))

We accomplish this by using the argument start_point_for_optim. Note that both the reduced (regimes = strok.reg) and full (regimes = regimesFitch$branch_regimes) regime specifications retrieve Atype, Syytype, and diagA objects from the outputs of the wrapper function (OUs and OUf, respectively). Let us explore the likelihood for both models:

OptOUs1$FinalFound$LogLik
## [1] 170.0573
OptOUf1$FinalFound$LogLik
## [1] 172.7866

These models give us new \(\mathbf{A}\) and \(\mathbf{\Sigma}_{yy}\) matrices that we will use in turn as starting points for new analyses under each regime:

OptOUs2<-mvSLOUCH::ouchModel(mvStree, mvData, regimes = strok.reg, root.regime = "generalist", predictors = c(3), Atype = OUs$BestModel$model$Atype, Syytype = OUs$BestModel$model$Syytype, diagA = OUs$BestModel$model$diagA, start_point_for_optim=list(A = OptOUs1$FinalFound$ParamsInModel$A, Syy = OptOUs1$FinalFound$ParamsInModel$Syy))
OptOUf2<-mvSLOUCH::ouchModel(mvStree, mvData, regimes = regimesFitch$branch_regimes, root.regime = "generalist", predictors = c(3), Atype = OUf$BestModel$model$Atype, Syytype = OUf$BestModel$model$Syytype, diagA = OUf$BestModel$model$diagA, start_point_for_optim=list(A = OptOUf1$FinalFound$ParamsInModel$A, Syy = OptOUf1$FinalFound$ParamsInModel$Syy))

Note that this time we did not use OUOUstart (the simple model) to retrieve the \(\mathbf{A}\) and \(\mathbf{\Sigma}_{yy}\) matrices for both regimes, but the specific outputs obtained in the previous analyses (OptOUs1 for the reduced regime, and OptOUf1 for the full regime). Let us compare the likelihood with the previous analyses:

OptOUs2$FinalFound$LogLik
## [1] 170.0573
OptOUf2$FinalFound$LogLik
## [1] 172.7866

We can see that the likelihood remained the same. We repeat this procedure a few times (i.e. by the fifth time you use OptOUs2 and OptOUf2 as starting points) and you will notice that the likelihood does not change, implying that we are at a maximum:

OptOUs2$FinalFound$LogLik
## [1] 170.0573
OptOUf2$FinalFound$LogLik
## [1] 172.7866

However, if we compare these values with the likelihoods of the top models of the wrapper function:

OUs$BestModel$BestModel$LogLik
## [1] 167.8627
OUf$BestModel$BestModel$LogLik
## [1] 171.4091

It is easy to see that the likelihood values of the former are higher. We experiment what will happen if we try from a new starting point:

OUOUreStart<-mvSLOUCH::ouchModel(mvStree, mvData, predictors = c(3), Atype = "Diagonal", diagA = NULL)

The \(\mathbf{A}\) and \(\mathbf{\Sigma}_{yy}\) matrices are not strikingly different from the first attempt:

OUOUreStart$FinalFound$ParamsInModel$A
##          HuPCL     RaL      HuL
## HuPCL 1.370563 0.00000 0.000000
## RaL   0.000000 1.14815 0.000000
## HuL   0.000000 0.00000 1.112349
OUOUreStart$FinalFound$ParamsInModel$Syy
##           HuPCL         RaL       HuL
## HuPCL 0.3140709 -0.09064659 0.9107321
## RaL   0.0000000  0.15317158 0.8229434
## HuL   0.0000000  0.00000000 0.7574673

But let us see what happens with the likelihoods when we use these matrices as starting points for the analyses using selective regimes:

FinalOUs1<-mvSLOUCH::ouchModel(mvStree, mvData, regimes = strok.reg, root.regime = "generalist", predictors = c(3), Atype = OUs$BestModel$model$Atype, Syytype = OUs$BestModel$model$Syytype, diagA = OUs$BestModel$model$diagA, start_point_for_optim=list(A = OUOUreStart$FinalFound$ParamsInModel$A, Syy = OUOUreStart$FinalFound$ParamsInModel$Syy))
FinalOUf1<-mvSLOUCH::ouchModel(mvStree, mvData, regimes = regimesFitch$branch_regimes, root.regime = "generalist", predictors = c(3), Atype = OUf$BestModel$model$Atype, Syytype = OUf$BestModel$model$Syytype, diagA = OUf$BestModel$model$diagA, start_point_for_optim=list(A = OUOUreStart$FinalFound$ParamsInModel$A, Syy = OUOUreStart$FinalFound$ParamsInModel$Syy))

The likelihoods are very similar, in the first case slightly higher in the second slightly lower:

FinalOUs1$FinalFound$LogLik
## [1] 170.3498
FinalOUf1$FinalFound$LogLik
## [1] 172.5287

And if we recursively use the same matrices as starting points for new analyses:

FinalOUs2<-mvSLOUCH::ouchModel(mvStree, mvData, regimes = strok.reg, root.regime = "generalist", predictors = c(3), Atype = OUs$BestModel$model$Atype, Syytype =OUs$BestModel$model$Syytype, diagA = OUs$BestModel$model$diagA, start_point_for_optim=list(A = FinalOUs1$FinalFound$ParamsInModel$A, Syy = FinalOUs1$FinalFound$ParamsInModel$Syy))
FinalOUf2<-mvSLOUCH::ouchModel(mvStree, mvData, regimes = regimesFitch$branch_regimes, root.regime = "generalist", predictors = c(3), Atype = OUf$BestModel$model$Atype, Syytype = OUf$BestModel$model$Syytype, diagA = OUf$BestModel$model$diagA, start_point_for_optim=list(A = FinalOUf1$FinalFound$ParamsInModel$A, Syy = FinalOUf1$FinalFound$ParamsInModel$Syy))

we can see that no further improvements are possible:

FinalOUs2$FinalFound$LogLik
## [1] 170.3498
FinalOUf2$FinalFound$LogLik
## [1] 172.5287

As before the customized procedure attained a higher likelihood than the wrapper function. By comparing the AICc values of these optimized models, you can notice that the reduced model is more decisively preferred this time:

FinalOUf2$FinalFound$ParamSummary$aic.c - FinalOUs2$FinalFound$ParamSummary$aic.c
## [1] 2.772658

Confirming that the reduced model is doing a better job than the full model, and that we can focus our attention on the former for interpretation. In regards to strokers, it seems like the type of motion works more effectively as selective regime than the particular medium being pulled (water or soil). \(\mathrm{R}^{2}\) is low, indicating that there is, as often in phylogenetic comparative studies, considerable variation in the data left to be still explained:

FinalOUs2$FinalFound$ParamSummary$RSS$R2
## [1] 0.1140503

So now we can go back and explore the phenomena described at the beginning using this model (Data section). Let us start with the adaptive significance of limb morphology in terms of locomotor types. This can be explored by comparing the estimated primary optima for the different niches. mvSLOUCH not only reports the estimates, but also \(95\%\) generalized least squares (GLS) confidence intervals conditional on \(\mathbf{A}\) and diffusion matrix parameters:

FinalOUs2$FinalFound$ParamSummary$confidence.interval$regression.summary$mPsi.regression.confidence.interval
## $Lower.end
##       arboreal cursorial generalist scansorial  stroker
## HuPCL 2.991739  3.602147   3.550271   2.958096 3.049015
## RaL   3.465628  4.930851   4.187941   3.552501 3.451865
## HuL   3.804471  4.731460   4.332871   3.638942 3.625662
## 
## $Estimated.Point
##       arboreal cursorial generalist scansorial  stroker
## HuPCL 3.647465  4.957143   3.895716   3.846989 3.830358
## RaL   4.066752  6.172462   4.497850   4.368488 4.171164
## HuL   4.400474  5.961315   4.631353   4.448900 4.341989
## 
## $Upper.end
##       arboreal cursorial generalist scansorial  stroker
## HuPCL 4.303191  6.312140   4.241160   4.735882 4.611702
## RaL   4.667876  7.414074   4.807759   5.184474 4.890464
## HuL   4.996477  7.191171   4.929834   5.258859 5.058317

This object lists the values of the primary optima (field Estimated.Point) as well as the lower (field Lower.end) and upper (field Upper.end) bounds of the \(95\%\) confidence interval. Nothing conclusive can be said about dectopectoral crest (HuPCL) and humerus lengths (HuL), as the confidence intervals for the primary optima of the different niches overlap to some extent. The same does not happen with radius length (RaL) though, where the confidence intervals reveal more distinguishable differences among niches. To see this more clearly, let’s plot these estimates. First we will use the above object (listing the primary optimum estimates with confidence intervals) to create a data frame (DFs) indicating the locomotor habits (Ecology) as a factor, and the radius length primary optimum (RaL) with the lower (lower) and upper (upper) bounds of the confidence intervals:

DFs<-data.frame(
 Ecology=factor(colnames(FinalOUs2$FinalFound$ParamSummary$confidence.interval$regression.summary$mPsi.regression.confidence.interval$Estimated.Point)),
  RaL=FinalOUs2$FinalFound$ParamSummary$confidence.interval$regression.summary$mPsi.regression.confidence.interval$Estimated.Point["RaL",],
  upper=FinalOUs2$FinalFound$ParamSummary$confidence.interval$regression.summary$mPsi.regression.confidence.interval$Upper.end["RaL",],
  lower=FinalOUs2$FinalFound$ParamSummary$confidence.interval$regression.summary$mPsi.regression.confidence.interval$Lower.end["RaL",]
)

We can now create the plot with ggplot2 using this dataframe and keeping the same color coding that was specified earlier to map the locomotor niches on the tree (Lumped regimes section):

ggplot2::ggplot(DFs, ggplot2::aes(Ecology, RaL))+
  ggplot2::geom_point(size=4, colour=c("green","red","purple","orange","blue")) +
    ggplot2::geom_errorbar(ggplot2::aes(ymin=lower,ymax=upper),width=0.1,lwd=1.5, colour=c("green","red","purple","orange","blue"))+
  ggplot2::xlab("Locomotor habits")+
  ggplot2::ylab("RaL(log)") + ggplot2::coord_flip()

The primary optimum for radius length is significantly higher for cursorials when compared to arboreals, generalists, and strokers. Although scansorials also exhibit a low optimum, its differentiation with the cursorial locomotion is less clear as their confidence intervals overlap (albeit not extensively). Recall that radius length is informative of output lever dynamics and overall limb elongation (Data section). The high value of the primary optimum for cursorial carnivorans reflects large output levers and limbs that tend to be long, increasing relative velocity transmissions and maximizing the distance covered by each stride (e.g. cheetahs, grey wolves, spotted hyenas). Arboreal, generalist, and stroker carnivorans have smaller output levers that result in higher mechanical advantage, indicative of higher force outputs when compared to cursorials. Not unexpectedly, arboreal and scansorial carnivorans group together in terms of the primary optima, as both locomotor types involve climbing abilities. But as indicated above, the scansorial locomotor type is the least differentiated with the cursorial one. This makes sense considering that some scansorial species rely on speed to some degree for hunting (e.g. pumas, snow leopards), similar to cursorial carnivorans.

It seems then that the trade-off between speed and strength is well reflected in the radius of carnivorans with different locomotor habits. But what happens with dectopectoral crest? Although this morphological feature is not distinguishing locomotor types as well as radius length, we would expect the two variables to be at odds with each other under the trade-off scenario, considering that relatively large values of the former favor high mechanical advantage, while relatively large values of the latter favor high velocity transmissions. We can explore this aspect of the trade-off by inspecting the association between the two variables in more detail. The evolutionary regression and overall correlations indicate that these variables are positively associated:

FinalOUs2$FinalFound$ParamSummary$evolutionary.regression
##            HuL
## HuPCL 1.140198
## RaL   1.047589
FinalOUs2$FinalFound$ParamSummary$corr.matrix
##           HuPCL       RaL       HuL
## HuPCL 1.0000000 0.9126224 0.9423476
## RaL   0.9126224 1.0000000 0.9851982
## HuL   0.9423476 0.9851982 1.0000000

Most likely though, this positive association reflects an absolute effect due to scaling. Basically, this strong positive association would indicate that as animals get larger, their limb measurements increase too. The story changes when we explore the conditional regression and correlation coefficients:

FinalOUs2$FinalFound$ParamSummary$trait.regression
## [[1]]
##              RaL      HuL
## HuPCL -0.6109474 1.780221
## 
## [[2]]
##          HuPCL      HuL
## RaL -0.1238151 1.188763
FinalOUs2$FinalFound$ParamSummary$conditional.corr.matrix
##            HuPCL        RaL
## HuPCL  1.0000000 -0.2750355
## RaL   -0.2750355  1.0000000

Note that the association with the scaling factor (HuL) remains positive, but now the association of the responses (HuPCL and RaL) is negative. This is, relative dectopectoral crest and humerus lengths are inversely associated, consistent with the trade-off hypothesis. We can assess the strength of this pattern with confidence intervals, but these are trickier to obtain for the above parameters than it was for the primary optimum. mvSLOUCH offers an alternative way of estimating confidence intervals for these and many other parameters, as presented in the next section.

Parametric bootstrap

The optimized speed of the new mvSLOUCH version allows obtaining confidence intervals for a variety of parameters under parametric bootstrapping. Here we will focus on the estimates linked to the trade-off pattern explored at the end of the previous section (evolutionary.regression, corr.matrix, trait.regression, conditional.corr.matrix):

BT<-mvSLOUCH::parametric.bootstrap(estimated.model = FinalOUs2, phyltree = mvStree, 
values.to.bootstrap = c("evolutionary.regression", "corr.matrix", "trait.regression", "conditional.corr.matrix"), 
regimes = strok.reg, root.regime = "generalist", predictors = c(3), numboot = 1000, 
Atype = OUs$BestModel$model$Atype, 
Syytype = OUs$BestModel$model$Syytype, 
diagA = OUs$BestModel$model$diagA)

The bootstrap function (parametric.bootstrap) uses the estimates of a fitted evolutionary model (estimated.model), which in our case results from the optimizing procedure of the previous section (FinalOUs2). The function also requires the specification of the phylogenetic tree (phyltree), the number of bootstrap samples (numboot), and the list of estimates for which we want obtain confidence intervals (values.to.bootstrap). Other arguments for this function have been introduced before and deal with the specification of the regime (regimes, root.regime), explanatory variable (predictors), and model setting (Atype, Syytype, diagA). The bootstrap procedure will take some time and mvSLOUCH will print the current iteration and the elapsed time so that you can monitor its progress. When it is completed, the resulting object (BT) will hold the simulated data and model outputs of each bootstrap replicate. Now we need to retrieve the confidence intervals from this object. We will do so by computing the \(2.5\%\) values on both tails to determine the \(95\%\) confidence region from the collection of estimates. Let us start with the evolutionary regression:

BTU.EvoReg <- FinalOUs2$FinalFound$ParamSummary$evolutionary.regression
BTU.EvoReg[] <- 0L
BTL.EvoReg <- BTU.EvoReg
for(i in 1:nrow(FinalOUs2$FinalFound$ParamSummary$evolutionary.regression)){
  BT.EvoReg<-quantile(sapply(BT$bootstrapped.parameters$evolutionary.regression,function(x) 
    x[i]), c(0.025, 0.975))
  BTL.EvoReg[i,] <- BT.EvoReg[1]
  BTU.EvoReg[i,] <- BT.EvoReg[2]
}

The lower (BTL.EvoReg) and upper (BTU.EvoReg) bounds of the \(95\%\) confidence interval for the evolutionary regression:

BTL.EvoReg
##                HuL
## HuPCL 5.570356e-05
## RaL   1.641088e-04
BTU.EvoReg
##            HuL
## HuPCL 1.317613
## RaL   1.316561

We can see that the confidence regions for both coefficients concentrate on positive values. The bootstrap procedure tends to produce wide confidence intervals, and thus is conservative. Overall, the positive general association between the responses and the explanatory variable is supported by the confidence intervals. Let’s see if the strength of this positive association is confirmed by the correlation matrix:

BTU.CorrMat <- rep(NA,length(as.vector(FinalOUs2$FinalFound$ParamSummary$corr.matrix)))
BTL.CorrMat<-BTU.CorrMat
for(i in 1:length(as.vector(FinalOUs2$FinalFound$ParamSummary$corr.matrix))){
  BT.CorrMat<-quantile(sapply(BT$bootstrapped.parameters$corr.matrix,function(x) x[i]),c(0.025,0.975))
  BTL.CorrMat[i] <- BT.CorrMat[1]
  BTU.CorrMat[i] <- BT.CorrMat[2]
}
BTL.CorrMat <- matrix(BTL.CorrMat, nrow =
  nrow(FinalOUs2$FinalFound$ParamSummary$corr.matrix))
BTU.CorrMat <- matrix(BTU.CorrMat, nrow =
  nrow(FinalOUs2$FinalFound$ParamSummary$corr.matrix))
dimnames(BTL.CorrMat) <- dimnames(BTU.CorrMat)<-
  list(row.names(FinalOUs2$FinalFound$ParamSummary$corr.matrix),
      colnames(FinalOUs2$FinalFound$ParamSummary$corr.matrix))

The lower (BTL.CorrMat) and upper (BTU.CorrMat) bounds of the \(95\%\) confidence interval for the correlation matrix:

BTL.CorrMat
##              HuPCL          RaL          HuL
## HuPCL 1.000000e+00 0.0006004333 6.343599e-05
## RaL   6.004333e-04 1.0000000000 2.802997e-04
## HuL   6.343599e-05 0.0002802997 1.000000e+00
BTU.CorrMat
##           HuPCL       RaL       HuL
## HuPCL 1.0000000 0.9553552 0.9693332
## RaL   0.9553552 1.0000000 0.9919731
## HuL   0.9693332 0.9919731 1.0000000

The estimates confirm the evolutionary regression results, with positive correlation coefficients. These results point towards a significant positive association among morphological variables, consistent with the absolute effects due to scaling discussed earlier. The trade-off can be addressed by studying the conditional regression coefficients:

NA.TrtReg<-lapply(1:length(BT$bootstrapped.parameters$trait.regression), function(x) 
  rep(NA,length(unlist(FinalOUs2$FinalFound$ParamSummary$trait.regression))))
BTU.TrtReg <- rep(NA, length(unlist(FinalOUs2$FinalFound$ParamSummary$trait.regression)))
BTL.TrtReg <- BTU.TrtReg
for(i in 1:length(unlist(FinalOUs2$FinalFound$ParamSummary$trait.regression))){
  BT.TrtReg<-quantile(sapply(relist(unlist(BT$bootstrapped.parameters$trait.regression),
                                     NA.TrtReg), function(x) x[i]), c(0.025, 0.975))
  BTL.TrtReg[i] <- BT.TrtReg[1]
  BTU.TrtReg[i] <- BT.TrtReg[2]
}
BTL.TrtReg <- relist(BTL.TrtReg, FinalOUs2$FinalFound$ParamSummary$trait.regression)
BTU.TrtReg <- relist(BTU.TrtReg, FinalOUs2$FinalFound$ParamSummary$trait.regression)

The lower (BTL.TrtReg) and upper (BTU.TrtReg) bounds of the \(95\%\) confidence interval for the conditional regression coefficients:

BTL.TrtReg
## [[1]]
##             RaL        HuL
## HuPCL -1.323227 -0.1773817
## 
## [[2]]
##          HuPCL        HuL
## RaL -0.4799471 -0.0142687
BTU.TrtReg
## [[1]]
##            RaL      HuL
## HuPCL 1.142619 2.521926
## 
## [[2]]
##         HuPCL     HuL
## RaL 0.8413354 1.67723

The negative association of the responses (HuPCL and RaL) cannot be considered significant as the confidence intervals include a large range of positive values. Let us see if the conditional correlation matrix confirms this result:

BTU.CondCorr <-
  rep(NA,length(as.vector(FinalOUs2$FinalFound$ParamSummary$conditional.corr.matrix)))
BTL.CondCorr<-BTU.CondCorr 
for(i in 1:length(as.vector(FinalOUs2$FinalFound$ParamSummary$conditional.corr.matrix))){
  BT.CondCorr<-quantile(sapply(BT$bootstrapped.parameters$conditional.corr.matrix,function(x) x[i]),c(0.025,0.975))
  BTL.CondCorr[i] <- BT.CondCorr[1]
  BTU.CondCorr[i] <- BT.CondCorr[2]
}
BTL.CondCorr <- matrix(BTL.CondCorr, nrow = 
  nrow(FinalOUs2$FinalFound$ParamSummary$conditional.corr.matrix))
BTU.CondCorr<-matrix(BTU.CondCorr,nrow =
  nrow(FinalOUs2$FinalFound$ParamSummary$conditional.corr.matrix))
dimnames(BTL.CondCorr)<-dimnames(BTU.CondCorr)<-list(row.names(FinalOUs2$FinalFound$ParamSummary$conditional.corr.matrix), 
  colnames(FinalOUs2$FinalFound$ParamSummary$conditional.corr.matrix))

The lower (BTL.CondCorr) and upper (BTU.CondCorr) bounds of the \(95\%\) confidence interval for the conditional correlation matrix:

BTL.CondCorr
##            HuPCL        RaL
## HuPCL  1.0000000 -0.6219967
## RaL   -0.6219967  1.0000000
BTU.CondCorr
##           HuPCL       RaL
## HuPCL 1.0000000 0.9205928
## RaL   0.9205928 1.0000000

Indeed, the conditional correlation between response variables includes a wide range of positive values in its confidence interval, and thus the negative association cannot be considered significant. Although this rejects the trade-off hypothesis between dectopectoral crest and radius lengths, it is important to keep in mind that these bootstrap confidence intervals are conservative. The trade-off might be present in the data, albeit weakly. The weakness of the trade-off could be also associated with a suboptimal selective regime specification. The locomotor categories assigned here are convenient because they belong to the same body of data as the morphological variables (Samuels, Meachen, and Sakai 2013). However, the assignation of Samuels, Meachen, and Sakai (2013), which is based on the proportion of time spent using different locomotor modes, generates some groupings that do not necessarily reflect similarity of motion in the forelimbs. For example, the polar bear (Ursus maritimus) is classified as semiaquatic, grouping this species with the otters, despite differences in limb usage by these taxa both on water and on land. Similarly, the generalist classification of the wolverine (Gulo gulo) might not reflect the digging capabilities of this species very well, particularly in the snow. Although the proportion of time spent in various locomotor modes is informative, a selective regime specification more tuned to the role of the forelimbs during motion might have the potential of revealing a clearer trade-off between relative dectopectoral crests and radius lengths. For now, though, all we can say is that the trade-off is better reflected by contrasting the radius of different locomotor types than by contrasting this structure with the dectopectoral crest.

References

Bartoszek, K., J. Pienaar, P. Mostad, S. Andersson, and T. F. Hansen. 2012. “A Phylogenetic Comparative Method for Studying Multivariate Adaptation.” J. Theor. Biol. 314: 204–15.

Burnham, K. P., and D. R. Anderson. 2002. Model Selection and Multimodel Inference: A Practical Information–Theoretic Approach. New York: Springer.

Burnham, K. P., D. R. Anderson, and K. P. Huyvaert. 2011. “AIC Model Selection and Multimodel Inference in Behavioral Ecology: Some Background, Observations, and Comparisons.” Behav. Ecol. Sociobiol. 65: 23–35.

Hedges, S. B., J. Marin, M. Suleski, M. Paymer, and S. Kumar. 2015. “Tree of Life Reveals Clock–Like Speciation and Diversification.” Mol. Biol. Evol. 32: 835–45.

Johnson, W. E., E. Eizirik, J. Pecon-Slattery, W. J. Murphy, A. Antunes, E. Teeling, and S. J. O’Brien. 2006. “The Late Miocene Radiation of Modern Felidae: A Genetic Assessment.” Science 311: 73–77.

Mitov, V., K. Bartoszek, G. Asimomitis, and T. Stadler. 2020. “Fast Likelihood Calculation for Multivariate Gaussian Phylogenetic Models with Shifts.” Theor. Pop. Biol. 131: 66–78.

R Core Team. 2020. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org.

Samuels, J. X., J. A. Meachen, and S. A. Sakai. 2013. “Postcranial Morphology and the Locomotor Habits of Living and Extinct Carnivorans.” J. Morphol. 274 (2): 121–46.

Wickham, H. 2016. ggplot2: Elegant Graphics for Data Analysis. New York: Springer.

Wilson, D. E., and D. M. Reeder. 2005. Mammal Species of the World: A Taxonomic and Geographic Reference. Baltimore, Maryland: Johns Hopkins University Press.