[Rd] achieved.alpha calculated wrong in wilcox.test (PR#8557)
charlie@stat.umn.edu
charlie at stat.umn.edu
Thu Feb 2 23:26:28 CET 2006
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
More information about the R-devel
mailing list