Back to blog

Orbit orientation isn't correlated with anything: a reanalysis of Heesy (2008)

In this post, I re-analyze data from a paper, “Ecomorphology of orbit orientation and the adaptive significance of binocular vision in primates and other mammals” by Dr. Christopher Heesy (2008). I demonstrate that correcting for the statistical non-independence of species data undermines the paper’s major conclusions, and I highlight directions for future work.

The ecomorphology of orbit orientation (Heesy 2008)

“Forward-facing orbits”, or the tendancy of the orbits of the skull to face forwards, have long been included in the suite of characteristics that differentiate the first primates from other early mammals. Other differentiating characteristics include hand and foot adaptations for grasping, hindlimb and ankle modifications for leaping, and dental features associated with eating more fruit than their predecessors.

Several hypotheses propose that orbit convergence was functionally important for early primates. The Grasp-Leaping Hypothesis posits that convergent orbits evolved to facilitate depth perception during rapid leaping in an arboreal environment (Szalay and Dagosto 1980). Several versions of the Angiosperm Diversification Hypothesis propose that convergent orbits facilitated visually directed feeding on fruits or insects on the terminal branches of angiosperms (e.g., Rasmussen 1990; Sussman 1991). Finally, the Nocturnal Visual Predation Hypothesis proposes that convergent orbits evolved to facilitate visually guided predation in a low-light environment (Cartmill 1992).

In a 2008 study, Chris Heesy analyzed a large dataset on orbit orientation and ecological variables in mammals. From his analysis, he concluded that “These data are entirely consistent with the nocturnal visual predation hypothesis of primate origins.” Below, I describe and reproduce Heesy’s (2008) study using data that was published along with the paper. I discuss Heesy’s results, then show how incorporating phylogeny into the statistical analysis undermines the paper’s major conclusions.

Heesy (2008) considered three ecological variables as potential predictors of orbit orientation across mammals: activity pattern (nocturnal = 1, other = 0), diet (insectivory = 1, other = 0), and substrate use (arboreal = 1, other = 0). He captured orbit orientation with three angles measured on the skull: convergence, frontation, and verticality. Important lines and planes used for computing these angles are diagrammed below (the Figure was copied from Heesy 2008):

Technical notes on quantifying orbit orientation: Orbit convergence is measured by the angle between the saggital and orbital planes, the latter being defined by the 3 points: OA- anterior orbit margin, OI- inferior orbit margin, and OS- superior orbit margin (this is shown in figure a above). Orbit frontation is measured by the angle between the line formed by the intersection of the saggital and orbital planes (shown in figure b above), and the inion-nasion line (not shown), which essentially runs along the midline from the back of the braincase through the forehead. Orbit verticality is the angle between the orbital plane and the palatal plane (formed by the prosthion and the midpoints of the M1 alveolar borders). All three of these measures capture the extent to which the orbits face ‘forwards’.

Heesy (2008) assessed the relationship between orbit orientation and ecological variables using MANOVA, with orbit convergence, frontation, and verticality as a multivariate response variable, and nocturnality, insectivory, and arboreality as binary predictor variables (plus interaction terms). This is the full model:

[convergence, frontation, verticality] ~ nocturnality * insectivory * arboreality

In addition to MANOVA, one-way ANOVAs were performed for pairs of independent and dependent variables.

To address the concern that certain taxonomic groups might be driving the results, Heesy repeated all analyses on two subsets of the data: 1) mammals excluding marsupials and anthropoid primates, and 2) mammals excluding marsupials and all primates. These particular groups were removed because it was argued that marsupials and primates (particularly anthropoids) have either special morphological constraints or highly derived orbit morphology (see the paper for details).

Now I will replicate Heesy’s analysis using the data from the paper, with the slight difference that I dropped 35 of Heesy’s original 331 taxa because they lacked phylogenetic data. I did this in order to facilitate a direct comparison between the results of phylogenetic and non-phylogenetic analyses. As the R code below demonstrates, my non-phylogenetic results are highly consistent with what Heesy (2008) reported in Table 1 of his paper, despite the removal of 35 taxa:

# import data
data <- read.csv("https://raw.githubusercontent.com/rgriff23/Heesy_2008_reanalysis/master/data/HeesyData.csv", header=TRUE, row.names=1)

