---
title: "Searching for temporal trends in phenotypes or rates"
author: "Silvia Castiglione, Carmela Serio, Pasquale Raia"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{search.trend}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

require(RRphylo)
options(knitr.kable.NA = '',rmarkdown.html_vignette.check_title = FALSE)
```

## Index
1. [search.trend basics](#basics)
2. [Temporal trends on the entire tree](#tree)
3. [Temporal trends at clade level](#nodes)
4. [Multivariate data](#multi)
5. [Guided examples](#examples)

## search.trend basics{#basics}
The branch-wise estimation of phenotypic evolutionary rates and the computation of ancestral states at each node make [RRphylo](RRphylo.html) especially suitable to study temporal trends in phenotypic means and absolute rates applying to the entire phylogeny or just on a part of it. This is particularly true as the RRphylo method is specifically meant to work with phylogenies of extinct species, taking full advantage from fossil information. The function `search.trend` (Castiglione et al. 2019) is designed to explore the domain of macroevolution, addressing patterns such as Cope's rule, or to test other specific research questions as to whether the rate of evolution increased/decreased in a particular tree or clade within it. Although possible in principle, running `search.trend` without extinct species in the tree makes little sense.   

`search.trend` takes an object produced by `RRphylo`, which is necessary to retrieve branch-wise rate estimates and ancestral characters at nodes. By default, the function searches for significant temporal trends in phenotypic means and absolute evolutionary rates as applying to the entire phylogeny. If a specific hypothesis about one or more clades is tested, the clades presumed to experience trended evolution, in either phenotypes, rates or both, must be indicated by the argument `node`. If more than one clade is under testing, the function performs a pairwise comparison of trends also between the clades.

## Temporal trends on the entire tree{#tree}
To search for temporal trends in either phenotypes or absolute rates occurring on the entire phylogeny, the function computes the regression slopes between phenotypic (absolute rate) values at nodes and tips and their age, meant as the distance from the tree root. Since a significant relationship between phenotypes (absolute rates) and age is possible even under the Brownian Motion (BM from now on), significance is assessed comparing the real slope to a family of 100 slopes (the number can be specified by setting the argument `nsim` within the function) generated by simulating phenotypes according to the BM process (*BMslopes*). 

```{r,echo=FALSE,message=FALSE,warning=FALSE,fig.dim=c(8,4),fig.align='center',out.width='100%',dpi=220}
require(ape)
require(geiger)
require(phytools)
require(ddpcr)

RRphylo:::range01->range01

sim.bdtree(b = .5, d = 0.2,seed=14)->tree
rt<-0
s2m<-1
es=2

set.seed(14)
fastBM(tree,a=rt,sig2=s2m)->y

### TREND ###
yt <- (diag(vcv(tree))^es)/(diag(vcv(tree))) * y
data.frame(tree$edge[,2],nodeHeights(tree)[,2])->hei
tree$tip.label[hei[which(hei[,1]<(Ntip(tree)+1)),1]]->hei[which(hei[,1]<(Ntip(tree)+1)),1]
hei[,2]->rootD
c(0.0001,rootD)->rootD
names(rootD)<-c((Ntip(tree)+1),hei[,1])

makeL(tree)->L
makeL1(tree)->L1
RRphylo(tree,yt)->rr

quiet(search.trend(rr,yt,clus=2/parallel::detectCores(),foldername=tempdir())->st.rates)

c(rr$aces,yt)->phen

rr$rates->rts
abs(rts)->rts
as.matrix(as.data.frame(rts[match(names(rootD),rownames(rts)),]))->rts
range01(rts)->rts
rootD->rootC
tend<-2

rr$ace[1]->a
rootV<-a
rr$lambda->lambda
BTS<-list()
for(i in 1:100){
  rootC->rootB
  fastBM(tree,sig2=1,a=a,bounds=c(min(yt),max(yt)))->yb
  betas <- (solve(t(L) %*% L + lambda * diag(ncol(L))) %*% 
              t(L)) %*% (as.matrix(yb)-rootV)
  bts <- abs(betas)
  RRphylo:::range01(bts)->BTS[[i]]
  as.matrix(as.data.frame(bts[match(names(rootB),rownames(bts)),]))->bts
  bts
}

