R devel - Feb 2006 - achieved.alpha calculated wrong in wilcox.test (PR#8557)

In R-2.2.1 stable the file wilcox.test.R line 86 has

    achieved.alpha<-2*psignrank(trunc(qu),n)

and should have

    achieved.alpha<-2*psignrank(trunc(qu)-1,n)

this is apparently a thinko not a typo so similar statements are probably
wrong too (line 97, line 109, line 293, line 304, line 316).

Reference: Hollander and Wolfe or

    http://www.stat.umn.edu/geyer/5601/examp/signrank.html
    http://www.stat.umn.edu/geyer/5601/examp/ranksum.html

If the signed rank c. i. is diffs[k] to diffs[length(diffs) + 1 - k],
then the probability of failure to cover for a two-sided interval is

    2 * Pr(T < k)

where T is a random variable having the null distribution of the test
statistic.  It is easy to check this is true in the k = 1 case.  The
confidence interval fails to cover if and only if ALL of the data points
are to one side of the true unknown parameter value, and the probability
of this happening is Pr(T = 0), not Pr(T <= 1) as the code on line 86 has it.

This leads to bizarre behavior.  Consider the following homework problem.

---------- begin R output ----------
> X <-
read.table(url("http://www.stat.umn.edu/geyer/5601/hwdata/t3-3.txt"),
+     header = TRUE)
> attach(X)
> wilcox.test(y, x, paired = TRUE, conf.int = TRUE)
        Wilcoxon signed rank test

data:  y and x
V = 24, p-value = 0.1094
alternative hypothesis: true mu is not equal to 0
92.2 percent confidence interval:
 -25 605
sample estimates:
(pseudo)median
            80

Warning message:
Requested conf.level not achievable in: switch(alternative, two.sided = {
----------- end R output -----------

when the correct behavior is given by the code on

    http://www.stat.umn.edu/geyer/5601/examp/signrank.html#conf

---------- begin R output ----------
> conf.level <- 0.95
> z <- sort(y - x)
> n <- length(z)
> walsh <- outer(z, z, "+") / 2
> walsh <- sort(walsh[lower.tri(walsh, diag = TRUE)])
> m <- length(walsh)
> alpha <- 1 - conf.level
> k <- qsignrank(alpha / 2, n)
> if (k == 0) k <- k + 1
> print(k)
[1] 3
> cat("achieved confidence level:", 1 - 2 * psignrank(k - 1, n),
"\n")
achieved confidence level: 0.953125
> c(walsh[k], walsh[m + 1 - k])
[1] -25 605
----------- end R output -----------

Note that wilcox.test reports the correct interval, which
is diffs[k] to diffs[length(diffs) + 1 - k] with k == 3.
Its mistake is to claim that this is only
a 92.2 percent confidence interval when it is
actually a 95.3125 percent confidence interval.

IMHO it is also a bug to give no option to get the actual achieved confidence
level unless the level requested is not achievable and then only in the text
of a warning.  But that is debatable.

-- 
Charles Geyer
Professor, School of Statistics
University of Minnesota
charlie at stat.umn.edu