# MANOVA/ANOVA (all taxa)
heesy.model.1 <- manova(cbind(Convergence, Frontition, Verticality) ~ Noccode*Fauncode*Arbcode, data=data)
summary.manova(heesy.model.1, test="Wilks")
>                            Df   Wilks approx F num Df den Df    Pr(>F)    
>  Noccode                    1 0.58495   65.989      3    279 < 2.2e-16 ***
>  Fauncode                   1 0.90222   10.079      3    279 2.516e-06 ***
>  Arbcode                    1 0.76720   28.220      3    279 5.740e-16 ***
>  Noccode:Fauncode           1 0.90691    9.546      3    279 5.072e-06 ***
>  Noccode:Arbcode            1 0.90855    9.361      3    279 6.467e-06 ***
>  Fauncode:Arbcode           1 0.96236    3.638      3    279   0.01331 *  
>  Noccode:Fauncode:Arbcode   1 0.97735    2.155      3    279   0.09357 .  
>  Residuals                281                                             
>  ---
>  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary.aov(heesy.model.1)
>   Response Convergence :
>                            Df Sum Sq Mean Sq F value    Pr(>F)    
>  Noccode                    1   2127  2127.1 11.0809 0.0009886 ***
>  Fauncode                   1   2329  2328.6 12.1307 0.0005750 ***
>  Arbcode                    1  14800 14800.3 77.1018 < 2.2e-16 ***
>  Noccode:Fauncode           1   2211  2211.0 11.5183 0.0007883 ***
>  Noccode:Arbcode            1   3528  3528.3 18.3807 2.489e-05 ***
>  Fauncode:Arbcode           1    467   467.0  2.4330 0.1199305    
>  Noccode:Fauncode:Arbcode   1   1154  1154.4  6.0138 0.0148029 *  
>  Residuals                281  53940   192.0                      
>  ---
>  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>  
>   Response Frontition :
>                            Df Sum Sq Mean Sq F value    Pr(>F)    
>  Noccode                    1   6481  6480.9 22.8581 2.817e-06 ***
>  Fauncode                   1   8092  8092.4 28.5416 1.897e-07 ***
>  Arbcode                    1   5290  5290.5 18.6592 2.170e-05 ***
>  Noccode:Fauncode           1   3245  3244.8 11.4443 0.0008191 ***
>  Noccode:Arbcode            1   1772  1772.1  6.2503 0.0129875 *  
>  Fauncode:Arbcode           1    330   329.6  1.1624 0.2818935    
>  Noccode:Fauncode:Arbcode   1    684   683.8  2.4118 0.1215460    
>  Residuals                281  79672   283.5                      
>  ---
>  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>  
>   Response Verticality :
>                            Df  Sum Sq Mean Sq  F value    Pr(>F)    
>  Noccode                    1 13352.4 13352.4 146.8381 < 2.2e-16 ***
>  Fauncode                   1  1095.5  1095.5  12.0479  0.000600 ***
>  Arbcode                    1   900.4   900.4   9.9013  0.001829 ** 
>  Noccode:Fauncode           1  1676.3  1676.3  18.4341 2.424e-05 ***
>  Noccode:Arbcode            1   831.8   831.8   9.1478  0.002720 ** 
>  Fauncode:Arbcode           1     8.6     8.6   0.0948  0.758401    
>  Noccode:Fauncode:Arbcode   1    72.8    72.8   0.8002  0.371798    
>  Residuals                281 25552.1    90.9                       
>  ---
>  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>  
>  7 observations deleted due to missingness
# MANOVA/ANOVA (drop marsupials and anthropoids)
data.eu1 <- data[!(data$ORDER%in%c("Didelphimorphia","Diprotodontia")),]
data.eu1 <- data.eu1[data$GROUP != "Anthropoidea",]
heesy.model.2 <- manova(cbind(Convergence, Frontition, Verticality) ~ Noccode*Fauncode*Arbcode, data=data.eu1)
summary.manova(heesy.model.2, test="Wilks")
>                            Df   Wilks approx F num Df den Df    Pr(>F)    
>  Noccode                    1 0.70945  27.7130      3    203 4.595e-15 ***
>  Fauncode                   1 0.86746  10.3386      3    203 2.313e-06 ***
>  Arbcode                    1 0.83768  13.1117      3    203 7.338e-08 ***
>  Noccode:Fauncode           1 0.94455   3.9727      3    203  0.008847 ** 
>  Noccode:Arbcode            1 0.96681   2.3230      3    203  0.076193 .  
>  Fauncode:Arbcode           1 0.95980   2.8343      3    203  0.039307 *  
>  Noccode:Fauncode:Arbcode   1 0.98915   0.7420      3    203  0.528149    
>  Residuals                205                                             
>  ---
>  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary.aov(heesy.model.2)
>   Response Convergence :
>                            Df Sum Sq Mean Sq F value    Pr(>F)    
>  Noccode                    1    374   374.3  1.7172   0.19152    
>  Fauncode                   1   1160  1159.9  5.3218   0.02206 *  
>  Arbcode                    1   7445  7445.0 34.1602 1.975e-08 ***
>  Noccode:Fauncode           1   1312  1311.9  6.0195   0.01498 *  
>  Noccode:Arbcode            1   1080  1080.0  4.9553   0.02710 *  
>  Fauncode:Arbcode           1    409   408.7  1.8751   0.17239    
>  Noccode:Fauncode:Arbcode   1    485   484.9  2.2250   0.13733    
>  Residuals                205  44679   217.9                      
>  ---
>  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>  
>   Response Frontition :
>                            Df Sum Sq Mean Sq F value    Pr(>F)    
>  Noccode                    1    842   841.6  2.9913 0.0852158 .  
>  Fauncode                   1   8689  8688.6 30.8838 8.458e-08 ***
>  Arbcode                    1   4343  4342.6 15.4359 0.0001166 ***
>  Noccode:Fauncode           1   1032  1032.2  3.6688 0.0568304 .  
>  Noccode:Arbcode            1     89    88.9  0.3159 0.5746996    
>  Fauncode:Arbcode           1    180   180.3  0.6407 0.4243747    
>  Noccode:Fauncode:Arbcode   1    228   228.2  0.8112 0.3688148    
>  Residuals                205  57673   281.3                      
>  ---
>  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>  
>   Response Verticality :
>                            Df  Sum Sq Mean Sq F value    Pr(>F)    
>  Noccode                    1  5021.6  5021.6 59.3361 5.581e-13 ***
>  Fauncode                   1  1116.0  1116.0 13.1873 0.0003561 ***
>  Arbcode                    1   806.1   806.1  9.5247 0.0023073 ** 
>  Noccode:Fauncode           1   563.8   563.8  6.6620 0.0105469 *  
>  Noccode:Arbcode            1    56.4    56.4  0.6663 0.4152830    
>  Fauncode:Arbcode           1    18.7    18.7  0.2214 0.6384377    
>  Noccode:Fauncode:Arbcode   1     8.4     8.4  0.0989 0.7534657    
>  Residuals                205 17349.1    84.6                      
>  ---
>  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>  
>  30 observations deleted due to missingness
# MANOVA/ANOVA (drop marsupials and all primates)
data.eu2 <- data[!(data$ORDER%in%c("Didelphimorphia","Diprotodontia","Primates")),]
heesy.model.3 <- manova(cbind(Convergence, Frontition, Verticality) ~ Noccode*Fauncode*Arbcode, data=data.eu2)
summary.manova(heesy.model.3, test="Wilks")
>                            Df   Wilks approx F num Df den Df    Pr(>F)    
>  Noccode                    1 0.74583  18.6298      3    164 1.896e-10 ***
>  Fauncode                   1 0.79615  13.9969      3    164 3.610e-08 ***
>  Arbcode                    1 0.91804   4.8804      3    164  0.002816 ** 
>  Noccode:Fauncode           1 0.98138   1.0374      3    164  0.377606    
>  Noccode:Arbcode            1 0.95960   2.3013      3    164  0.079143 .  
>  Fauncode:Arbcode           1 0.95074   2.8325      3    164  0.040024 *  
>  Noccode:Fauncode:Arbcode   1 0.99085   0.5049      3    164  0.679449    
>  Residuals                166                                             
>  ---
>  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary.aov(heesy.model.3)
>   Response Convergence :
>                            Df  Sum Sq Mean Sq F value    Pr(>F)    
>  Noccode                    1  4128.8  4128.8 40.6434 1.738e-09 ***
>  Fauncode                   1  1219.4  1219.4 12.0035 0.0006758 ***
>  Arbcode                    1     4.8     4.8  0.0471 0.8283762    
>  Noccode:Fauncode           1   134.8   134.8  1.3273 0.2509418    
>  Noccode:Arbcode            1   307.5   307.5  3.0270 0.0837449 .  
>  Fauncode:Arbcode           1    78.2    78.2  0.7702 0.3814214    
>  Noccode:Fauncode:Arbcode   1     3.7     3.7  0.0368 0.8480499    
>  Residuals                166 16863.2   101.6                      
>  ---
>  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>  
>   Response Frontition :
>                            Df Sum Sq Mean Sq F value   Pr(>F)   
>  Noccode                    1   1770 1769.56  6.7530 0.010201 * 
>  Fauncode                   1   2337 2336.91  8.9182 0.003251 **
>  Arbcode                    1    167  166.72  0.6363 0.426210   
>  Noccode:Fauncode           1    134  134.20  0.5121 0.475219   
>  Noccode:Arbcode            1    411  411.20  1.5692 0.212079   
>  Fauncode:Arbcode           1    733  733.22  2.7981 0.096259 . 
>  Noccode:Fauncode:Arbcode   1    167  167.28  0.6384 0.425434   
>  Residuals                166  43499  262.04                    
>  ---
>  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>  
>   Response Verticality :
>                            Df  Sum Sq Mean Sq F value    Pr(>F)    
>  Noccode                    1  1074.2 1074.19 16.6869 6.846e-05 ***
>  Fauncode                   1   301.2  301.16  4.6783   0.03198 *  
>  Arbcode                    1   105.8  105.76  1.6430   0.20171    
>  Noccode:Fauncode           1    48.5   48.50  0.7534   0.38667    
>  Noccode:Arbcode            1    37.6   37.61  0.5843   0.44571    
>  Fauncode:Arbcode           1    22.2   22.16  0.3442   0.55820    
>  Noccode:Fauncode:Arbcode   1    59.2   59.22  0.9199   0.33890    
>  Residuals                166 10685.9   64.37                      
>  ---
>  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>  
>  4 observations deleted due to missingness

