The MCCR test for Pybus & Harvey's (2000)
γ statistic from incompletely sampled phylogenies used to be implemented in the
no-longer-available-on-CRAN
laser package of
Dan Rabosky. The method
is pretty simple, though, & since we use it in teaching I figured it would be easy enough to
add to *phytools*.

Here is a function to do the test, as well as `print`

& `plot`

methods
for the resultant object:

```
mccr<-function(obj,rho=1,nsim=100,plot=TRUE){
N<-round(Ntip(obj$tree)/rho)
tt<-pbtree(n=N,nsim=nsim)
foo<-function(x,m) drop.tip(x,sample(x$tip.label,m))
tt<-lapply(tt,foo,m=N-Ntip(obj$tree))
g<-sapply(tt,function(x) ltt(x,plot=FALSE)$gamma)
P<-if(obj$gamma>median(g)) 2*mean(g>=obj$gamma) else 2*mean(g<=obj$gamma)
result<-list(gamma=obj$gamma,"P(two-tailed)"=P,null.gamma=g)
class(result)<-"mccr"
result
}
print.mccr<-function(x,digits=4,...){
cat("Object of class \"mccr\" consisting of:\n\n")
cat(paste("(1) A value for Pybus & Harvey's \"gamma\"",
" statistic of ",round(x$gamma,digits),".\n\n",sep=""))
cat(paste("(2) A two-tailed p-value from the MCCR test of ",
round(x$'P(two-tailed)',digits),".\n\n", sep = ""))
cat(paste("(3) A simulated null-distribution of gamma from ",
length(x$null.gamma)," simulations.\n\n",sep=""))
}
plot.mccr<-function(x,...){
hist(x$null.gamma,breaks=min(c(max(12,round(length(x$null.gamma)/10)),20)),
xlim=range(c(x$gamma,x$null.gamma)),
main=expression(paste("null distribution of ",
gamma)),xlab=expression(gamma),col="lightgrey")
arrows(x0=x$gamma,y0=par()$usr[4],y1=0,length=0.12,
col=make.transparent("blue",0.5),lwd=2)
text(x$gamma,0.98*par()$usr[4],
expression(paste("observed value of ",gamma)),
pos=if(x$gamma>mean(x$null.gamma)) 2 else 4)
}
```

Now let's see how to use it. First, I'll run it on the completely sampled tree:

```
## pure birth tree with complete sampling:
pb<-pbtree(n=100)
lineages<-ltt(pb)
```

```
lineages
```

```
## Object of class "ltt" containing:
##
## (1) A phylogenetic tree with 100 tips and 99 internal nodes.
##
## (2) Vectors containing the number of lineages (ltt) and branching times (times) on the tree.
##
## (3) A value for Pybus & Harvey's "gamma" statistic of -0.2688, p-value = 0.7881.
```

```
## MCCR test assuming complete sampling:
object<-mccr(lineages,rho=1,nsim=200)
object
```

```
## Object of class "mccr" consisting of:
##
## (1) A value for Pybus & Harvey's "gamma" statistic of -0.2688.
##
## (2) A two-tailed p-value from the MCCR test of 0.67.
##
## (3) A simulated null-distribution of gamma from 200 simulations.
```

```
plot(object)
```

A neat thing to show is that the MCCR test (based on simulation) will give a P-value that converges on the P-value from the parametric test if you run enough simulations. Let's see that:

```
object<-mccr(lineages,rho=1,nsim=10000)
lineages
```

```
## Object of class "ltt" containing:
##
## (1) A phylogenetic tree with 100 tips and 99 internal nodes.
##
## (2) Vectors containing the number of lineages (ltt) and branching times (times) on the tree.
##
## (3) A value for Pybus & Harvey's "gamma" statistic of -0.2688, p-value = 0.7881.
```

```
object
```

```
## Object of class "mccr" consisting of:
##
## (1) A value for Pybus & Harvey's "gamma" statistic of -0.2688.
##
## (2) A two-tailed p-value from the MCCR test of 0.7942.
##
## (3) A simulated null-distribution of gamma from 10000 simulations.
```

Next, we can simulate incomplete sampling, see the effect on the P-value from `ltt`

,
and then compare that to a P-value obtained from the MCCR test:

```
incomplete.pb<-drop.tip(pb,sample(pb$tip.label,round(0.5*Ntip(pb))))
lineages<-ltt(incomplete.pb,plot=FALSE)
lineages
```

```
## Object of class "ltt" containing:
##
## (1) A phylogenetic tree with 50 tips and 49 internal nodes.
##
## (2) Vectors containing the number of lineages (ltt) and branching times (times) on the tree.
##
## (3) A value for Pybus & Harvey's "gamma" statistic of -0.9876, p-value = 0.3234.
```

```
object<-mccr(lineages,rho=0.5,nsim=200)
object
```

```
## Object of class "mccr" consisting of:
##
## (1) A value for Pybus & Harvey's "gamma" statistic of -0.9876.
##
## (2) A two-tailed p-value from the MCCR test of 0.76.
##
## (3) A simulated null-distribution of gamma from 200 simulations.
```

Finally, for fun we can simulate a tree with a shape that should correspond to a positive
value of γ. Here, I'll do it using the *phytools* internal function
`ebTree`

. We can then subsample the tree to simulate incomplete sampling.

```
eb<-phytools:::ebTree(pbtree(n=200,scale=1),-1.2)
ltt(eb)
```

```
## Object of class "ltt" containing:
##
## (1) A phylogenetic tree with 200 tips and 199 internal nodes.
##
## (2) Vectors containing the number of lineages (ltt) and branching times (times) on the tree.
##
## (3) A value for Pybus & Harvey's "gamma" statistic of 2.0067, p-value = 0.0448.
```

```
incomplete.eb<-drop.tip(eb,sample(eb$tip.label,round(0.5*Ntip(eb))))
ltt(incomplete.eb)
```

```
## Object of class "ltt" containing:
##
## (1) A phylogenetic tree with 100 tips and 99 internal nodes.
##
## (2) Vectors containing the number of lineages (ltt) and branching times (times) on the tree.
##
## (3) A value for Pybus & Harvey's "gamma" statistic of -0.2614, p-value = 0.7938.
```

```
object<-mccr(ltt(incomplete.eb,plot=FALSE),rho=0.5,nsim=200)
object
```

```
## Object of class "mccr" consisting of:
##
## (1) A value for Pybus & Harvey's "gamma" statistic of -0.2614.
##
## (2) A two-tailed p-value from the MCCR test of 0.06.
##
## (3) A simulated null-distribution of gamma from 200 simulations.
```

```
plot(object)
```

Neat.

This is not yet in *phytools*, but I will add documentation & put it in a future release.

Cheap Watches At Walmart

ReplyDeleteBest Watches Mvmt

Hublot Watches Images

Audemars Piguet Watches Toronto

Panerai Watches Fake

Rolex Watches 41mm

Copy Watches Rolex

Richard Mille Watches Price In Pakistan

Swiss Watches Second Hand

replica Watches Turkey Online

My name is JOHN ROBET am from SA I want to share a testimony of how Dr.Olu herbal mixture cream saves me from shame and disgrace, my penis was a big problem to me as the size was really so embarrassing,and i was also having weak erection problem. I can't make love to my wife and my penis was just too small a full grown man like me having 4 inches penis and to worsen it i don't last in sex i cant even last two minutes it was really a thing of shame to me. My wife was really tired of me because my sex life was very poor,she never enjoyed sex,i was always thinking and searching for solutions everywhere until when i saw a testimony of how Dr.Olu herbal mixture cream have been helping people regarding their sex life, so i decided to give him a try and to my greatest surprise in less than one week of taking the herbs my penis grow to 8 inches i couldn't believe my eyes and as i speak now my penis is now 8 inches and i do not have week erection again. I can make love to my wife longer in bed. And my marriage is now stable,my wife now enjoy me very well in bed. can contact him drolusoltionhome@gmail.com {) or call or what-apps him through +2348140654426 and also his you-tube page https://youtu.be/VMTo3gqbO08 . .he also specialize on the following things

ReplyDeleteBREAST ENLARGEMENT

BREAST LIFT & REDUCTION

PENILE/PENIS ENLARGEMENT

GENERAL BODY SLIMMING

HIPS & THIGH REDUCTION

REGAIN VIRGINITY BACK

Thanks for the Enlarging my penis sir, you indeed save my marriage...I am really grateful sir,