Orbit orientation isn't correlated with anything: a reanalysis of Heesy (2008)
11 Jun 2017In this post, I reanalyze 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 nonindependence of species data undermines the paper’s major conclusions, and I highlight directions for future work.
The ecomorphology of orbit orientation (Heesy 2008)
“Forwardfacing 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 GraspLeaping 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 lowlight 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 inionnasion 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, oneway 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 nonphylogenetic analyses. As the R code below demonstrates, my nonphylogenetic 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.2e16 ***
> Fauncode 1 0.90222 10.079 3 279 2.516e06 ***
> Arbcode 1 0.76720 28.220 3 279 5.740e16 ***
> Noccode:Fauncode 1 0.90691 9.546 3 279 5.072e06 ***
> Noccode:Arbcode 1 0.90855 9.361 3 279 6.467e06 ***
> 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.2e16 ***
> Noccode:Fauncode 1 2211 2211.0 11.5183 0.0007883 ***
> Noccode:Arbcode 1 3528 3528.3 18.3807 2.489e05 ***
> 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.817e06 ***
> Fauncode 1 8092 8092.4 28.5416 1.897e07 ***
> Arbcode 1 5290 5290.5 18.6592 2.170e05 ***
> 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.2e16 ***
> 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.424e05 ***
> 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.595e15 ***
> Fauncode 1 0.86746 10.3386 3 203 2.313e06 ***
> Arbcode 1 0.83768 13.1117 3 203 7.338e08 ***
> 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.975e08 ***
> 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.458e08 ***
> 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.581e13 ***
> 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.896e10 ***
> Fauncode 1 0.79615 13.9969 3 164 3.610e08 ***
> 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.738e09 ***
> 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.846e05 ***
> 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 rundown of the results (focusing on main effects):
 Main effects are always in the expected direction (nocturnality, faunivory, and arboreality are associated with more forwardfacing orbits)
 All main effects are significant in all three MANOVAs.
 In oneway ANOVAs for the full dataset, all main effects are significant
 In oneway ANOVAs for the subset without marsupials and anthropoids, all main effects are significant except for nocturnality, which is a nonsignificant predictor of convergence and frontition. This is the biggest discrepency between my results and Heesy (2008), who found that nocturnality was significant in all of the oneway ANOVAs for this group.
 In oneway ANOVAs for the subset without marsupials and primates, only the main effects for nocturnality and insectivory are significant (arboreality is not)
This last finding that arboreality is nonsignificant in the oneway 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 nonsignificant 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 pvalues in Heesy’s (2008) original results.
Problem of phylogenetic autocorrelation
Phylogenetic autocorrelation, also called phylogenetic pseudoreplication or phylogenetic nonindependence, is a wellknown 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 BinindaEmonds 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 BinindaEmonds 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 approxF numDf denDf Pr(>F) Pr(>F) given phy
> group 1 0.6222 58.697 3 290 1.1091e29 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 approxF numDf denDf Pr(>F) Pr(>F) given phy
> group 1 0.90743 9.895 3 291 3.1196e06 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 approxF numDf denDf Pr(>F) Pr(>F) given phy
> group 1 0.75702 30.599 3 286 3.4778e17 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 (nonphylogenetic) and phylogenetic pvalues in the last two columns of the MANOVA table. We can see that all of the standard pvalues are highly significant, while the phylogenetic pvalues 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 preprocessing 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 Crossproducts
> 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 Crossproducts
> 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 Crossproducts
> 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.

First is expanding and refining the dataset. There is room for more complete sampling of extant taxa, and for the addition of fossil data. It also would be good to include other aspects of skull morphology in the models, such as the size and shape of the braincase and basicranium, perhaps in a geometric morphometrics framework. Additionally, the ecological variables could use some work. Although activity pattern is relatively straightforward to code as a discrete variable, diet and substrate use are much more fraught, and predictor variables should be selected to reflect the hypotheses being tested as directly as possible. For example, if “leaping” is a factor that is thought to select for convergent orbits, then “leaping” is a much better predictor variable than “arboreality”, which includes things like lorises and sloths that creep slowly along branches. Similarly, if “visual predation”” is a factor thought to select for convergent orbits, then “visual predation” is a much better predictor variable than “faunivory”, which includes things like ayeayes that find their prey by tapping on trees like woodpeckers.

Second is fitting more sophisticated models that aim to detect shifts in the tempo or mode of evolution (e.g., like this). This is a rather different approach that aims to identify parts of a phylogeny where the rules of evolution seem to change. For instance, it would be interesting to know if there is a statistically significant ‘shift’ in the mean orbit orientation at the base of the primate clade. That wouldn’t tell us why the shift occurred, but it would at least support our assumption that something happened there that warrants explanation. I think the Nocturnal Visual Predation hypothesis would predict a mean shift at the base of both the primates and the carnivores, as the ancestors of both groups are thought to have been nocturnal visual predators.
These two directions are complementary, since adding more data (particularly fossils) increases statistical power to detect deviations from Brownian evolution.
References:

BinindaEmonds, O.R., Cardillo, M., Jones, K.E., MacPhee, R.D., Beck, R.M., Grenyer, R., Price, S.A., Vos, R.A., Gittleman, J.L. and Purvis, A., 2007. The delayed rise of presentday mammals. Nature 446: 507512.

Cartmill, M. 1992. New views on primate origins. Evolutionary Anthropology 1(3):10511.

Heesy, C.P., 2008. Ecomorphology of orbit orientation and the adaptive significance of binocular vision in primates and other mammals. Brain, Behavior and Evolution 71:5467.

Martins, E.P., and T. Garland. 1991. Phylogenetic analyses of the correlated evolution of continuous character: a simulation study. Evolution 45:53457.

Rasmussen. D.T. 1990. Primate origins: lessons from a neotropical marsupial. American Journal of Primatology 22:26377.

Sussman, R.W. 1991. Primate origins and the evolution of angiosperms. American Journal of Primatology 23:209223.

Szalay, F.S., and M. Dagosto. 1980. Locomotor adaptations as reflected on the humerus of Paleogene primates. Folia Primatologica 34:145.