Type I error & power for ratebytree with data type="discrete"

Yesterday I posted about an extension of the phytools function ratebytree to fit different rates of discrete character evolution to different trees, and then to compare this model to one in which the rates of evolutionary change between states have been forced to be the same.

I thought it might be neat to do some very simple analyses of the type I error rate of the method - as well as of its power to detect a difference in rate when one has been simulated. This is something that we have done (more thoroughtly) in our manuscript describing the statistical properties of ratebytree for continuous traits.

To accomplish this, I will first need to simulate under the null & measure the fraction of times amongst my simulations for which the null is rejected (type I errors). Next, I will simulate under various differences in evolutionary rate between trees and measure the fraction of times that the null hypothesis is rejected when false (power).

Here I go. First, with type I error:

library(phytools)
packageVersion("phytools") ## check we have 0.6.25 or above

## [1] '0.6.25'

## simulate pairs of trees
trees<-replicate(200,pbtree(n=100,nsim=2),simplify=FALSE)
head(trees,n=3)

## [[1]]
## 2 phylogenetic trees
## 
## [[2]]
## 2 phylogenetic trees
## 
## [[3]]
## 2 phylogenetic trees

## now simulate pairs of datasets on each of these trees
Q<-matrix(c(-1,1,1,-1),2,2,dimnames=list(letters[1:2],
    letters[1:2]))
foo<-function(t,Q)
    list(as.factor(sim.history(t[[1]],Q,message=F)$states),
        as.factor(sim.history(t[[2]],Q,message=F)$states))
x<-lapply(trees,foo,Q=Q)
head(x,n=3)