sapply(BTS,function(x) x[match(names(rootD),rownames(x)),])->fu

apply(fu,1,function(x) quantile(x,.975))->ful
apply(fu,1,function(x) quantile(x,.025))->mins
predict(loess(ful~rootD))->fk
predict(loess(mins~rootD))->fk.min
names(fk)<-names(ful)
names(fk.min)<-names(ful)

##abline(lm(rts~rootD))
data.frame(rootD,fk,fk.min)->dfk
dfk[order(dfk[,1]),]->dfk

data.frame(rts=rts,col=rep("red",nrow(rts)))->rts
as.character(rts[,2])->rts[,2]
rts[which(rownames(rts)%in%tree$tip.label),2]<-"green"

par(mfrow=c(1,2),mar=c(4,3,3,1))
plot(NA,ylim=range(rts[,1]),xlim=range(rootD),
     mgp=c(1.8,0.5,0),
     xlab="distance from the root",ylab="rescaled rates",main="Simulation for trend\nin absolute rates")
polygon(c(dfk$rootD,rev(dfk$rootD)),c(dfk$fk,rev(dfk$fk.min)),col = rgb(0.5, 0.5, 0.5,
                                                                        0.3), border = NA)
points(rootD,rts[,1],cex=1.5,bg=rts$col,
       col="black",pch=21)
abline(lm(rts[,1]~rootD),lwd=2,col="#ff00ff")

legend("topleft",legend=c("nodes","tips","brownian range"),fill=c("red","green","#dadad9"),bty="n")

#### DRIFT
ds=1
yt1<-y+diag(vcv(tree))*ds

RRphylo:::range01(yt1)->yd


data.frame(tree$edge[,2],nodeHeights(tree)[,2])->hei
tree$tip.label[hei[which(hei[,1]<(Ntip(tree)+1)),1]]->hei[which(hei[,1]<(Ntip(tree)+1)),1]
hei[,2]->rootD
c(0.0001,rootD)->rootD
names(rootD)<-c((Ntip(tree)+1),hei[,1])

makeL(tree)->L
makeL1(tree)->L1
RRphylo(tree,yd)->rr
quiet(search.trend(rr,yd,clus=2/parallel::detectCores(),foldername=tempdir())->st.phen)

c(rr$aces,yd)->phen
names(phen)<-c(rownames(rr$aces),names(yd))
rootD[match(names(phen),names(rootD))]->rootD

rootD->rootP

phen[[1]]->a->rootV
lm(phen~rootD)->regr

#ratematrix(tree,yt1)[1]->S2[m] 
phenD<-list()
yl<-list()
for(i in 1:100)
{
  fastBM(tree,sig2=s2m,a=rootV,bounds=c(min(yd),max(yd)))->yc
  yc->yl[[i]]
  RRphylo:::range01(yc)->yc
  lambda <- rr$lambda
  betas <- (solve(t(L) %*% L + lambda * diag(ncol(L))) %*% 
              t(L)) %*% (as.matrix(yc) - rootV)
  aceRR <- (L1 %*% betas[1:Nnode(tree), ]) + rootV
  c(aceRR,yc)->phenC
  names(phenC)<-c(rownames(aceRR),names(yc))
  phenC->phenD[[i]]
  #coef(summary(lm(phenC~rootP)))[2]->slopec[i]
  
}

sapply(phenD,function(x) x[match(names(rootD),names(x))])->fu

apply(fu,1,function(x) quantile(x,.975))->ful
apply(fu,1,function(x) quantile(x,.025))->mins
predict(loess(ful~rootD))->fk
predict(loess(mins~rootD))->fk.min
names(fk)<-names(ful)
names(fk.min)<-names(ful)

