Martin Maechler
on Thu, 12 Dec 2019 17:20:47 +0100 writes:
Karolis Koncevičius
on Mon, 9 Dec 2019 23:43:36 +0200 writes:
>> So I tried adding Infinity support for all cases. And it
>> is (as could be expected) more complicated than I
>> thought.
> "Of course !" Thank you, Karolis, in any case!
>> It is easy to add Inf support for the test. The problems
>> start with conf.int=TRUE.
>> Currently confidence intervals are computed via
>> `uniroot()` and, in the case of infinities, we are
>> computationally looking for roots over infinite interval
>> which results in an error. I suspect this is the reason
>> Inf values were removed in the first place.
> Maybe. It's still wrong to be done "up front". I'm sure
> 98% (or so ;-) of all calls to wilcox.test() do *not* use
> conf.int = TRUE
>> Just a note, I found a few more errors/inconsistencies
>> when requesting confidence intervals with paired=TRUE
>> (due to Infinities being left in).
>> Current error in Inf-Inf scenario:
>> wilcox.test(c(1,2,Inf), c(4,8,Inf), paired=TRUE,
>> conf.int=TRUE) Error in if (ZEROES) x <- x[x != 0] :
>> missing value where TRUE/FALSE needed
> Good catch .. notably as it also happens when
> conf.int=FALSE as by default. My version of wilcox.test()
> now does give the same as when the to 'Inf' are left away.
>> NaN confidence intervals:
>> wilcox.test(c(1:9,Inf), c(21:28,Inf,30), paired=TRUE,
>> conf.int=TRUE)
>> Wilcoxon signed rank test with continuity correction
>> data: c(1:9, Inf) and c(21:28, Inf, 30) V = 9.5, p-value
>> = 0.0586 alternative hypothesis: true location shift is
>> not equal to 0 0 percent confidence interval: NaN NaN
>> sample estimates: midrange NaN
> I don't see a big problem here. The NaN's in some sense
> show the best that can be computed with this data.
> Statistic and P-value, but no conf.int.
>> The easiest "fix" for consistency would be to simply
>> remove Infinity support from the paired=TRUE case.
> I strongly disagree. We are not pruning good
> functionality just for some definition of consistency.
>> But going with the more desirable approach of adding
>> Infinity support for non-paired cases - it is currently
>> not clear to me what confidence intervals and
>> pseudomedian should be. Specially when Infinities are on
>> both sides.
> I deem that not to be a big deal. They are not defined
> given the default formulas and that is reflected by NA /
> NaN in those parts of the result.
>> Regards, Karolis Koncevičius.
> But I have also spent a few hours now on
> wilcox.test.default() behavior notably also looking at the
> "rounding" / "machine precision" situation, and also on
> your remark that the 'method: ...' does not indicate well
> enough what was computed.
> In my (not yet committed) but hereby proposed enhancement
> of wilcox.test(), I have a new argument, 'digits.rank =
> Inf' (the default 'Inf' corresponding to the current
> behavior) with help page documentation:
> digits.rank: a number; if finite, ‘rank(signif(r,
> digits.rank))’ will be used to compute ranks for the test
> statistic instead of (the default) ‘rank(r)’.
> and then in 'Details :'
> For stability reasons, it may be advisable to use
> rounded data or to set ‘digits.rank = 7’, say, such that
> determination of ties does not depend on very small
> numeric differences (see the example).
> and then in 'Examples: '
> ## accuracy in ties determination via 'digits.rank':
> wilcox.test( 4:2, 3:1, paired=TRUE) # Warning: cannot
> compute exact p-value with ties wilcox.test((4:2)/10,
> (3:1)/10, paired=TRUE) # no ties => *no* warning
> wilcox.test((4:2)/10, (3:1)/10, paired=TRUE, digits.rank =
> 9) # same ties as (4:2, 3:1)
> ----------------------
> Lastly, I propose to replace "test" by "exact test" in the
> 'method' component (and print out) in case exact
> computations were used. This information should be part
> of the returned "htest" object, and not only visible from
> the arguments and warnings that are printed during the
> computations. This last change is in some sense the "most
> back-incompatible" change of these, because many
> wilcox.test() printouts would slightly change, e.g.,
>> w0 <- wilcox.test( 1:5, 4*(0:4), paired=TRUE)
> Wilcoxon signed rank exact test
> data: 1:5 and 4 * (0:4) V = 1, p-value = 0.125
> alternative hypothesis: true location shift is not equal
> to 0
> where before (in R <= 3.6.x) it is just
> Wilcoxon signed rank test
> data: ......... ............... ...............
> but I think R 4.0.0 is a good occasion for such a change.
> Constructive feedback on all this is very welcome! Martin
... none ... I "assume" this means everybody likes the idea ;-)
Anyway, now comitted to R-devel (for R 4.0.0), svn rev 77569
(in 'NEW FEATURES').
Martin
>> On 2019-12-07 23:18, Karolis Koncevičius wrote:
>>> Thank you for a fast response. Nice to see this mailing
>>> list being so alive.
>>>
>>> Regarding Inf issue: I agree with your assessment that
>>> Inf should not be removed. The code gave me an
>>> impression that Inf values were intentionally removed
>>> (since is.finite() was used everywhere, except for
>>> paired case). I will try to adjust my patch according to
>>> your feedback.
>>>
>>> One more thing: it seems like you assumed that issues
>>> 2:4 are all related to machine precision, which is not
>>> the case - only 2nd issue is. Just wanted to draw this
>>> to your attention, in case you might have some feedback
>>> and guidelines about issues 3 and 4 as well.
>>>
>>>
>>>
>>> On 2019-12-07 21:59, Martin Maechler wrote:
>>>>>>>>> Karolis Koncevičius on Sat, 7 Dec 2019 20:55:36
>>>>>>>>> +0200 writes:
>>>>
>>>> > Hello, > Writing to share some things I've found
>>>> about wilcox.test() that seem a > a bit inconsistent.
>>>>
>>>> > 1. Inf values are not removed if paired=TRUE
>>>>
>>>> > # returns different results (Inf is removed): >
>>>> wilcox.test(c(1,2,3,4), c(0,9,8,7)) >
>>>> wilcox.test(c(1,2,3,4), c(0,9,8,Inf))
>>>>
>>>> > # returns the same result (Inf is left as value with
>>>> highest rank): > wilcox.test(c(1,2,3,4), c(0,9,8,7),
>>>> paired=TRUE) > wilcox.test(c(1,2,3,4), c(0,9,8,Inf),
>>>> paired=TRUE)
>>>>
>>>> Now which of the two cases do you consider correct ?
>>>>
>>>> IHMO, the 2nd one is correct: it is exactly one
>>>> property of the *robustness* of the wilcoxon test and
>>>> very desirable that any (positive) outlier is treated
>>>> the same as just "the largest value" and the test
>>>> statistic (and hence the p-value) should not change.
>>>>
>>>> So I think the first case is wrong, notably if
>>>> modified, (such that the last y is the overall maximal
>>>> value (slightly larger sample):
>>>>
>>>>> wilcox.test(1:7, 1/8+ c(9:4, 12))
>>>>
>>>> Wilcoxon rank sum test
>>>>
>>>> data: 1:7 and 1/8 + c(9:4, 12) W = 6, p-value = 0.01748
>>>> alternative hypothesis: true location shift is not
>>>> equal to 0
>>>>
>>>>> wilcox.test(1:7, 1/8+ c(9:4, 10000))
>>>>
>>>> Wilcoxon rank sum test
>>>>
>>>> data: 1:7 and 1/8 + c(9:4, 10000) W = 6, p-value =
>>>> 0.01748 alternative hypothesis: true location shift is
>>>> not equal to 0
>>>>
>>>>> wilcox.test(1:7, 1/8+ c(9:4, Inf))
>>>>
>>>> Wilcoxon rank sum test
>>>>
>>>> data: 1:7 and 1/8 + c(9:4, Inf) W = 6, p-value =
>>>> 0.03497 alternative hypothesis: true location shift is
>>>> not equal to 0
>>>>
>>>> The Inf case should definitely give the same as the
>>>> 10'000 case. That's exactly one property of a robust
>>>> statistic.
>>>>
>>>> Thank you, Karolis, this is pretty embarrassing to only
>>>> be detected now after 25+ years of R in use ...
>>>>
>>>> The correct fix starts with replacing the is.finite()
>>>> by !is.na() and keep the 'Inf' in the rank
>>>> computations... (but then probably also deal with the
>>>> case of more than one Inf, notably the Inf - Inf
>>>> "exception" which is not triggered by your example...)
>>>>
>>>>
>>>> ---
>>>>
>>>> Ben addressed the "rounding" / numerical issues
>>>> unavoidable for the other problems.
>>>>
>>>> > 2. tolerance issues with paired=TRUE.
>>>>
>>>> > wilcox.test(c(4, 3, 2), c(3, 2, 1), paired=TRUE) > #
>>>> ... > # Warning: cannot compute exact p-value with
>>>> ties
>>>>
>>>> > wilcox.test(c(0.4,0.3,0.2), c(0.3,0.2,0.1),
>>>> paired=TRUE) > # ... > # no warning
>>>>
>>>> > 3. Always 'x observations are missing' when
>>>> paired=TRUE
>>>>
>>>> > wilcox.test(c(1,2), c(NA_integer_,NA_integer_),
>>>> paired=TRUE) > # ... > # Error: not enough (finite)
>>>> 'x' observations
>>>>
>>>> > 4. No indication if normal approximation was used:
>>>>
>>>> > # different numbers, but same "method" name >
>>>> wilcox.test(rnorm(10), exact=FALSE, correct=FALSE) >
>>>> wilcox.test(rnorm(10), exact=TRUE, correct=FALSE)
>>>>
>>>>
>>>> > From all of these I am pretty sure the 1st one is
>>>> likely unintended, > so attaching a small patch to
>>>> adjust it. Can also try patching others if > consensus
>>>> is reached that the behavioiur has to be modified.
>>>>
>>>> > Kind regards, > Karolis Koncevičius.
>>>>