## [[1]]
## [[1]][[1]]
##   t6  t52  t99 t100  t87  t88  t67  t68  t20  t95  t96  t81  t82  t21  t22 
##    a    a    a    a    b    b    b    a    b    a    a    a    a    b    b 
##  t15   t1   t2  t47  t48  t75  t76  t69  t70  t91  t92  t73  t74  t46  t12 
##    b    b    b    b    b    b    b    a    a    a    a    b    a    a    a 
##  t16  t17  t44  t97  t98  t63  t32  t24  t27  t77  t78   t3  t18  t19  t64 
##    b    a    a    b    b    a    a    b    a    a    a    b    a    a    b 
##  t83  t84   t8  t10  t56  t57  t42  t43  t23  t38  t39  t40  t41  t33  t71 
##    b    a    b    a    b    a    b    a    a    b    b    a    b    a    a 
##  t72  t28  t25  t13  t30  t31   t5  t61  t62  t29  t89  t90  t85  t49  t36 
##    a    a    a    a    a    a    b    b    a    a    b    a    b    a    b 
##  t37  t26  t79  t80  t53  t54  t93  t94  t86  t45  t50  t51  t14  t58  t59 
##    b    b    a    a    b    a    a    a    a    b    b    b    a    a    a 
##   t7   t4   t9  t34  t35  t11  t55  t60  t65  t66 
##    a    a    a    a    a    b    a    a    b    a 
## Levels: a b
## 
## [[1]][[2]]
##  t26  t51  t52  t16  t91  t92  t28  t24  t43  t44  t67  t99 t100  t45  t46 
##    a    a    a    a    b    b    b    b    b    b    b    b    b    b    b 
##  t76  t83  t84  t74  t75  t20  t55  t56   t5   t6  t63  t64  t48  t49  t72 
##    a    a    a    a    a    a    a    a    b    b    a    a    a    a    a 
##  t73  t25  t36  t59  t60  t40  t41  t68  t69  t54  t47   t4  t13  t17  t18 
##    a    b    a    a    b    a    a    a    a    a    a    a    b    a    a 
##  t79  t87  t88  t37  t35   t8  t29  t30  t89  t90  t42  t12   t9  t38  t39 
##    b    a    a    a    a    b    b    b    a    a    b    a    b    b    b 
##  t14  t95  t96  t11  t31  t32  t85  t86  t19  t97  t98  t82  t53  t61  t62 
##    a    b    b    b    b    a    b    b    b    a    a    a    b    b    a 
##  t10   t2  t57  t58   t7  t70  t71   t3  t33  t34  t50  t77  t78  t65  t66 
##    a    b    a    a    b    b    b    b    b    b    a    b    b    b    a 
##  t27  t80  t81   t1  t21  t93  t94  t22  t23  t15 
##    b    b    b    a    b    b    b    b    a    a 
## Levels: a b
## 
## 
## [[2]]
## [[2]][[1]]
##  t40  t97  t98  t66  t67  t53  t43  t80  t81  t27  t54  t55  t10  t71  t72 
##    b    a    a    a    a    a    a    a    a    b    b    b    a    a    a 
##  t62  t63  t41  t42  t50  t75  t85  t86  t12   t7  t60  t61  t58  t59  t19 
##    a    b    b    b    b    b    b    b    a    a    a    a    a    a    a 
##  t26  t46  t47  t35  t37  t91  t92  t24   t1  t20  t21  t38  t39  t87  t88 
##    a    a    a    a    a    a    a    b    a    a    a    b    b    a    a 
##  t28   t5   t3   t4  t73  t74   t9   t2  t36  t48  t49  t82  t83  t84  t18 
##    a    a    b    a    a    a    a    a    b    a    a    a    a    a    a 
##  t29  t30  t17  t51  t52  t64  t65  t22  t23  t11  t70  t76  t77  t78  t79 
##    a    b    b    b    b    a    a    a    b    b    a    b    a    a    a 
##  t93  t94  t14  t56  t57  t68  t69  t25  t89  t90  t34  t13  t15  t16  t31 
##    a    b    a    a    a    a    a    a    a    b    a    a    b    b    a 
##  t32  t33  t99 t100  t95  t96   t6  t44  t45   t8 
##    a    a    a    a    a    a    a    b    b    a 
## Levels: a b
## 
## [[2]][[2]]
##  t53  t54  t61  t62  t11   t8   t3  t51  t52  t88  t89  t50  t99 t100  t27 
##    a    a    a    a    b    b    b    b    b    a    a    b    b    b    b 
##  t17  t18  t80  t81  t69  t70  t55  t12   t6  t56  t57  t75  t76  t68  t82 
##    a    b    b    b    b    b    b    a    a    a    a    a    a    b    a 
##  t83  t49  t93  t94  t90  t16  t23  t59  t60  t34  t28  t29  t32  t33  t24 
##    a    b    a    a    a    b    b    a    a    a    b    b    b    b    a 
##  t47  t48  t21   t4   t2   t9  t25  t26  t58  t95  t96  t42  t43   t7  t66 
##    b    b    a    a    a    b    b    b    a    a    a    b    a    a    b 
##  t67  t77  t97  t98  t35  t36   t1  t37  t38  t22  t41  t71  t72  t46  t86 
##    a    b    a    a    b    a    a    b    a    b    b    b    a    b    b 
##  t87  t73  t74  t44  t45  t39  t40  t13  t91  t92  t14  t15  t64  t65  t63 
##    b    b    b    b    b    b    b    a    a    a    a    b    b    b    b 
##  t84  t85   t5  t10  t30  t31  t78  t79  t19  t20 
##    b    b    a    a    b    b    a    a    a    a 
## Levels: a b
## 
## 
## [[3]]
## [[3]][[1]]
##  t56  t57  t71  t72  t97  t98  t77  t78   t4   t2  t59  t60  t39  t20  t26 
##    a    a    b    a    b    b    b    b    b    a    a    a    a    b    b 
##  t27  t10  t13  t14   t9  t30  t89  t90  t67  t68  t51  t15  t36  t37  t11 
##    a    b    a    b    b    a    a    a    b    b    b    a    b    b    a 
##  t49  t75  t76  t55  t31  t44  t45   t8  t53  t54  t48  t91  t92  t40  t41 
##    b    a    a    b    b    a    b    b    b    b    b    b    b    a    b 
##  t17  t46  t47  t95  t96  t29  t25  t18  t19   t6   t1  t23  t24  t73  t74 
##    b    b    a    a    a    a    a    b    a    a    a    a    b    b    b 
##   t7  t22  t58  t79  t80  t81  t82  t16  t52  t69  t70  t61  t62  t85  t86 
##    a    b    b    b    b    b    b    b    b    a    a    b    b    a    a 
##  t83  t84  t33  t34  t99 t100  t50  t32  t28  t93  t94  t35  t38  t63  t64 
##    b    b    a    a    b    b    b    a    b    a    a    a    a    a    a 
##  t87  t88  t42  t43  t21  t12   t5   t3  t65  t66 
##    a    a    b    b    b    b    b    a    a    a 
## Levels: a b
## 
## [[3]][[2]]
##   t8   t9  t58  t63  t64  t73  t92  t93  t69  t56  t57  t25  t60  t85  t86 
##    b    b    a    a    b    b    b    a    b    b    b    b    b    b    b 
##   t7  t65  t66  t37  t38   t2  t12  t49  t54  t55  t18  t15  t16  t36  t83 
##    a    a    b    b    b    a    b    a    a    a    a    b    a    a    a 
##  t84  t81  t82  t27  t89  t90   t6  t74  t75  t59  t76  t77  t47  t48  t98 
##    a    b    a    a    a    a    b    a    a    a    b    b    b    b    b 
##  t99 t100  t70  t94  t95  t91  t79  t80  t29  t23  t78  t96  t97   t1  t21 
##    b    b    a    b    b    b    a    a    a    b    b    b    b    b    a 
##  t22  t10  t43  t44  t14  t52  t53  t19  t32  t33  t17  t28  t30  t31  t39 
##    a    b    a    b    b    b    a    b    b    a    a    a    b    a    b 
##  t40  t71  t72  t61  t62  t20  t11   t4  t50  t51   t5  t45  t46  t26  t34 
##    b    b    b    a    b    b    a    a    a    b    b    b    a    b    a 
##  t35   t3  t13  t24  t41  t42  t87  t88  t67  t68 
##    a    b    b    a    b    b    b    b    a    b 
## Levels: a b