data.frame(rootD,fk,fk.min)->dfk
dfk[order(dfk[,1]),]->dfk

data.frame(phen=phen,col=rep("red",length(phen)))->phen
as.character(phen[,2])->phen[,2]
phen[which(rownames(phen)%in%tree$tip.label),2]<-"green"

plot(NA,ylim=range(phen[,1]),xlim=range(rootD),
     mgp=c(1.8,0.5,0),
     xlab="distance from the root",ylab="rescaled phenotypes",main="Simulation for trend\nin mean phenotypes")
polygon(c(dfk$rootD,rev(dfk$rootD)),c(dfk$fk,rev(dfk$fk.min)),col = rgb(0.5, 0.5, 0.5,
                                                                        0.3), border = NA)
points(rootD,phen[,1],cex=1.5,bg=phen$col,
       col="black",pch=21)
abline(lm(phen[,1]~rootD),lwd=2,col="#ff00ff")



as.data.frame(rbind(c(st.rates[[4]],NA),c(st.phen[[3]][1:3],NA,st.phen[[3]][4])))->res

colnames(res)[4:5]<-c("spread","dev")
rownames(res)<-c("absolute rate regression","phenotypic regression")
```
```{r,eval=FALSE,message=FALSE,warning=FALSE}
search.trend(RR=RR,y=y,foldername=getwd())->ST
```
```{r,echo=FALSE,message=FALSE,warning=FALSE}
require(kableExtra)

knitr::kable(res,digits=3,align="c") %>%
kable_styling(full_width = F, position = "center")  

```

As in the table above, `search.trend` results for phenotypic (`$phenotypic.regression`) and absolute rate (`$rate.regression`) regressions both include: the actual regression slope and its p-value (**p.real**), and the p-value derived by contrasting the real slope to the *BMslopes* (**p.random**). The user wants to look at **p.random** to see if a real trend applies. 

The output of `search.trend` also includes two metrics specifically designed to quantify the magnitude of the phenotypic/absolute rates deviation from BM. The **spread** represents the deviation from BM produced by a trend in absolute rates. It is calculated as the ratio between the range of phenotypic values and the range of such values halfway along the tree height, divided to the same metric value generated under BM. **spread** is 1 under BM.
The **dev** quantifies the deviation from BM produced by a trend in phenotypic means. It represents the deviation of the phenotypic mean from the root value in terms of the number of standard deviations of the trait distribution. **dev** is zero under the BM.

Therefore, the example above produced significant results for both phenotypic and absolute rates temporal trends (**p.random** are both < 0.05). Absolute rates increase over time is significant (**p.real** < 0.05) and twice as much (**dev** = 1.941) as expected under BM. The increase in phenotypic means through time is significant (**p.real** < 0.05) and significantly different from the BM simulation, describing a deviation in mean phenotype some one half of the standard deviation of the phenotype (**spread** = 0.542).


## Temporal trends at clade level{#nodes}
A peculiarity of `search.trend` is that it can test whether individual clades follow a different trend in phenotypic means (absolute rates) over time, as compared to the rest of the tree. 

In the case of phenotypic trend, individual clades are tested for the hypothesis that trend intensity (slope) does not depart from BM and whether the marginal means of the regression values differ from the means calculated for the rest of the tree. In the case of a trend in absolute rates the actual regression slope depends on the relative position (age) of the focal nodes respective to the root, given the exponential nature of phenotypic variance change in time. Because of this, estimated marginal means and regression slopes for the rate versus age regression of the focal clade are contrasted directly to the same metrics generated for the rest of the tree. The comparisons for both phenotypes (estimated marginal means) and rates regression (estimated marignal means and regression slopes) are implemented by using the functions `emmeans` and `emtrends`  in package emmeans (Lenth 2020).

If more than one clade is under testing, the function performs the pairwise comparisons of trend intensity between clades. For absolute rates regression, such comparison is performed by computing differences in estimated marginal means and regression slopes as described above. Trends in phenotypic means are compared by computing the difference in regression slopes between the clades and contrasting such value to a random distribution of differences generated under BM. The difference in estimated marginal means of phenotype versus age regression between clades is also returned. 

```{r,echo=FALSE,message=FALSE,warning=FALSE,fig.dim=c(8,4),fig.align='center',out.width='100%',dpi=220}
require(plotrix)

