[Rd] p.adjust(<NA>s), was "Re: [BioC] limma and p-values"

[Gordon Smyth]

Martin,

Thanks for this.  I was planning to code suggestions for p.adjust() myself 
today, but you got in before me.

I don't think that the new way of handling NAs is quite correct.  The 
trouble is that all of the adjustment methods other than "none" and 
"bonferonni" are step-down or step-up methods, meaning that you need to 
know all the p-values before you can adjust any of them.

Suppose for example that you have a single p-value 'p=0.05' and you want to 
adjust this using Holm's method with 'n=10'.  Holm's method is a step-down 
Bonferonni method so the adjusted p-value can be anywhere between 0.05 and 
0.5 depending on what the other nine p-values are.  So the correct adjusted 
p-value is indeterminate.

The value returned by 'p.adjust(0.05,method="holm",n=10)' is 0.5.  The same 
value is returned if 'p' is a vector with nine NAs:

   p <- c(0.05,rep(NA,9))
   p.adjust(p,method="holm")

In effect the nine NA p-values are treated as if they are equal to one.  So 
rather than accommodating unknown data, p.adjust() is in effect making data 
up.  I think it would be best not to do this but, if you do, there should 
at least be a warning message like

"warning: p-values indeterminate in the presence of NAs, most conservative 
possible values being returned"

I think that the correct treatment is that implemented for 'na.rm=TRUE', 
i.e., the non-NA adjusted p-values should be the same as if there were no 
NAs in 'p'.  This is because NAs in 'p' virtually always indicate, not an 
existing but unknown p-value, but a p-value which could not be computed 
because there was insufficient data to do so.  This means that there is 
zero probability of rejecting the null hypothesis for these 
cases.  Therefore it is unnecessary to adjust the other p-values for these 
cases.

The argument 'n' to 'p.adjust()' is problematic for the same reasons that 
NAs in 'p' are.  Setting 'n' to be less than than sum(!is.na(p)) should I 
think return an error.  Setting 'n' to be more than sum(!is.na(p)) makes 
all the adjusted p-values indeterminate unless method="none" or 
method="bonferroni".  And setting 'n' to be between sum(!is.na(p)) and 
length(p) when these are different is just horrible.  Would anything be 
lost if this argument was removed?

Gordon

At 07:29 AM 9/01/2005, Martin Maechler wrote:
>I've thought more and made experiements with R code versions
>and just now committed a new version of  p.adjust()  to R-devel
>--> https://svn.r-project.org/R/trunk/src/library/stats/R/p.adjust.R
>which does sensible NA handling by default and
>*additionally* has an "na.rm" argument (set to FALSE by
>default).  The extended 'Examples' secion on the help page
>     https://svn.r-project.org/R/trunk/src/library/stats/man/p.adjust.Rd
>shows how the new NA handling is typically much more sensible
>than using "na.rm = TRUE".
>
>Martin
>
>
> >>>>> "MM" == Martin Maechler <maechler at stat.math.ethz.ch>
> >>>>>     on Sat, 8 Jan 2005 17:19:23 +0100 writes:
>
> >>>>> "GS" == Gordon K Smyth <smyth at wehi.edu.au>
> >>>>>     on Sat, 8 Jan 2005 01:11:30 +1100 (EST) writes:
>
>     MM>     <.............>
>
>     GS> p.adjust() unfortunately gives incorrect results when
>     GS> 'p' includes NAs.  The results from topTable are
>     GS> correct.  topTable() takes care to remove NAs before
>     GS> passing the values to p.adjust().
>
>     MM> There's at least one bug in p.adjust(): The "hommel"
>     MM> method currently does not work at all with NAs (and I
>     MM> have an uncommitted fix ready for this bug).  OTOH, the
>     MM> current version of p.adjust() ``works'' with NA's, apart
>     MM> from Hommel's method, but by using "n = length(p)" in
>     MM> the correction formulae, i.e. *including* the NAs for
>     MM> determining sample size `n' {my fix to "hommel" would do
>     MM> this as well}.
>
>     MM> My question is what p.adjust() should do when there are
>     MM> NA's more generally, or more specifically which `n' to
>     MM> use in the correction formula. Your proposal amounts to
>     MM> ``drop NA's and forget about them till the very end''
>     MM> (where they are wanted in the result), i.e., your sample
>     MM> size `n' would be sum(!is.na(p)) instead of length(p).
>
>     MM> To me it doesn't seem obvious that this setting "n =
>     MM> #{non-NA observations}" is desirable for all P-value
>     MM> adjustment methods. One argument for keeping ``n = #{all
>     MM> observations}'' at least for some correction methods is
>     MM> the following "continuity" one:
>
>     MM> If only a few ``irrelevant'' (let's say > 0.5) P-values
>     MM> are replaced by NA, the adjusted relevant small P-values
>     MM> shouldn't change much, ideally not at all.  I'm really
>     MM> no scholar on this topic, but e.g. for "holm" I think I
>     MM> would want to keep ``full n'' because of the above
>     MM> continuity argument.  BTW, for "fdr", I don't see a
>     MM> straightforward way to achieve the desired continuity.
>     MM> 5D Of course, p.adjust() could adopt the possibility of
>     MM> chosing how NA's should be treated e.g. by another
>     MM> argument ``use.na = TRUE/FALSE'' and hence allow both
>     MM> versions.
>
>     MM> Feedback very welcome, particularly from ``P-value
>     MM> experts'' ;-)
>
>     MM> Martin Maechler, ETH Zurich