On Jul 20, 2013, at 10:37 AM, Scott Robinson wrote:
> Dear List,
>
> I have been trying to use p.adjust() to do BH multiple test correction and
> have gotten some unexpected results. I thought that the equation for this was:
>
> pBH = p*n/i
Looking at the code for `p.adjust`, you see that the method is picked from a
switch function
lp <- length(p)
BH = {
i <- lp:1L
o <- order(p, decreasing = TRUE)
ro <- order(o)
pmin(1, cummin(n/i * p[o]))[ro]
}
You may not have sorted the p-values in pList.
>
> where p is the original p value, n is the number of tests and i is the rank
> of the p value. However when I try and recreate the corrected p from my most
> significant value it does not match up to the one computed by the method
> p.adjust:
>
>> setwd("C:/work/Methylation/IMA/GM/siteLists")
>>
>> hypTable <- read.delim("hypernormal vs others.txt")
>> pList <- hypTable$p
>> names(pList) <- hypTable$site
>>
>> adjusted <- p.adjust(pList, method="BH")
>> adjusted[1]
> cg27433479
> 0.05030589
>>
>> pList[1]*nrow(hypTable)/1
> cg27433479
> 0.09269194
>
No data provided, so unable to pursue this further.
> I tried to recreate this is a small example of a vector of 5 p values but
> everything worked as expected there. I was wondering if there is some subtle
> difference about how p.adjust operates? Is there something more complicated
> about how to calculate 'n' or 'i' - perhaps due to identical p values being
> assigned the same rank or something? Does anyone have an idea what might be
> going on here?
--
David Winsemius
Alameda, CA, USA
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
