In response to my post yesterday showing how
to create a customized plot of a fitted M*k* model with *phytools* one user
asked:

*“Very cool! can we scale the arrow width according to the transition rate?”*

This is not yet an option in the plot method, per se… but, in R, it's usually not a question of 'if' but 'how'?

Let's see.

Going back to our example from yesterday of the
evolution of ecomorph state in *Anolis* lizards, let's re-load our phylogeny and dataset, as follows:

```
library(phytools)
data(anoletree)
anolis.tree<-as.phylo(anoletree)
ecomorph<-as.factor(getStates(anoletree,"tips"))
```

Now we can proceed to fit our model. We'll focus on a slightly *easier* case by imagining that we're fitting a symmetric
(`"SYM"`

) transition model.

```
fitSYM<-fitMk(anolis.tree,ecomorph,model="SYM",
pi="fitzjohn")
print(fitSYM,digits=4)
```

```
## Object of class "fitMk".
##
## Fitted (or set) value of Q:
## CG GB TC TG Tr Tw
## CG -0.0652 0.0000 0.0318 0.0000 0.0000 0.0334
## GB 0.0000 -0.0835 0.0000 0.0326 0.0000 0.0508
## TC 0.0318 0.0000 -0.1223 0.0000 0.0156 0.0750
## TG 0.0000 0.0326 0.0000 -0.0736 0.0000 0.0409
## Tr 0.0000 0.0000 0.0156 0.0000 -0.0391 0.0235
## Tw 0.0334 0.0508 0.0750 0.0409 0.0235 -0.2237
##
## Fitted (or set) value of pi:
## CG GB TC TG Tr Tw
## 0.0011 0.0015 0.0079 0.0020 0.0002 0.9874
## due to treating the root prior as (a) nuisance.
##
## Log-likelihood: -72.0732
##
## Optimization method used was "nlminb"
```

Now, let's plot our model in the ordinary way - however, we'll save the object that's returned invisibly by our plotting method:

```
par(mar=rep(0.1,4))
obj<-plot(fitSYM,show.zeros=FALSE)
```

Now let's look at our object:

```
obj
```

```
## states x y
## 1 CG 1.0 0.000000e+00
## 2 GB 0.5 8.660254e-01
## 3 TC -0.5 8.660254e-01
## 4 TG -1.0 1.224606e-16
## 5 Tr -0.5 -8.660254e-01
## 6 Tw 0.5 -8.660254e-01
```

We can see that it consists of the coordinates of all the vertices of our plot.

To start, let's draw lines between each vertex - but this time adjusting for the line with depending on the transition rate.
We'll also only draw a line if the transition rate is *larger* than a particular, small tolerance value - let's say 1e-4.

As such, we can begin by figuring out our desired line width. Let's say that our line width should range from 1 to 10. (Some R graphical devices support non-integer line widths - others don't, so what we see will depend a little bit on which device we use to render our plot!)

```
Qhat<-as.Qmatrix(fitSYM)
Qhat
```

```
## Estimated Q matrix:
## CG GB TC TG Tr Tw
## CG -0.06524880 0.00000000 0.03181507 0.00000000 0.00000000 0.03343374
## GB 0.00000000 -0.08346690 0.00000000 0.03262396 0.00000000 0.05084294
## TC 0.03181507 0.00000000 -0.12233268 0.00000000 0.01555222 0.07496539
## TG 0.00000000 0.03262396 0.00000000 -0.07356945 0.00000000 0.04094548
## Tr 0.00000000 0.00000000 0.01555222 0.00000000 -0.03909540 0.02354318
## Tw 0.03343374 0.05084294 0.07496539 0.04094548 0.02354318 -0.22373073
```

```
tol<-1e-4
q.range<-range(Qhat[Qhat>tol])
lwd.range<-c(1,10)
Q.lwd<-matrix(0,nrow(Qhat),ncol(Qhat),
dimnames=dimnames(Qhat))
for(i in 1:nrow(Q.lwd)) for(j in 1:ncol(Q.lwd)){
if(Qhat[i,j]>tol){
Q.lwd[i,j]<-(Qhat[i,j]-q.range[1])/diff(q.range)*
diff(lwd.range)+lwd.range[1]
}
}
Q.lwd
```

