Wilcoxon秩-和检验结果的解读中遇到问题

areg

下面示例数据来源自?wilcox.test的帮助文件中的示例数据：

> x <- c(0.80, 0.83, 1.89, 1.04, 1.45, 1.38, 1.91, 1.64, 0.73, 1.46)

> y <- c(1.15, 0.88, 0.90, 0.74, 1.21)

> wilcox.test(x, y)

        Wilcoxon rank sum test

data:  x and y

W = 35, p-value = 0.2544

alternative hypothesis: true location shift is not equal to 0

######################################################

为理解W到底表示什么，按Wilcoxon秩-和检验公式，一步步分解如下，对x及y进行混和编秩，得较小例数秩为：2,5,6,8,9

（1）在上述结果中，W=35表示什么？怎么来的？

（2）Ry<-sum(2,5,6,8,9)        #求较例数的秩和，是30

（3）ERy<-5*(5+10+1)/2        #表示秩和期望，为40

（4）VarRy<-5*10*(5+10+1)/12                         #方差，为66.66667

（5）T<-(abs(Ry-ERy)-1/2)/VarRy^(1/2)    #等于 1.163508

（6）p_v<-2*(1-pnorm(T))    #求精确p-值，等于0.2446236

（7）最后结果与使用wilcox.test(x,y)求的p-值近似一致，按理使用wilcox.test后求得的W应该是较少例数的样本的秩和，本例中它应该是30，问题出在此处。

minhaosheng

[s:14] 这个问题应该是由于R中采用的W统计量的计算方法，并不是你上面所述的那个标准方法，我在rseek搜到wilcox.test函数的代码，

在这里svn.r-project.org/R/trunk/src/library/stats/R/wilcox.test.R

其中有这么一段