## now fit our models
fitsER<-mapply(ratebytree,trees=trees,x=x,MoreArgs=list(model="ER"),
    SIMPLIFY=FALSE)
head(fitsER,n=3)

## [[1]]
## ML common-rate model:
##  q   k   logL
## value    0.7339  1   -125.016
## 
## Model fitting method was "optimize".
## 
## ML multi-rate model:
##  q   k   logL
## tree1    1.2666  2   -124.1744
## tree2    0.5602
## 
## Model fitting method was "nlminb".
## 
## Likelihood ratio: 1.6831 
## P-value (based on X^2): 0.1945 
## 
## 
## [[2]]
## ML common-rate model:
##  q   k   logL
## value    0.5811  1   -117.8442
## 
## Model fitting method was "optimize".
## 
## ML multi-rate model:
##  q   k   logL
## tree1    0.4412  2   -117.5178
## tree2    0.687
## 
## Model fitting method was "nlminb".
## 
## Likelihood ratio: 0.6527 
## P-value (based on X^2): 0.4191 
## 
## 
## [[3]]
## ML common-rate model:
##  q   k   logL
## value    1.0253  1   -123.4968
## 
## Model fitting method was "optimize".
## 
## ML multi-rate model:
##  q   k   logL
## tree1    0.8676  2   -123.228
## tree2    1.3832
## 
## Model fitting method was "nlminb".
## 
## Likelihood ratio: 0.5376 
## P-value (based on X^2): 0.4634

## finally, pull out the P-values for each analysis:
P<-sapply(fitsER,function(x) x$P.chisq)
obj<-hist(P,breaks=seq(0,1,by=0.05),plot=F)
plot(obj$mids,obj$counts/sum(obj$counts),type="h",
    lwd=18,col=make.transparent("blue",0.4),lend=1,
    xlab="P-value",ylab="relative frequency",ylim=c(0,0.2))
title(main="P-values for data simulated under the null")
abline(h=0.05,col=make.transparent("red",0.5),lwd=4,
    lty="dashed",lend=1)

mean(P<=0.05)

## [1] 0.06

We can see that the type I error rate is elevated above its nominal level, although just barely in this case. This should not be especially surprising as our expectation that the likelihood-ratio has a Χ² distribution under the null is asymptotic with N, meaning that we only expect it to hold true for relatively large sample sizes. Conversely, of course, we should expect the type I error to increase for tests conducted using smaller trees (which it does - though I have not shown that here).

Next, we can examine the power of the method to detect a difference in rate between trees if one is present. Here I will simulate a transition rate q₁ for tree 1 fixed at 1.0, and vary the transition rate q₂ for the second tree between 1.0 and 4.0 in intervals of 0.2. Since this is a whole bunch of simulations, I'm only going to use the first 50 sets of trees from before:

trees<-trees[1:50]
q1<-1
q2<-seq(1,4,by=0.2)
foo<-function(q2,q1,trees){
    Q1<-matrix(c(-q1,q1,q1,-q1),2,2,dimnames=list(letters[1:2],
        letters[1:2]))
    Q2<-matrix(c(-q2,q2,q2,-q2),2,2,dimnames=list(letters[1:2],
        letters[1:2]))
    x<-lapply(trees,function(t,Q1,Q2)
        list(as.factor(sim.history(t[[1]],Q1,message=F)$states),
        as.factor(sim.history(t[[2]],Q2,message=F)$states)),
        Q1=Q1,Q2=Q2)
    fits<-mapply(ratebytree,trees=trees,x=x,MoreArgs=list(model="ER"),
        SIMPLIFY=FALSE)
    sapply(fits,function(x) x$P.chisq)
}
P<-sapply(q2,foo,q1=q1,trees=trees)
power<-apply(P,2,function(x) mean(x<=0.05))
plot(power~q2,cex=1.8,bg="grey",pch=21,xlab=expression(q[2]/q[1]))

So, power increases for an increase in the difference in rate between trees 1 & 2 - just as we'd expect. Note that for this smaller set of simulations we see a more dramatically increased type I error rate under the null hypothesis of no difference in rate (q₂ / q₁ = 1.0).

More later.