An astute reader of our book, and follower of my blog, recently
pointed out that *Ateles*
(the spider monkeys) were incorrectly coded as cathemeral in our
dataset for phylogenetic ANCOVA of this chapter, when they should’ve
properly been coded as unambiguously *diurnal*.

Honestly, I cannot account for where this error originates.
Fortunately, it does not qualitatively affect any of the analyses given
in the chapter. Nonetheless, the results are *quantitatively*
affected (in the sense that the values of coefficients and test
statistics change slightly) when this error is corrected, so I thought
it would be worthwhile to post the *corrected* results here. To
ensure that readers of the book obtain quantitatively equivalent results
to those of the published version, we will *not* replace the
posted data file; however, I have now included the corrected primates
dataset of the book chapter with *phytools* (>= 1.7-8), so
readers with that version of the package can follow along below.

```
## load nlme
library(nlme)
## load phytools
library(phytools)
## check package version
packageVersion("phytools")
```

`## [1] '1.7.8'`

```
## load primate dataset and tree
data(primate.tree)
data(primate.data)
head(primate.data[,5:6],10)
```

```
## Activity_pattern Activity_pattern_code
## Allenopithecus_nigroviridis Diurnal 0
## Alouatta_palliata Diurnal 0
## Alouatta_seniculus Diurnal 0
## Aotus_trivirgatus Nocturnal 2
## Arctocebus_aureus Nocturnal 2
## Arctocebus_calabarensis Nocturnal 2
## Ateles_fusciceps Diurnal 0
## Ateles_geoffroyi Diurnal 0
## Ateles_paniscus Diurnal 0
## Avahi_laniger Nocturnal 2
```

The diel activity pattern of all *Ateles* species can now be
seen to be correctly coded as `Diurnal`

rather than
`Cathemeral`

.

Now let’s re-run our analyses.

```
## create correlation structure
spp<-rownames(primate.data)
corBM<-corBrownian(phy=primate.tree,form=~spp)
## fit ANCOVA model
primate.ancova<-gls(log(Orbit_area)~log(Skull_length)+
Activity_pattern,data=primate.data,
correlation=corBM)
anova(primate.ancova)
```

```
## Denom. DF: 86
## numDF F-value p-value
## (Intercept) 1 2045.7504 <.0001
## log(Skull_length) 1 382.0018 <.0001
## Activity_pattern 2 9.8960 1e-04
```

The *F*-values and *p*-values of our fitted model
coefficients change slightly, but not in a substantive way compared to
the analysis that was published in the book.

```
## set the margins of our plot using par
par(mar=c(5.1,5.1,2.1,2.1))
## set the point colors for the different levels
## of our factor
pt.cols<-setNames(c("#87CEEB","#FAC358","black"),
levels(primate.data$Activity_pattern))
## plot the data
plot(Orbit_area~Skull_length,data=primate.data,pch=21,
bg=pt.cols[primate.data$Activity_pattern],
log="xy",bty="n",xlab="skull length (cm)",
ylab=expression(paste("orbit area (",mm^2,")")),
cex=1.2,cex.axis=0.7,cex.lab=0.8)
## add a legend
legend("bottomright",names(pt.cols),pch=21,pt.cex=1.2,
pt.bg=pt.cols,cex=0.8)
## create a common set of x values to plot our
## different lines for each level of the factor
xx<-seq(min(primate.data$Skull_length),
max(primate.data$Skull_length),length.out=100)
## add lines for each level of the factor
lines(xx,exp(predict(primate.ancova,
newdata=data.frame(Skull_length=xx,
Activity_pattern=as.factor(rep("Cathemeral",100))))),
lwd=2,col=pt.cols["Cathemeral"])
lines(xx,exp(predict(primate.ancova,
newdata=data.frame(Skull_length=xx,
Activity_pattern=as.factor(rep("Diurnal",100))))),
lwd=2,col=pt.cols["Diurnal"])
lines(xx,exp(predict(primate.ancova,
newdata=data.frame(Skull_length=xx,
Activity_pattern=as.factor(rep("Nocturnal",100))))),
lwd=2,col=pt.cols["Nocturnal"])
```