sim.bdtree(b = .5, d = 0.2,seed=14)->T

set.seed(93)
sizedsubtree(T,20)->n
repeat({
  c(getMommy(T,n),getDescendants(T,n))->outs
  sizedsubtree(T,20)->n2
  if(n2%in%outs==FALSE & getSis(T,n2,printZoom = FALSE)!=n & n!=n2) break
})

max(nodeHeights(T))/nodeHeights(T)[,2][which(T$edge[,2]==n)]->dN1
max(nodeHeights(T))/nodeHeights(T)[,2][which(T$edge[,2]==n2)]->dN2
ds1<-.5
ds2<- -.5
es1=.01
es2=1.9
runif(1,-10,10)->rt
S2<-sample(seq(0.1,10,0.01),1)
sample(c(es1,es2),1)->esx1
sample(c(es1,es2),1)->esx2

fastBM(T,a=rt,sig2=S2)->yB
yB->y

### Phenotypes ###
y[match(tips(T,n),names(y))]->yDR
diag(vcv(T)[names(yDR),names(yDR)])->timeDR
coef(lm(yDR~timeDR))[2]->tend
if(tend*ds1>0) ds1->ds1 else ds1*(-1)->ds1
ds1*dN1->dsx1
yD <- (timeDR*dsx1)+ yDR
y[match(tips(T,n),names(y))]<-yD

y[match(tips(T,n2),names(y))]->yDR
diag(vcv(T)[names(yDR),names(yDR)])->timeDR
coef(lm(yDR~timeDR))[2]->tend
if(tend*ds2>0) ds2->ds2 else ds2*(-1)->ds2
ds2*dN2->dsx2
yD <- (timeDR*dsx2)+ yDR
y[match(tips(T,n2),names(y))]<-yD

RRphylo(T,y)->RR
quiet(search.trend(RR,y,node=c(n,n2),clus=2/parallel::detectCores(),foldername=tempdir())->STphen)
c(RR$aces[,1],y)->phen

c(n,getDescendants(T,n))->desn
T$tip.label[desn[which(desn<=Ntip(T))]]->desn[which(desn<=Ntip(T))]

c(n2,getDescendants(T,n2))->desn2
T$tip.label[desn2[which(desn2<=Ntip(T))]]->desn2[which(desn2<=Ntip(T))]

dist.nodes(T)[Ntip(T)+1,]->ages
T$tip.label[as.numeric(names(ages)[which(as.numeric(names(ages))<=Ntip(T))])]->names(ages)[which(as.numeric(names(ages))<=Ntip(T))]
ages[match(names(phen),names(ages))]->ages


### Rates ###
es1=.01
es2=1.9
sample(c(es1,es2),1)->esx1
sample(c(es1,es2),1)->esx2

yB->y
y[match(tips(T,n),names(y))]->yT
diag(vcv(T)[names(yT),names(yT)])->timeT
yT <- ((timeT^esx1)/timeT)*yT

y[match(tips(T,n2),names(y))]->yT2
diag(vcv(T)[names(yT2),names(yT2)])->timeT2
yT2 <- ((timeT2^esx2)/timeT2)*yT2
y[match(tips(T,n),names(y))]<-yT
y[match(tips(T,n2),names(y))]<-yT2

RRphylo(T,y)->RR
quiet(search.trend(RR,y,node=c(n,n2),clus=2/parallel::detectCores(),foldername=tempdir())->STrates)
abs(RR$rates[,1])->rats