Here is a run-down of the results (focusing on main effects):

This last finding- that arboreality is non-significant in the one-way ANOVAs when marsupials and all primates are dropped from the analysis- is the key result that Heesy (2008) interprets as support for the Nocturnal Visual Predation hypothesis and against hypotheses that focus on arboreality as a driving factor behind orbit evolution.

I have to say, I find it highly suspect to base this conclusion on a non-significant statistical result, particularly one which only appears after much of the data (and thus statistical power) has been removed. I also think it is noteworthy that the interaction between nocturnality and faunivory was not significant in the final model, since I would think that the Nocturnal Visual Predation Hypothesis should predict a significant interaction between these variables (not just independent effects). That said, I will now demonstrate that more appropriate phylogenetic models eliminate nearly all of the significant statistical results, rendering moot the interpretation of the particular arrangement of p-values in Heesy’s (2008) original results.

Problem of phylogenetic autocorrelation

Phylogenetic autocorrelation, also called phylogenetic pseudoreplication or phylogenetic non-independence, is a well-known problem in comparative evolutionary biology. The problem arises when researchers attempt to apply statistical methods that assume independent observations (e.g., linear regression) to evolutionary data that is not statistically independent. For example, a human, a chimpanzee, and a mouse cannot be considered 3 independent entities because a human and a chimpanzee are much more closely related to each other than they are to a mouse, thus they have had less time to evolve independently. Failure to account for phylogenetic autocorrelation in statistical models leads to elevated risk of Type I error, or false positives, because relationships that are driven by phylogenetic patterns may be incorrectly attributed to predictor variables in the model (Martins & Garland 1991).

