R-beta: Re: Post-hoc tests

[Matthew Kay]

Jim,

  Thanks for the code! I'm new to R so now I have to figure out how to
'install' it. ;-) Thomas Lumley has also fowarded me some of his code
(bonferroni, holm, and hochberg), which I have included below for your interest
(you may have already received it via the mailing list).

  Since my inquiry about post-hoc tests I have learned that my data requires a
more power method than the Bonferroni. I'm looking for a Fisher routine now
(also known as LSD), which from what I am told is quite powerful in detecting
significant differences among groups. My objective is to determine "which
groups are significantly different from every other group?" and this involves
comparing each groups average to the average of every other group, recursively.
For example, here is a sample output from SYSTAT:

1=group A
2=group B
3=group C

    1      2      3
1  0.000
2  0.988  0.000
3  0.001  0.111  0.000

The numbers in the table are p values, indicating significance between group
means. I'm not sure if Thomas Lumley's code does this, I haven't had time to
try it out yet. From what you have told me I assume that your code doesn't
compare group averages to each other, it compares each group's average to the
average of the entire dataset. That could be useful to me too and I do
appreciate you sending it to me. How hard would it be to modify your code to
compare each groups average to the average of every other group?

I don't think I'm familiar with rank-order tests but I'm interested in taking a
look at what you have written.

Matt

> On Jul 16, 10:12am, Jim Lemon wrote:
> Subject: R-beta: Re: Post-hoc tests
>
> Matt,
>
> Here's a Bonferroni-corrected multiple one-sample t-test that I wrote
> some years ago.  It took a while to get it into R, as na.omit doesn't
> seem to handle vectors and I had to write a quick kludge (na.remove).
> Another more general point was that I discovered that the help page for
> t.test gives the name "parameters" for the degrees of freedom, as in S.
> However, the name is really "parameter" as I found after quite a bit of
> head-scratching and messing around.
>
> I have some rank-order tests which are similar, and I will see if I can
> dig them up as well.  These (and the above) were all written to test
> comparisons of the form "which groups are significantly different from
> the average of all groups?", which may not be what you are looking for,
> but I would be glad to modify them if I can.
>
> Jim
>

---------------------------------------------------------
From:    Thomas Lumley <thomas at biostat.washington.edu>
Date:    Tue Jul 14,  5:33pm -0700
To:      Matthew Kay <mwkay>

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())
}

-- 
Matthew W. Kay                       Office: (314) 935-7562
Washington University, St. Louis        Fax: (314) 935-4014
E-mail: mwkay at spud.me.wustl.edu
WWW: http://userfs.cec.wustl.edu/~mwk2/HOMEPAGE/index.html
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help 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-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._