```
## CG GB TC TG Tr Tw
## CG 0.000000 0.000000 3.463522 0.000000 0.000000 3.708720
## GB 0.000000 0.000000 0.000000 3.586054 0.000000 6.345893
## TC 3.463522 0.000000 0.000000 0.000000 1.000000 10.000000
## TG 0.000000 3.586054 0.000000 0.000000 0.000000 4.846611
## Tr 0.000000 0.000000 1.000000 0.000000 0.000000 2.210484
## Tw 3.708720 6.345893 10.000000 4.846611 2.210484 0.000000
```

Since our matrix is *symmetric* (as per our fitted model), we only have to traverse the upper triangular half of our matrix.

Just so the plot doesn't look really dumb, I'm going to use *plotrix* to overlay circles with our character states at all the
vertices of our graph.

```
library(plotrix)
par(mar=rep(0.1,4))
plot(NA,xlim=c(-1.2,1.2),ylim=c(-1.2,1.2),
axes=FALSE,bty="n",xlab="",ylab="")
for(i in 1:nrow(Qhat)){
for(j in i:ncol(Qhat)){
if(Q.lwd[i,j]>0){
segments(obj$x[i],obj$y[i],obj$x[j],obj$y[j],
lwd=Q.lwd[i,j])
}
}
}
invisible(mapply(draw.circle,obj$x,obj$y,col="white",
MoreArgs=list(radius=0.1)))
text(obj$x,obj$y,obj$states,cex=1.3,font=2)
```

OK, that's kind of close. The problem we still have is that our lines are not arrows - and they go *all* the way to each
vertex of the graph, instead of leaving space for our character state labels.

To figure out where to draw the lines so that they end where we want them, we're going to need to use a bit of trigonometry!

Let's do that & try again:

```
par(mar=rep(0.1,4))
plot(NA,xlim=c(-1.2,1.2),ylim=c(-1.2,1.2),
axes=FALSE,bty="n",xlab="",ylab="")
for(i in 1:nrow(Qhat)){
for(j in i:ncol(Qhat)){
if(Q.lwd[i,j]>0){
slope<-abs(diff(obj$y[c(i,j)])/diff(obj$x[c(i,j)]))
x0<-obj$x[i]+0.1*cos(atan(slope))*sign(diff(obj$x[c(i,j)]))
x1<-obj$x[j]+0.1*cos(atan(slope))*sign(diff(obj$x[c(j,i)]))
y0<-obj$y[i]+0.1*sin(atan(slope))*sign(diff(obj$y[c(i,j)]))
y1<-obj$y[j]+0.1*sin(atan(slope))*sign(diff(obj$y[c(j,i)]))
arrows(x0,y0,x1,y1,length=0.1,code=3)
}
}
}
text(obj$x,obj$y,obj$states,cex=1.3,font=2)
```

Cool. That worked.

Finally, let's combine the two plots:

```
par(mar=rep(0.1,4))
plot(NA,xlim=c(-1.2,1.2),ylim=c(-1.2,1.2),
axes=FALSE,bty="n",xlab="",ylab="")
for(i in 1:nrow(Qhat)){
for(j in i:ncol(Qhat)){
if(Q.lwd[i,j]>0){
slope<-abs(diff(obj$y[c(i,j)])/diff(obj$x[c(i,j)]))
x0<-obj$x[i]+0.1*cos(atan(slope))*sign(diff(obj$x[c(i,j)]))
x1<-obj$x[j]+0.1*cos(atan(slope))*sign(diff(obj$x[c(j,i)]))
y0<-obj$y[i]+0.1*sin(atan(slope))*sign(diff(obj$y[c(i,j)]))
y1<-obj$y[j]+0.1*sin(atan(slope))*sign(diff(obj$y[c(j,i)]))
arrows(x0,y0,x1,y1,length=0.1,code=3,lwd=Q.lwd[i,j])
}
}
}
text(obj$x,obj$y,obj$states,cex=1.3,font=2)
par(lend=1)
legend("topleft",
c(paste("q = ",round(q.range[1],3)),
paste("q = ",round(mean(q.range),3)),
paste("q = ",round(q.range[2],3))),
lwd=c(lwd.range[1],mean(lwd.range),lwd.range[2]),
bty="n")
```

That's OK. Phew!

Alright, for non-symmetric transition models I think it would be easier to program it into the `plot.Qmatrix`

function…. Look
for that in a future *phytools* update!

## No comments:

## Post a Comment

Note: due to the very large amount of spam, all comments are now automatically submitted for moderation.