Heesy was aware of this potential issue, but cited the absence of software for implementing phylogenetic MANOVA as the reason for using standard MANOVA. Lack of software might have been an issue when the paper was published in 2008, but fortunately it isn’t anymore. In the following section I use functions from the R packages geiger and geomorph to correct for phylogenetic autocorrelation in Heesy’s (2008) models, and discuss the impact it has on the results.

Reanalysis of Heesy (2008) with phylogeny

If you want to follow along with the rest of my analysis, install the following R packages (see info on installing geomorph):

# install packages
install.packages(ape)
install.packages(geiger)
install.packages(geomorph)

Next, import the phylogeny from GitHub. This is a trimmed down version of the mammal supertree published by Bininda-Emonds et al. (2007). I already did the work of cleaning of taxa labels and dropping taxa from Heesy’s data that were not present in the Bininda-Emonds phylogeny, so taxa labels will match up perfectly between the tree and data.

# load package
library(ape)

# import phylogenetic tree
tree <- read.nexus("https://raw.githubusercontent.com/rgriff23/Heesy_2008_reanalysis/master/data/HeesyTree.nexus")

# arbitrarily resolve polytomies
tree <- multi2di(tree)

We can perform a simple phylogenetic MANOVA using the aov.phylo function in the geiger package. This function has some limitations: it does not accept a data argument, and it only allows one predictor variable at a time. Still, we can use it to look at the MANOVA results for individual predictor variables:

# load package
library(geiger)

# define multidimensional response variable
angles <- data[,c("Convergence","Frontition","Verticality")]

# define univariate predictors (must be factors)
nocturnal <- as.factor(setNames(data$Noccode, row.names(data)))
insectivore <- as.factor(setNames(data$Fauncode, row.names(data)))
arboreal <- as.factor(setNames(data$Arbcode, row.names(data)))

# models
aov.phylo(angles ~ nocturnal, phy=tree)
>  Multivariate Analysis of Variance Table
>  
>  Response: dat
>             Df  Wilks approx-F num-Df den-Df     Pr(>F) Pr(>F) given phy
>  group       1 0.6222   58.697      3    290 1.1091e-29         0.000999
>  Residuals 292                                                          
>               
>  group     ***
>  Residuals    
>  ---
>  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>  Call:
>     manova(dat ~ group)
>  
>  Terms:
>                      group Residuals
>  resp 1            2431.92  83207.57
>  resp 2            6858.48 102631.62
>  resp 3           13954.52  30886.77
>  Deg. of Freedom         1       292
>  
>  Residual standard errors: 16.88068 18.74775 10.28478
>  Estimated effects may be unbalanced
>  2 observations deleted due to missingness
aov.phylo(angles ~ insectivore, phy=tree)
>  Multivariate Analysis of Variance Table
>  
>  Response: dat
>             Df   Wilks approx-F num-Df den-Df     Pr(>F) Pr(>F) given phy
>  group       1 0.90743    9.895      3    291 3.1196e-06           0.5574
>  Residuals 293
>  Call:
>     manova(dat ~ group)
>  
>  Terms:
>                      group Residuals
>  resp 1            3137.52  82574.85
>  resp 2            9815.18 100131.76
>  resp 3            1488.76  43697.67
>  Deg. of Freedom         1       293
>  
>  Residual standard errors: 16.78766 18.48639 12.21224
>  Estimated effects may be unbalanced
>  1 observation deleted due to missingness
aov.phylo(angles ~ arboreal, phy=tree)
>  Multivariate Analysis of Variance Table
>  
>  Response: dat
>             Df   Wilks approx-F num-Df den-Df     Pr(>F) Pr(>F) given phy
>  group       1 0.75702   30.599      3    286 3.4778e-17           0.1309
>  Residuals 288
>  Call:
>     manova(dat ~ group)
>  
>  Terms:
>                     group Residuals
>  resp 1          18162.23  62470.75
>  resp 2          11076.96  94938.79
>  resp 3           3086.58  40740.90
>  Deg. of Freedom        1       288
>  
>  Residual standard errors: 14.72794 18.15623 11.89376
>  Estimated effects may be unbalanced
>  6 observations deleted due to missingness

The output provides both standard (non-phylogenetic) and phylogenetic p-values in the last two columns of the MANOVA table. We can see that all of the standard p-values are highly significant, while the phylogenetic p-values are are much larger. Only nocturnality is significant here.

