From: <charlie_at_stat.umn.edu>

Date: Thu 02 Feb 2006 - 22:26:28 GMT

Date: Thu 02 Feb 2006 - 22:26:28 GMT

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@stat.umn.edu ______________________________________________ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-develReceived on Fri Feb 03 09:47:10 2006

*
This archive was generated by hypermail 2.1.8
: Fri 03 Feb 2006 - 07:01:47 GMT
*