par(mfrow=c(1,2),mar=c(4,3,3,1))
plot(max(ages)-ages[-which(names(ages)%in%c(desn,desn2))],rats[-which(names(rats)%in%c(desn,desn2))],
     ylim=range(rats),xlim=c(max(ages),min(ages)),
     pch=21,bg="gray80",cex=1.5,mgp=c(1.8,0.5,0),
     xlab="age",ylab="absolute rates",main="Simulation for trend\nin absolute rates")
abline(lm(rats~ages),lwd=3,col="gray20")
points(max(ages)-ages[which(names(ages)%in%desn)],rats[which(names(rats)%in%desn)],
       xlim=c(max(ages),min(ages)),
       pch=21,bg="slateblue2",cex=1.5)
ablineclip(lm(rats[which(names(rats)%in%desn)]~I(max(ages)-ages[which(names(ages)%in%desn)])),
           x1=max(max(ages)-ages[which(names(ages)%in%desn)]),
           x2=min(max(ages)-ages[which(names(ages)%in%desn)]),
           col="black",lwd=4.5)
ablineclip(lm(rats[which(names(rats)%in%desn)]~I(max(ages)-ages[which(names(ages)%in%desn)])),
           x1=max(max(ages)-ages[which(names(ages)%in%desn)]),
           x2=min(max(ages)-ages[which(names(ages)%in%desn)]),
           col="slateblue2",lwd=3)
points(max(ages)-ages[which(names(ages)%in%desn2)],rats[which(names(rats)%in%desn2)],
       xlim=c(max(ages),min(ages)),
       pch=21,bg="#ae9eff",cex=1.5)
ablineclip(lm(rats[which(names(rats)%in%desn2)]~I(max(ages)-ages[which(names(ages)%in%desn2)])),
           x1=max(max(ages)-ages[which(names(ages)%in%desn2)]),
           x2=min(max(ages)-ages[which(names(ages)%in%desn2)]),
           col="black",lwd=4.5)
ablineclip(lm(rats[which(names(rats)%in%desn2)]~I(max(ages)-ages[which(names(ages)%in%desn2)])),
           x1=max(max(ages)-ages[which(names(ages)%in%desn2)]),
           x2=min(max(ages)-ages[which(names(ages)%in%desn2)]),
           col="#ae9eff",lwd=3)
legend("topleft",legend=c("node 275","node 198","entire tree"),fill=c("slateblue2","#ae9eff","gray20"),bty="n",x.intersp = .5)

plot(max(ages)-ages[-which(names(ages)%in%c(desn,desn2))],phen[-which(names(phen)%in%c(desn,desn2))],
     ylim=range(phen),xlim=c(max(ages),min(ages)),
     pch=21,bg="gray80",cex=1.5,mgp=c(1.8,0.5,0),
     xlab="age",ylab="phenotype",main="Simulation for trend\nin mean phenotypes")
abline(lm(phen~ages),lwd=3,col="gray20")
points(max(ages)-ages[which(names(ages)%in%desn)],phen[which(names(phen)%in%desn)],
       xlim=c(max(ages),min(ages)),
       pch=21,bg="aquamarine3",cex=1.5)
ablineclip(lm(phen[which(names(phen)%in%desn)]~I(max(ages)-ages[which(names(ages)%in%desn)])),
           x1=max(max(ages)-ages[which(names(ages)%in%desn)]),
           x2=min(max(ages)-ages[which(names(ages)%in%desn)]),
           col="black",lwd=4.5)
ablineclip(lm(phen[which(names(phen)%in%desn)]~I(max(ages)-ages[which(names(ages)%in%desn)])),
           x1=max(max(ages)-ages[which(names(ages)%in%desn)]),
           x2=min(max(ages)-ages[which(names(ages)%in%desn)]),
           col="aquamarine3",lwd=3)

points(max(ages)-ages[which(names(ages)%in%desn2)],phen[which(names(phen)%in%desn2)],
       xlim=c(max(ages),min(ages)),
       pch=21,bg="aquamarine",cex=1.5)