It would be better to fit models that estimate the effects of all three predictor variables simultaneously (as in Heesy’s MANOVAs). The procD.pgls function in the geomorph package allows us to perform phylogenetic regression with a multivariate response variable and multiple predictor variables. The function does not handle missing data, so we have to do some pre-processing before fitting the model. The model takes a few moments to run:

# load package
library("geomorph")

# subset data (no missing values allowed)
data.cc <- data[,c("Convergence","Frontition","Verticality","Noccode","Fauncode","Arbcode","GROUP","ORDER")]
data.cc <- data.cc[complete.cases(data.cc),]
angles2 <- data.cc[,1:3]
vars <- data.cc[,4:6]
droptax <- row.names(data)[!(row.names(data) %in% row.names(data.cc))]
tree2 <- multi2di(drop.tip(tree, droptax))

# create new 'geomorph data frame' with response, predictors, and tree
df <- geomorph.data.frame(angles2, vars, tree2)

# multivariate MANOVA with procD.pgls in geomorph
procD.pgls(angles2 ~ Noccode * Fauncode * Arbcode, phy=tree2, data=df, print.progress=F)
>  
>  Call:
>  procD.pgls(f1 = angles2 ~ Noccode * Fauncode * Arbcode, phy = tree2,  
>      data = df, print.progress = F) 
>  
>  
>  
>  Type I (Sequential) Sums of Squares and Cross-products
>  Randomized Residual Permutation Procedure Used
>  1000 Permutations
>  
>                            Df      SS      MS       Rsq      F       Z
>  Noccode                    1   18.82 18.8184 0.0084944 2.4288 0.54575
>  Fauncode                   1    0.76  0.7620 0.0003439 0.0983 0.02918
>  Arbcode                    1    9.81  9.8081 0.0044272 1.2659 0.95895
>  Noccode:Fauncode           1    0.81  0.8068 0.0003642 0.1041 0.06241
>  Noccode:Arbcode            1    4.12  4.1157 0.0018578 0.5312 0.40842
>  Fauncode:Arbcode           1    2.12  2.1159 0.0009551 0.2731 0.21054
>  Noccode:Fauncode:Arbcode   1    1.76  1.7566 0.0007929 0.2267 0.17210
>  Residuals                281 2177.21  7.7481                         
>  Total                    288 2215.39                                 
>                           Pr(>F)
>  Noccode                   0.420
>  Fauncode                  0.978
>  Arbcode                   0.235
>  Noccode:Fauncode          0.941
>  Noccode:Arbcode           0.539
>  Fauncode:Arbcode          0.757
>  Noccode:Fauncode:Arbcode  0.792
>  Residuals                      
>  Total

It is worth fitting these models to subsets of the data if we are concerned that certain groups of mammals (e.g., marsupials and anthropoids/primates) are not subject to the same evolutionary ‘rules’ as the other groups. More complex models could actually try to model these differences, but my goal here is simply to replicate Heesy (2008) with the addition of a phylogeny, so we’ll fit the MANOVAs using the same subsets of the data as Heesy (2008):

