Comparing two objects of class "densityMap"

A phytools user contacted me recently about comparing character reconstructions obtained using stochastic character mapping by the use of offset trees as in here. Although this, of course, can be done - I also have another suggestion. That is, to compare the correlation of two different densityMap style trees. A densityMap tree has the posterior density from stochastic mapping of a binary trait mapped finely on the edges and nodes of a phylogeny. For more information about this plotting method, just check out my blog.

The way we go about computing the correlation is by traversing the edges of the tree and computing the correlation as 1 if the probabilities of our character state are both 0 or both 1 on the two trees; -1 if the probability is 0 in one tree and 1 in the other; and somewhere in between for intermediate probabilities, depending on how high or low, and how similar, they are.

Here's a demo using simulated data:

## load packages
library(phytools)
## simulate a tree & two character data vectors
tree<-pbtree(n=26,tip.label=LETTERS,scale=1)
Q<-matrix(c(-1,1,1,-1),2,2)
colnames(Q)<-rownames(Q)<-letters[1:2]
x<-sim.history(tree,Q)$states

## Done simulation(s).

x

##   A   B   C   D   E   F   G   H   I   J   K   L   M   N   O   P   Q   R 
## "b" "b" "b" "b" "b" "b" "b" "b" "a" "b" "b" "b" "a" "b" "a" "b" "b" "b" 
##   S   T   U   V   W   X   Y   Z 
## "b" "b" "b" "b" "a" "a" "b" "a"

y<-sim.history(tree,Q)$states

## Done simulation(s).

y

##   A   B   C   D   E   F   G   H   I   J   K   L   M   N   O   P   Q   R 
## "b" "b" "b" "b" "b" "a" "b" "b" "b" "b" "a" "a" "a" "a" "b" "b" "b" "a" 
##   S   T   U   V   W   X   Y   Z 
## "a" "b" "b" "b" "b" "b" "a" "a"

## perform stochastic mapping
mx<-make.simmap(tree,x,nsim=100)

## make.simmap is sampling character histories conditioned on the transition matrix
## Q =

##           a         b
## a -1.258205  1.258205
## b  1.258205 -1.258205

## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =

##   a   b 
## 0.5 0.5

## Done.

my<-make.simmap(tree,y,nsim=100)

## make.simmap is sampling character histories conditioned on the transition matrix
## Q =

##           a         b
## a -1.157145  1.157145
## b  1.157145 -1.157145

## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =

##   a   b 
## 0.5 0.5

## Done.

## compute density maps
dmapx<-densityMap(mx,plot=FALSE)

## sorry - this might take a while; please be patient

dmapy<-densityMap(my,plot=FALSE)

## sorry - this might take a while; please be patient

## compare them visually
par(mfcol=c(1,2))
plot(dmapx)
plot(dmapy,direction="leftwards")

## now compute the correlation between maps
obj<-dmapx
nl<-length(obj$cols)-1
## here I recenter the probability on zero & 
## compute the vector correlation
for(i in 1:length(obj$tree$maps)){
    nx<-2*(as.numeric(names(dmapx$tree$maps[[i]]))-nl/2)/nl
    ny<-2*(as.numeric(names(dmapy$tree$maps[[i]]))-nl/2)/nl
    names(obj$tree$maps[[i]])<-round((nx*ny+1)/2*1000)
}
## change to an object of class "contMap"
obj$lims<-c(-1,1)
class(obj)<-"contMap"
obj<-setMap(obj,c("white","black"))
par(mfcol=c(1,1))
plot(obj)

That's all there is to it. We see that regions of high certainty show up with high or low correlation - depending on the similarity between maps. Regions of low certainty will always show up with near zero correlation, which is precisely what we want.

That's it.

PGLS with measurement or sampling error in the dependent variable, y

In the following demo I show one way to do phylogenetic generalized least squares (PGLS) assuming a Brownian motion model for the variance-covariance structure of the residual error, for conditions in which there is measurement error in the estimation of y, the dependent variable in our model, but not our xs. This was motivated by a question from a UMass-Boston graduate student based on the phylogenetic (generalized) ANOVA - so this is a case in which there might be easily quantifiable error in y but not in our factor or factors.

The way I have done this is not by writing a new corStruct - although one would perhaps ideally like to this.

First, this is the likelihood functin that we are going to optimize. It also returns the fitted model & likelihood when opt=FALSE.