ablineclip(lm(phen[which(names(phen)%in%desn2)]~I(max(ages)-ages[which(names(ages)%in%desn2)])),
           x1=max(max(ages)-ages[which(names(ages)%in%desn2)]),
           x2=min(max(ages)-ages[which(names(ages)%in%desn2)]),
           col="black",lwd=4.5)
ablineclip(lm(phen[which(names(phen)%in%desn2)]~I(max(ages)-ages[which(names(ages)%in%desn2)])),
           x1=max(max(ages)-ages[which(names(ages)%in%desn2)]),
           x2=min(max(ages)-ages[which(names(ages)%in%desn2)]),
           col="aquamarine",lwd=3)
legend("topleft",legend=c("node 275","node 198","entire tree"),fill=c("aquamarine3","aquamarine","gray20"),bty="n",x.intersp = .5)

cbind(data.frame(node=names(STrates[[6]]),do.call(rbind,STrates[[6]])),
      data.frame(node=names(STphen[[5]]),do.call(rbind,STphen[[5]])))->res
```
```{r,eval=FALSE,message=FALSE,warning=FALSE}
search.trend(RR=RR,y=y,node=c(275,198),foldername=getwd())->ST
```
```{r,echo=FALSE,message=FALSE,warning=FALSE}
rownames(res)<-NULL
knitr::kable(res,digits=3,align="c") %>%
kable_styling(full_width = F, position = "center") %>%
column_spec(5, border_right = TRUE) %>%
add_header_above(c("Trend in absolute rates" = 5, "Trend in phenotypic means" = 5))

```

Same as for the entire tree, `search.trend` returns the results for phenotypic (`$node.phenotypic.regression`) and absolute rate (`$node.rate.regression`) regressions over time at each node under testing. For trends in absolute rates, significance level for both differences in estimated marginal means (**emm.difference** and **p.emm**) and regression slopes (**slope.difference** and **p.slope**) between each clade and the rest of the tree are reported. Results for trend in phenotypic means through each node include significance for the comparison of the real slope to the random slopes (**slope** and **p.slope**), and the significance for the difference in estimated marginal means between each clade and the rest of the tree (**emm.difference** and **p.emm**).

In the example above, the absolute rates versus age regression through the clade subtended by node 275 (just node 275 from now on) produced no significant difference in the estimated marginal means of the clade as compared to the rest of the tree (**p.emm** > 0.05) but a significant and negative difference in regression slopes (**p.slope** < 0.05). Conversely, the estimated marginal means of regression through node 198 are significantly lower than the rest of the tree (**emm.difference** < 0 and **p.emm** < 0.05), while the difference in slopes is not significant (**p.slope** > 0.05). As for the trend in phenotypic means, this is positive and significantly different from BM for node 275 (**p.slope** < 0.05). Such clade shows significantly larger estimated marginal means of phenotypes versus age regression than the rest of the tree (**emm.difference** > 0 and **p.emm** < 0.05). The phenotypic trend through node 198 is not significantly different from BM (**p.slope** > 0.05), but the estimated marginal means are larger than for the rest of the tree (**p.emm** < 0.05). 

Please note, the dark gray line in both figures represents the whole tree, inclusive of both clades under testing.

```{r,echo=FALSE,message=FALSE,warning=FALSE}
knitr::kable(STrates[[7]][[2]],digits=3,align="c") %>%
kable_styling(full_width = F, position = "center") %>%
add_header_above(c("Comparison of trends in absolute rates" = 6))

knitr::kable(STphen[[7]][[1]],digits=3,align="c") %>%
kable_styling(full_width = F, position = "center") %>%
add_header_above(c("Comparison of trends in phenotypic means" = 6))