else { ##-------------------------- 2-sample case ---------------------------

        if(length(y) < 1L)

            stop("not enough 'y' observations")

        METHOD <- "Wilcoxon rank sum test"

        r <- rank(c(x - mu, y))

        n.x <- as.double(length(x))

        n.y <- as.double(length(y))

        if(is.null(exact))

            exact <- (n.x < 50) && (n.y < 50)

        STATISTIC <- sum(r[seq_along(x)]) - n.x * (n.x + 1) / 2

        names(STATISTIC) <- "W"

采用你的数据，再利用这个代码，可以计算出W＝35,但是为什么这么算，我还没想清楚 [s:16]

areg

谢谢

我也刚弄明白，回来一看，有你回复。

上面我是把统计量弄错了，我是按一般资料上的求T统计量，使用的是正态近似方法。

而在wilcox.test中使用的是W统计量

W统计量为：第1组的秩和减去第1组秩和的数学期望值。

x <- c(0.80, 0.83, 1.89, 1.04, 1.45, 1.38, 1.91, 1.64, 0.73, 1.46)

y <- c(1.15, 0.88, 0.90, 0.74, 1.21)

r <- rank(c(x, y))

n.x <- as.double(length(x))

n.y <- as.double(length(y))

STATISTIC <- sum(r[seq_along(x)]) - n.x * (n.x + 1) / 2

names(STATISTIC) <- "W"

if(STATISTIC >(n.x*n.y/2)) p<-pwilcox(STATISTIC-1, n.x, n.y,lower.tail=FALSE) else p<-pwilcox(STATISTIC, n.x, n.y)

min(2 * p, 1)

> sum(r[seq_along(x)])

[1] 90

> n.x * (n.x + 1) / 2

[1] 55

上面的35就是90-55得到的。

一般资料上的第一组期望秩和恒等于E（R1）＝n1*(n1+n2+1)/2，是两独立样本（即两组）；在没有“结”时，单个样本的数学期望是n(n+1)/4。

为什么上面减的是n.x * (n.x + 1) / 2，而不n.x * (n.x + n.y+1) / 2，这点还是没有弄明白。

[s:12]

foshuochanyu

同关注，我的想法是：

1.关于是否用正态估计

楼主用函数wilcox.test()进行wilcoxn秩和检验，使用

========================================

> wilcox.test(x, y)

========================================

进行的是精确检验，即没有正态估计，因此得出的p-value为

========================================

p-value=0.2544

========================================

与接下来验证的(且不说是否正确)

=============================================================

（5）T<-(abs(Ry-ERy)-1/2)/VarRy^(1/2)    #等于 1.163508

（6）p_v<-2*(1-pnorm(T))    #求精确p-值，等于0.2446236

=============================================================

的p-value不一致是正常的。

而采用标准正态估计为（参考R-help）：

> wilcox.test(x, y, exact = FALSE, correct = FALSE)$p.value

[1] 0.2206714

2 关于统计量U（或者W）

维基百科的介绍（参考网址<http://en.wikipedia.org/wiki/Mann-Whitney_U_test> --->尽管图被阉割了）

统计量U的计算：

U.x = R.x-n.x*(n.x-1)/2=90-45=35

＃R.x为第一组的秩和

这里也提出如果样本数较多，统计量U近似为正态分布

3 关于近似正态分布的参数

不清楚楼主是参考哪些资料（可能是wilcox的原始论文，因为据说wilcox和Mann-Whitney的两篇文章中使用的检验方法实质虽然相同，但细节有异。我也再仔细看看）得出期望E的数学表达式

＃可能楼主参考文献中此处不是统计量U而是其他统计量

Mann-Whiney文章中（参考文章On a Test of Whether one of Two Random Variables is Stochastically Larger than the Other. H. B. Mann and D. R. Whitney）指出其正态分布期望为 n.x*n.y/2, 方差为 n.x*n.y*（n.x+n.y+1）/12

因此可做修改楼主验证为：

> ERy <- 5*10/2    #方差

> VarRy <- 5*10*(5+10+1)/12   ＃期望

> T <- (35-ERy)/VarRy^(1/2)   #标准化,此处用第一组的秩和，也可用第二组的秩和为15

> 2*(1-pnorm(T))

[1] 0.2206714

结果与使用

> wilcox.test(x, y, exact = FALSE, correct = FALSE)$p.value

[1] 0.2206714

一致

foshuochanyu

再回自己的回复（MS我很无聊） [s:12]

对于期望的问题，有了些新的认识。

楼主开始使用的E = n.x*(n.x+n.y+1)/2对应的统计量是秩和R而不是统计量U（也就是R运算结果的出现的W=35），所以用的期望也不一样。

这也解决了一个统计检验方法两个名字的问题：

wilxcon秩和检验——用的统计量是秩和R；

Mann-Whitney U检验——用的统计量是U

对于一组数据R和U的换算关系是U=R-n.x*(n.x+1)/2（推测）

areg

谢谢foshuochanyu

由于平时很少用非参检验，现在在假期，就找来看看相关的资料，由于现在分析主要用R。一般资料上都是使用近似正态，使用“T”统计量，看到有的资料上使用的统计量为“W”（用SAS分析的，还有用minitab分析的结果中统计量也是用W表示，mini中其实那个“W”就是例数较少的组的秩和），自己就试用R中的函数来练一下手................

结果是自己看资料不细心，弄出不少笑话。

点开R中wilcox.test的帮助文件，用里面的示例数据学习，但竟然把Wilcoxon秩-和（双侧，正态近似）的条件都忽略了（要求：n1及n2均不少于10，观察值变量有连续性分布），而我使用帮助文件中示例数据y只有5例 [s:12]

这下明白了，（1）自己先把统计量弄糊涂了；（2）minitable等的wilcoxon检验结果中的W与R环境中的W不是一回事。

在R环境中W=Rix-n.x*(n.x+1)/2 #某组的秩和减去对应组样本量nx(nx+1)除以2。你的推测正确，R环境中wilcox.test函数源程序就是这样编写的（我使用的onpenSUSE系统，查看源程序方便）

对于用秩和差事去的那个量n.x*(n.x+1)/2而不是n.x * (n.x + n.y+1) / 2，（此为近似检验中用秩和去减的量）。

以上问题，可能是由于我看的资料比较杂，而不同资料的表述有所不同造成的。谢谢两位的指点，让我对非参检验有了较清晰的了解。

areg

谢谢楼上两位，刚才我提交了回复，使用了一个“个性表情”，回复不知道跑哪去了 [s:12]

现在基本明白了R环境中的wilcoxon检验，以后抽空再练习一下。 [s:11]

foshuochanyu

恩，我也看了一下午，把以前的资料也翻出来仔细看看。

以后多交流。

minhaosheng

[quote]引用第6楼areg于2010-02-05 19:42发表的回 4楼(foshuochanyu) 的帖子 :

谢谢楼上两位，刚才我提交了回复，使用了一个“个性表情”，回复不知道跑哪去了 [s:12]

现在基本明白了R环境中的wilcoxon检验，以后抽空再练习一下。 [s:11]

[/quote]

[s:12]和我前天一样，回复里有一个超链接，也不知道回复到哪儿去了。[s:14]

Ihavenothing

为了防止广告，超链接会在后台审核，所以可能不是很及时。

areg

哈哈，明白了，以后就用最简单的形式，就不会出现这类问题了。

好久没来论坛，各方面的管理有新改进啊，争取尽快适应。

谢谢

Ihavenothing

看来也是一位老前辈啊。[s:11]

areg

好久没弄了，都忘记得差不多啦，就算从头再来吧，还各位多帮一下老兄 [s:11]

Alluvion

请问我在采用wilcox.test进行两个成对数据的秩序是否相同时，是否需要剔除NA值的成对样本？先行谢过……[s:19]

funpipi

请问对于n1=3,n2=3的情况，是否可以使用R环境下的wilcoxon检验？

我没有找到帮助文件中的“要求：n1及n2均不少于10”

salen83123

我也是遇到n1<10,n2>200，直接使用wilcox.test()恐怕不对，应该怎么办呢？