lk<-function(sig2,y,X,C,v=NULL,opt=TRUE){
    n<-nrow(C)
    if(is.null(v)) v<-rep(0,n)
    V<-sig2*C+diag(v)
    beta<-solve(t(X)%*%solve(V)%*%X)%*%(t(X)%*%solve(V)%*%y)
    logL<--(1/2)*t(y-X%*%beta)%*%solve(V)%*%(y-X%*%beta)-
        (1/2)*determinant(V,logarithm=TRUE)$modulus[1]-
        (n/2)*log(2*pi)
    if(opt) -logL[1,1] else list(beta=beta[,1],sig2e=sig2,logL=logL[1,1])
}

Now, next, I start with an illustrative demo using multivariable regressin to show that this function will return the same fitted model (when no sampling error within species is assumed) as gls(...,correlation=corBrownian(...)). Here I go:

## load libraries
library(phytools)
library(nlme)
## simulate tree & data
tree<-pbtree(n=100,scale=1)
X<-fastBM(tree,nsim=2)
colnames(X)<-c("x1","x2")
y<-cbind(rep(1,Ntip(tree)),X)%*%c(1,2,3)+fastBM(tree)
## first, fit the model using gls
fit.gls<-gls(y~x1+x2,data=data.frame(y,X),correlation=corBrownian(1,tree),
    method="ML")
fit.gls

## Generalized least squares fit by maximum likelihood
##   Model: y ~ x1 + x2 
##   Data: data.frame(y, X) 
##   Log-likelihood: -47.8675
## 
## Coefficients:
## (Intercept)          x1          x2 
##   0.9532797   1.9886343   3.0334548 
## 
## Correlation Structure: corBrownian
##  Formula: ~1 
##  Parameter estimate(s):
## numeric(0)
## Degrees of freedom: 100 total; 97 residual
## Residual standard error: 0.9198236

## now fit it using our custom function
## the interval for optimize is specified arbitrarily
fit.lk<-optimize(lk,c(0,1000),y=y,X=cbind(rep(1,Ntip(tree)),X),C=vcv(tree))
fitted<-lk(fit.lk$minimum,y=y,X=cbind(rep(1,Ntip(tree)),X),C=vcv(tree),
    opt=FALSE)
fitted

## $beta
##                  x1        x2 
## 0.9532797 1.9886343 3.0334548 
## 
## $sig2e
## [1] 0.8460792
## 
## $logL
## [1] -47.8675

So, here we can see that it works when we do not have sampling error in y. Next, let's simulate known sampling error and try again:

## these will be our within-species sampling variances
v<-setNames(rexp(n=Ntip(tree)),tree$tip.label)
ye<-setNames(sampleFrom(xbar=y,xvar=v,n=rep(1,length(y))),rownames(y))
fit.lk<-optimize(lk,c(0,1000),y=ye,X=cbind(rep(1,Ntip(tree)),X),
    C=vcv(tree),v=v)
fitted<-lk(fit.lk$minimum,y=ye,X=cbind(rep(1,Ntip(tree)),X),C=vcv(tree),
    v=v,opt=FALSE)
fitted

## $beta
##                x1       x2 
## 0.863873 2.039439 2.969356 
## 
## $sig2e
## [1] 0.6355355
## 
## $logL
## [1] -138.4668

## compare to:
gls(ye~x1+x2,data=data.frame(ye,X),correlation=corBrownian(1,tree),method="ML")

## Generalized least squares fit by maximum likelihood
##   Model: ye ~ x1 + x2 
##   Data: data.frame(ye, X) 
##   Log-likelihood: -231.8525
## 
## Coefficients:
## (Intercept)          x1          x2 
##    1.508517    1.947806    4.554398 
## 
## Correlation Structure: corBrownian
##  Formula: ~1 
##  Parameter estimate(s):
## numeric(0)
## Degrees of freedom: 100 total; 97 residual
## Residual standard error: 5.790837

Now, it also turns out that after we have optimized sig2e we can actually coerce gls into giving us the correct fitted model & likelihood. Here, I do this by distorting the edge lengths of our tree to take into account the fitted sig2e and within-species errors in y.

tt<-tree
tt$edge.length<-tt$edge.length*fitted$sig2e
for(i in 1:length(v)){
    tip<-which(tt$tip.label==names(v)[i])
    ii<-which(tt$edge[,2]==tip)
    tt$edge.length[ii]<-tt$edge.length[ii]+v[i]
}
vv<-diag(vcv(tt))
w<-varFixed(~vv)
fit.gls<-gls(ye~x1+x2,data=data.frame(ye,X),correlation=corBrownian(1,tt),method="ML",weights=w)
fit.gls

