Permalink
Browse files

print method for phylANOVA

  • Loading branch information...
liamrevell committed Nov 11, 2017
1 parent e18453a commit 2c21e9ab67c3de0efca1beab65b5d24bde08501b
Showing with 38 additions and 4 deletions.
  1. +1 −0 NAMESPACE
  2. +37 −4 R/phylANOVA.R
View
@@ -124,6 +124,7 @@ S3method(logLik, rerootingMethod)
S3method(print, fitBayes)
S3method(plot, fitBayes)
S3method(print, ratebystate)
S3method(print, phylANOVA)
importFrom(animation, ani.options, ani.record, ani.replay, saveVideo)
importFrom(ape, .PlotPhyloEnv, .uncompressTipLabel, ace, all.equal.phylo, as.DNAbin, as.phylo, bind.tree, branching.times, collapse.singles)
View
@@ -1,6 +1,6 @@
# function conducts phylogenetic ANOVA & posthoc tests
# some code from phy.anova() in "geiger"
# written by Liam Revell 2011, 2015
# written by Liam Revell 2011, 2015, 2017
phylANOVA<-function(tree,x,y,nsim=1000,posthoc=TRUE,p.adj="holm"){
if(is.null(names(x))){
@@ -17,7 +17,8 @@ phylANOVA<-function(tree,x,y,nsim=1000,posthoc=TRUE,p.adj="holm"){
sig2<-mean(pic(y,multi2di(tree))^2) # compute BM rate for y
x<-as.factor(x) # change x to factor
m<-length(levels(x))
F.obs<-anova(lm(y~x))[1,4] # F on empirical data
aov<-anova(lm(y~x))
F.obs<-aov[1,4] # F on empirical data
if(posthoc) T.obs<-tTests(x,y) # empirical ts
sims<-fastBM(tree,sig2=sig2,nsim=(nsim-1)) # simulate
F.null<-vector()
@@ -39,8 +40,12 @@ phylANOVA<-function(tree,x,y,nsim=1000,posthoc=TRUE,p.adj="holm"){
# control for multiple tests (if p.adj!="none")
P.T[lower.tri(P.T)]<-p.adjust(P.T[lower.tri(P.T)],method=p.adj)
for(i in 1:m) for(j in i:m) P.T[i,j]<-P.T[j,i]
return(list(F=F.obs,Pf=P.F,T=T.obs,method=p.adj,Pt=P.T))
} else return(list(F=F.obs,Pf=P.F))
obj<-list(F=F.obs,Pf=P.F,T=T.obs,method=p.adj,Pt=P.T,
"Sum Sq"=aov$"Sum Sq","Mean Sq"=aov$"Mean Sq")
} else obj<-list(F=F.obs,Pf=P.F,"Sum Sq"=aov$"Sum Sq",
"Mean Sq"=aov$"Mean Sq")
class(obj)<-"phylANOVA"
obj
}
# computes pairwise t-statistics with a pooled SD
@@ -66,3 +71,31 @@ tTests<-function(x,y){
for(i in 1:m) for(j in 1:m) T[i,j]<-compare.levels(levels(x)[i],levels(x)[j])
return(T)
}
## S3 print method
print.phylANOVA<-function(x,digits=6,...){
cat("ANOVA table: Phylogenetic ANOVA\n\n")
cat("Response: y\n")
object<-data.frame(round(x$"Sum Sq",digits),
round(x$"Mean Sq",digits),
c(round(x$F,digits),""),
c(round(x$Pf,digits),""))
colnames(object)<-c("Sum Sq","Mean Sq","F value",
"Pr(>F)")
rownames(object)<-c("x","Residual")
print(object)
cat("\nP-value based on simulation.\n")
cat("---------\n")
cat("\n")
if(!is.null(obj$T)){
cat(paste("Pairwise posthoc test using method = \"",
x$method,"\"\n\n",sep=""))
cat("Pairwise t-values:\n")
print(round(x$T,digits))
cat("\nPairwise corrected P-values:\n")
print(round(x$Pt,digits))
cat("---------\n")
cat("\n")
}
}

0 comments on commit 2c21e9a

Please sign in to comment.