Here is code for three p-value based post hoc tests: bonferroni, holm, and
hochberg. Holm is better than Bonferroni under all circumtances. Hochberg
is more powerful than Holm, but under some slightly pathological
situations it can exceed the nominal Type I error rate.
I haven't looked at this code for a long time, but I think it works
correctly. Refeerences for the methods are in "Adjusted p-values for
simultaneous inference" by S. Paul Wright in Biometrics some time in
the early 90s (I'm in the wrong country at the moment so I don't have
the
full reference).
-thomas
-------------------------
p.adjust.holm <-function(p,n=length(p)) {
##n<-length(p)
r<-rank(p)
index<-order(p)
qi<-p*(n+1-r)
for (i in 2:length(p)) {
qi[index[i]]<-max(qi[index[i]], qi[index[i-1]])
}
list(adjp=pmin(qi,1),p=p,method="Holm")
}
p.adjust.hochberg <-function(p) {
n<-length(p)
r<-rank(p)
index<-order(p)
qi<-p*(n+1-r)
for (i in (n-1):1) {
qi[index[i]]<-min(qi[index[i]], qi[index[i+1]])
}
list(adjp=qi,p=p,method="Hochberg")
}
p.adjust.bonferroni<-function(p,n=length(p)){
list(adjp=pmin(p*n,1),p=p,method="Bonferroni")
}
p.adjust<-function(p,method=c("hochberg","holm","bonferroni"),...){
how<-pmatch(method[1],c("hochberg","holm","bonferroni"))
if (is.na(how)) stop(paste("Don't know method:",method))
m<-match.call()
m[[1]]<-as.name(paste("p.adjust",c("hochberg","holm","bonferroni")[how],sep="."))
m$method<-NULL
eval(m,sys.parent())
}
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To:
r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._