## Generalized least squares fit by maximum likelihood
##   Model: ye ~ x1 + x2 
##   Data: data.frame(ye, X) 
##   Log-likelihood: -138.4384
## 
## Coefficients:
## (Intercept)          x1          x2 
##    0.863873    2.039439    2.969356 
## 
## Correlation Structure: corBrownian
##  Formula: ~1 
##  Parameter estimate(s):
## numeric(0)
## Variance function:
##  Structure: fixed weights
##  Formula: ~vv 
## Degrees of freedom: 100 total; 97 residual
## Residual standard error: 1.0169

## compare to:
fitted

## $beta
##                x1       x2 
## 0.863873 2.039439 2.969356 
## 
## $sig2e
## [1] 0.6355355
## 
## $logL
## [1] -138.4668

Now, here's a similar demo using the phylogenetic (generalized) ANOVA:

## evolve a factor on the tree
Q<-matrix(c(-2,1,1,1,-2,1,1,1,-2),3,3)
colnames(Q)<-rownames(Q)<-letters[1:3]
x<-as.factor(sim.history(tree,Q)$states)

## Done simulation(s).

y<-fastBM(tree)
y[x=="a"]<-y[x=="a"]+1
y[x=="b"]<-y[x=="b"]+2
y[x=="c"]<-y[x=="c"]+3
fit.gls<-gls(y~x,data=data.frame(y,x),correlation=corBrownian(1,tree),method="ML")
fit.gls

## Generalized least squares fit by maximum likelihood
##   Model: y ~ x 
##   Data: data.frame(y, x) 
##   Log-likelihood: -58.60169
## 
## Coefficients:
## (Intercept)          xb          xc 
##   0.8929341   0.8188519   1.8667981 
## 
## Correlation Structure: corBrownian
##  Formula: ~1 
##  Parameter estimate(s):
## numeric(0)
## Degrees of freedom: 100 total; 97 residual
## Residual standard error: 1.024053

summary(fit.gls)

## Generalized least squares fit by maximum likelihood
##   Model: y ~ x 
##   Data: data.frame(y, x) 
##        AIC      BIC    logLik
##   125.2034 135.6241 -58.60169
## 
## Correlation Structure: corBrownian
##  Formula: ~1 
##  Parameter estimate(s):
## numeric(0)
## 
## Coefficients:
##                 Value Std.Error   t-value p-value
## (Intercept) 0.8929341 0.3801839  2.348690  0.0209
## xb          0.8188519 0.1411306  5.802084  0.0000
## xc          1.8667981 0.1658942 11.252944  0.0000
## 
##  Correlation: 
##    (Intr) xb    
## xb -0.181       
## xc -0.169  0.431
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -3.51162932 -1.05473759 -0.05870776  0.70967087  2.99993479 
## 
## Residual standard error: 1.024053 
## Degrees of freedom: 100 total; 97 residual

X<-model.matrix(~x)
fit.lk<-optimize(lk,c(0,1000),y=y,X=X,C=vcv(tree))
fitted<-lk(fit.lk$minimum,y=y,X=X,C=vcv(tree),opt=FALSE)
fitted

## $beta
## (Intercept)          xb          xc 
##   0.8929341   0.8188519   1.8667981 
## 
## $sig2e
## [1] 1.048674
## 
## $logL
## [1] -58.60169

## now with sampling error in y
ye<-setNames(sampleFrom(xbar=y,xvar=v,n=rep(1,length(y))),rownames(y))
## fit model
fit.lk<-optimize(lk,c(0,1000),y=ye,X=X,C=vcv(tree),v=v)
fitted<-lk(fit.lk$minimum,y=ye,X=X,C=vcv(tree),v=v,opt=FALSE)
fitted

## $beta
## (Intercept)          xb          xc 
##   0.7462831   1.2109657   1.9742749 
## 
## $sig2e
## [1] 1.338498
## 
## $logL
## [1] -136.757

## coerce tree again:
tt<-tree
tt$edge.length<-tt$edge.length*fitted$sig2e
for(i in 1:length(v)){
    tip<-which(tt$tip.label==names(v)[i])
    ii<-which(tt$edge[,2]==tip)
    tt$edge.length[ii]<-tt$edge.length[ii]+v[i]
}
vv<-diag(vcv(tt))
w<-varFixed(~vv)
fit.gls<-gls(ye~x,data=data.frame(ye,x),correlation=corBrownian(1,tt),method="ML",weights=w)
fit.gls