```

If more than one clade is indicated, the results also include pairwise comparison of trends between clades (`$group.comparison`). The comparison of trends in absolute rates is performed in the same way as the comparison of a single clade against the rest of the tree. Thus, results consist of differences in both estimated marginal means (**emm.difference**) and slopes (**slope.difference**) between clades with respective p.values (**p.emm** and **p.slope** respectively). For the comparison of phenotypic trends between clades, the differences in regression slopes (**slope.difference**) and estimated marginal means (**emm.difference**) between clades are returned. In this case the p-value for the slope difference (**p.slope**) is computed through randomization.

In the example above, the comparison of trends in absolute rates at nodes 275 and 198 is marginally significant for slope difference only (**p.emm** > 0.05 and **p.slope** < 0.05). Conversely, trends in phenotypic means are significantly different in both slope (**p.slope** < 0.05) and estimate marginal means (**p.emm** < 0.05).

The function stores plots of phenotype versus age and absolute rate versus age regression (similar to those in the figures above) as pdf files within the folder indicated by the argument `foldername`.

## Multivariate data {#multi}
When applied on multivariate data, `search.trend` treats each phenotypic and rate (derived by `RRphylo`) component independently. Also, it performs the whole set of analyses on a multivariate vector of phenotypes (rates), derived by computing the L2-norm of the individual phenotypes (rates) at each branch (the same as in multivariate RRphylo).

## Guided examples {#examples}

```{r,out.width='100%',fig.dim=c(8,7),message=FALSE,dpi=220}
# load the RRphylo example dataset including Cetaceans tree and data
data("DataCetaceans")
DataCetaceans$treecet->treecet # phylogenetic tree
DataCetaceans$masscet->masscet # logged body mass data
DataCetaceans$brainmasscet->brainmasscet # logged brain mass data
DataCetaceans$aceMyst->aceMyst # known phenotypic value for the most recent common ancestor of Mysticeti

require(geiger)
par(mar=c(0,0,0,1))
plot(ladderize(treecet,FALSE),show.tip.label = FALSE,edge.color = "gray40")
plotinfo<-get("last_plot.phylo",envir =ape::.PlotPhyloEnv)
nodelabels(text="",node=128,frame="circle",bg="red",cex=0.5)
nodelabels(text="Mystacodon",node=128,frame="n",bg="w",cex=1,adj=c(-0.1,0.5),font=2)
range(plotinfo$yy[which(treecet$tip.label%in%tips(treecet,128))])->yran128
rect(plotinfo$x.lim[2]+0.4,yran128[1],plotinfo$x.lim[2]+0.7,yran128[2],col="red",border="red")
range(plotinfo$yy[which(treecet$tip.label%in%tips(treecet,142))])->yran142
rect(plotinfo$x.lim[2]+0.4,yran142[1],plotinfo$x.lim[2]+0.7,yran142[2],col="blue",border="blue")
mtext(c("Mysticeti","Odontoceti"), side = 4,line=-0.5,at=c(sum(yran128)/2,sum(yran142)/2),
      cex=1.5,adj=0.5,col=c("red","blue"))
```

```{r,eval=FALSE,message=FALSE,dpi=220}
# check the order of your data: best if data vectors
# are sorted in the same order of the species on the phylogeny
masscet[match(treecet$tip.label,names(masscet))]->masscet

# Set the body mass of Mysticetes ancestor (Mystacodon selenensis) 
# as known value at node and perform RRphylo on the vector of (log) body mass
RRphylo(tree=treecet,y=masscet,aces=aceMyst)->RR

# Perform search.trend on the RR object and (log) body mass by indicating Mysticeti as focal clade
search.trend(RR=RR,y=masscet,node=as.numeric(names(aceMyst)),foldername=tempdir())

```


## References
Castiglione, S., Serio, C., Mondanaro, A., Di Febbraro, M., Profico, A., Girardi, G., & Raia, P. (2019). Simultaneous detection of macroevolutionary patterns in phenotypic means and rate of change with and within phylogenetic trees including extinct species. PloS one, 14: e0210101. doi.org/10.1371/journal.pone.0210101

Lenth, R. (2020). emmeans: Estimated Marginal Means, aka Least-Squares Means. R package version 1.4.4. https://CRAN.R-project.org/package=emmeans