# drop marsupials and anthropoids, then rerun phylogenetic MANOVA
data.cc.eu1 <- data.cc[!(data.cc$ORDER%in%c("Didelphimorphia","Diprotodontia")),]
data.cc.eu1 <- data.cc.eu1[data.cc.eu1$GROUP != "Anthropoidea",]
angles3 <- data.cc.eu1[,1:3]
vars3 <- data.cc.eu1[,4:6]
tree3 <- drop.tip(tree2, setdiff(tree2$tip.label,row.names(data.cc.eu1)))
df3 <- geomorph.data.frame(angles3, vars3, tree3)
procD.pgls(angles3 ~ Noccode * Fauncode * Arbcode, phy=tree3, data=df3, print.progress=F)
>  
>  Call:
>  procD.pgls(f1 = angles3 ~ Noccode * Fauncode * Arbcode, phy = tree3,  
>      data = df3, print.progress = F) 
>  
>  
>  
>  Type I (Sequential) Sums of Squares and Cross-products
>  Randomized Residual Permutation Procedure Used
>  1000 Permutations
>  
>                            Df      SS      MS       Rsq      F       Z
>  Noccode                    1   27.70 27.6995 0.0179847 3.7741 0.93823
>  Fauncode                   1    0.46  0.4588 0.0002979 0.0625 0.01914
>  Arbcode                    1    9.26  9.2595 0.0060120 1.2616 1.00097
>  Noccode:Fauncode           1    3.49  3.4867 0.0022638 0.4751 0.27845
>  Noccode:Arbcode            1    1.25  1.2540 0.0008142 0.1709 0.12721
>  Fauncode:Arbcode           1    0.40  0.4005 0.0002600 0.0546 0.03804
>  Noccode:Fauncode:Arbcode   1    7.72  7.7250 0.0050157 1.0525 0.76869
>  Residuals                203 1489.88  7.3393                         
>  Total                    210 1540.17                                 
>                           Pr(>F)
>  Noccode                   0.261
>  Fauncode                  0.988
>  Arbcode                   0.237
>  Noccode:Fauncode          0.664
>  Noccode:Arbcode           0.833
>  Fauncode:Arbcode          0.968
>  Noccode:Fauncode:Arbcode  0.322
>  Residuals                      
>  Total
# drop marsupials and primates, then rerun phylogenetic MANOVA
data.cc.eu2 <- data.cc.eu1[data.cc.eu1$ORDER != "Primates",]
angles4 <- data.cc.eu2[,1:3]
vars4 <- data.cc.eu2[,4:6]
tree4 <- drop.tip(tree3, setdiff(tree3$tip.label,row.names(data.cc.eu2)))
df4 <- geomorph.data.frame(angles4, vars4, tree4)
procD.pgls(angles4 ~ Noccode * Fauncode * Arbcode, phy=tree4, data=df4, print.progress=F)
>  
>  Call:
>  procD.pgls(f1 = angles4 ~ Noccode * Fauncode * Arbcode, phy = tree4,  
>      data = df4, print.progress = F) 
>  
>  
>  
>  Type I (Sequential) Sums of Squares and Cross-products
>  Randomized Residual Permutation Procedure Used
>  1000 Permutations
>  
>                            Df      SS     MS       Rsq      F       Z
>  Noccode                    1   33.24 33.238 0.0236845 4.1180 1.01256
>  Fauncode                   1    0.20  0.201 0.0001429 0.0248 0.00751
>  Arbcode                    1   10.88 10.885 0.0077561 1.3486 1.10420
>  Noccode:Fauncode           1    9.17  9.172 0.0065358 1.1364 0.71120
>  Noccode:Arbcode            1    5.35  5.346 0.0038095 0.6624 0.45790
>  Fauncode:Arbcode           1    0.71  0.715 0.0005094 0.0886 0.07315
>  Noccode:Fauncode:Arbcode   1    3.97  3.967 0.0028267 0.4915 0.40791
>  Residuals                166 1339.85  8.071                         
>  Total                    173 1403.38                                
>                           Pr(>F)
>  Noccode                   0.214
>  Fauncode                  0.998
>  Arbcode                   0.191
>  Noccode:Fauncode          0.347
>  Noccode:Arbcode           0.512
>  Fauncode:Arbcode          0.914
>  Noccode:Fauncode:Arbcode  0.565
>  Residuals                      
>  Total

Wow, nothing is significant! Thus, simply incorporating phylogeny into the analysis has transformed a study for which nearly every statistical test yielded a significant result into a study with NO significant results. It is hard to find a more perfect example of how important it can be to include phylogeny in analyses of interspecific data.

Concluding remarks

It is important to keep in mind what ‘negative’ or ‘null’ results such as these do not tell us… specifically, the fact that there are no significant relationships between the ecological variables and orbit orientation does not imply that evolutionary shifts in orbit orientation were not driven by the ecological variables in question. It could be that any or all of these variables played a role in directing the evolution of early primate orbits. However, unless our data captures numerous instances of correlated shifts in orbit orientation and ecology, we are unlikely to find statistical associations between these variables after controlling for phylogenetic autocorrelation. Inferences about such singular or very rare events in evolutionary history may have to depend more on circumstantial rather than statistical evidence.

That said, I think there is room to do more phylogenetic comparative work on the evolution of orbit orientation, and I see at least two major directions for advancing this project.

These two directions are complementary, since adding more data (particularly fossils) increases statistical power to detect deviations from Brownian evolution.

References:

Back to blog