## Generalized least squares fit by maximum likelihood
##   Model: ye ~ x 
##   Data: data.frame(ye, x) 
##   Log-likelihood: -134.6691
## 
## Coefficients:
## (Intercept)          xb          xc 
##   0.7462831   1.2109657   1.9742749 
## 
## Correlation Structure: corBrownian
##  Formula: ~1 
##  Parameter estimate(s):
## numeric(0)
## Variance function:
##  Structure: fixed weights
##  Formula: ~vv 
## Degrees of freedom: 100 total; 97 residual
## Residual standard error: 0.8591587

summary(fit.gls)

## Generalized least squares fit by maximum likelihood
##   Model: ye ~ x 
##   Data: data.frame(ye, x) 
##        AIC     BIC    logLik
##   277.3383 287.759 -134.6691
## 
## Correlation Structure: corBrownian
##  Formula: ~1 
##  Parameter estimate(s):
## numeric(0)
## Variance function:
##  Structure: fixed weights
##  Formula: ~vv 
## 
## Coefficients:
##                 Value Std.Error  t-value p-value
## (Intercept) 0.7462831 0.3900422 1.913339  0.0587
## xb          1.2109657 0.2408905 5.027037  0.0000
## xc          1.9742749 0.3019824 6.537715  0.0000
## 
##  Correlation: 
##    (Intr) xb    
## xb -0.295       
## xc -0.260  0.348
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -3.2076462 -0.8804235 -0.1460402  0.8224632  3.4883858 
## 
## Residual standard error: 0.8591587 
## Degrees of freedom: 100 total; 97 residual

OK, that's it.

How to get the colors of the tips in a plotted "contMap" object

Today I received an interesting question about the phytools function contMap. contMap, remember, is a function that maps the observed & reconstructed value of a continuous character on to the tree. The question is as follows:

“How can I recover the color codes that are used to paint the tips (or the last bit) of a phylogeny after a reconstruction with contMap? I have a character mapped on my phylogeny and I want to use the same tip colors in a scatter plot (first two components of a PCA in this case)… each species in the plot should have the same color as in the tree. I tried going through the code but I get lost at some point and can't figure out how the last bit of the terminal branches are indexed and how the color assigment is done…”

An object created by contMap is actually just the same object that is read & plotted by plotSimmap, which plot.contMap uses internally. The colors are encoded in the "contMap" element cols, and the corresponding indices for the colors mapped along each edge, in the object element tree$maps.

To demo how we can pull the colors of all mapped species from this object, see the following demo:

library(phytools)
## simulate tree & data
tree<-pbtree(n=26,tip.label=LETTERS)
X<-fastBM(tree,nsim=2)
colnames(X)<-c("trait 1","trait 2")
## contMap for the first character only
obj<-contMap(tree,X[,1])

## function to pull out the color of a given tip
foo<-function(obj,tip){
    jj<-which(obj$tree$tip.label==tip)
    kk<-which(obj$tree$edge[,2]==jj)
    setNames(obj$cols[names(obj$tree$maps[[kk]])[length(obj$tree$maps[[kk]])]],
        NULL)
}
colors<-sapply(obj$tree$tip.label,foo,obj=obj)
colors

##           A           B           C           D           E           F 
## "#004BFFFF" "#0023FFFF" "#00FF5FFF" "#0018FFFF" "#0EFF00FF" "#10FF00FF" 
##           G           H           I           J           K           L 
## "#3000FFFF" "#0062FFFF" "#0F00FFFF" "#00A1FFFF" "#00FF64FF" "#00FF4DFF" 
##           M           N           O           P           Q           R 
## "#00FF33FF" "#00FF49FF" "#00FF79FF" "#FF7F00FF" "#FF2400FF" "#FF0300FF" 
##           S           T           U           V           W           X 
## "#00FF03FF" "#3EFF00FF" "#D9FF00FF" "#FFAE00FF" "#FF4A00FF" "#FFF500FF" 
##           Y           Z 
## "#FFA400FF" "#FF7800FF"

OK, now how about using these to plot in a scatterplot? We can do that in our two dimensional simulated space:

plot(X,col=colors,pch=19,cex=1.5)
## just to position the labels a bit more nicely
ii<-X[,1]>(max(X[,1])+min(X[,1]))/2
text(X[ii,],rownames(X)[ii],pos=2)
text(X[!ii,],rownames(X)[!ii],pos=4)

That's all there is to it!