Thanks for the hint. I tested this generic example and had the same behavior also in a special example, that can be found below. This example does not involve continuity correction but exibits the same "unexpected" behavior:
GDS_example = function() { print("---> GDS EXAMPLE <---", quote=F) #For GDS1331!!! exp_all_200601_at = c(281.600, 209.300, 202.600, 263.700, 242.500, 281.300, 329.100, 268.200, 311.700, 230.600, 280.600, 233.900, 270.200, 274.000, 169.000, 154.500, 160.100, 196.900, 169.500, 102.800, 148.600, 111.700, 117.000, 173.000, 79.100, 175.700, 149.500, 151.600, 78.700, 113.900, 196.200) c1 = c(1:14) c2 = c(20:31) group1 = exp_all_200601_at[c1] group2 = exp_all_200601_at[c2] wilcox_approx = wilcox.test(group1, group2, correct=F) # Deviation estimate -> N.B. the TRUE deviation value can NOT be known as we would have to sample _ALL_ possible "tissue" etc. #wilcox.test(group1, group2, conf.int=T)$estimate g1f = rep(1, length(c1)) g2f = rep(2, length(c2)) gf = as.factor(c(g1f, g2f)) exp_groups_200601_at = exp_all_200601_at[c(c1, c2)] exp_groups_data_frame = data.frame( exp_groups_200601_at, gf ) wilcox_exact = wilcox_test( exp_groups_200601_at ~gf, data = exp_groups_data_frame, conf.int=T, distribution="exact" ) print(wilcox_approx, quote=F) print(wilcox_exact) print(paste("P-value manually:", 2*(1 - pnorm(statistic(wilcox_exact)), quote=F) )) } Best, Cedric On Tuesday, 17. August 2010 17:50:48 Thomas Lumley wrote: > After fixing the parentheses in your code so it does run, it seems that the > difference is that wilcox.test defaults to using a continuity correction > and your manual calculation does not. Use wilcox.test(big1, big2, > correct=FALSE). > > -thomas > > On Tue, 17 Aug 2010, Cedric Laczny wrote: > > Hi, > > > > I became a little bit confused when working with the Wilcoxon test in R. > > As far as I understood, there are mainly two versions: > > 1) wilcox.test{stats}, which is the default and an approximation, > > especially, when ties are involved > > 2) wilcox_test{coin}, which does calculate the distribution _exactly_ > > even, with ties. > > > > I have the following scenario: > > > > #---BeginCode--- > > # big example > > size = 60 > > big1 = rnorm(size, 0, 1) > > big2 = rnorm(size, 0.5, 1 > > > > g1f = rep(1, size) > > g2f = rep(2, size) > > big = c(big1, big2) > > data_frame = data.frame(big, gr=as.factor(c(g1f, g2f))) > > > > wilcox_approx = wilcox.test(big1, big2) > > wilcox_exact = wilcox_test(big ~ gr, data=data_frame, > > distribution="exact") #---EndCode--- > > > > I found here > > http://www-stat.stanford.edu/~susan/courses/s141/hononpara.pdf that > > wilcox.test (at least for the signed rank test) relies on exact > > (p-)values until N = n1 + n2 = 50. > > I can reproduce this, when using e.g. size = 15. The p-values then are > > the same, as I would expect it, having read the info from the link. > > > > #---BeginCode--- > > print(paste("Wilcox approx p-value:", wilcox_approx$p.value), quote=F) > > print(paste("Wilcox exact p-value:", pvalue(wilcox_exact)), quote=F) > > #---EndCode--- > > > > That said, if I set e.g. size = 60, then the p-values of wilcox.test and > > wilcox_test differ, as expected. > > > > What's puzzling me particularly is the differing results when wanting to > > calculate the p-value manually, for bigger sample sizes. > > > > So, if we get the W-score from wilcox.test: > > > > #---BeginCode--- > > W_big = wilcox.test(big1, big2))$statistic > > #---EndCode--- > > > > and "convert" it to a Z-score, like this: > > > > #---BeginCode--- > > mu_big = (size^2)/2 > > sd_big = sqrt(size*size*(size + size + 1)/12) > > N = size + size > > sd_big_corr = sqrt( (size * size) / (N * (N - 1)) * (N^3 - N) / 12 ) > > > > Z_big = (((W_big - mu_big)/sd_big) > > #---EndCode--- > > > > The Z-Score (Z_big) is equal to the statistic of wilcox_test. > > So far so good. And now comes the main problem. > > When I follow the documentation correctly, the p-value for a given > > W-score/- statistic ist calculated using the normal-approximation with > > the Z-score. However, when I do that, I get a different result than what > > I would expect. Because I would expect the p-value of wilcox.test to be > > equal to 2*pnorm(Z_big), which is in fact _not_ equal. Please see: > > > > #---BeginCode--- > > p_value_manual = 2 * pnorm(Z_big) > > > > print("--- Resulting pvalues --- ", quote=F) > > print(paste("Wilcox approx p-value:", wilcox_approx$p.value), quote=F) > > print(paste("Wilcox exact p-value:", pvalue(wilcox_exact)), quote=F) > > print(paste("P-value manual:", p_value_manual), quote=F) > > #---EndCode--- > > > > So how is the calculation of the p-value performed in wilcox.test, when > > the sample sizes are big? Because this might explain why the value > > differs from that being calculated manually. > > > > Best regards, > > > > Cedric > > > > ______________________________________________ > > R-help@r-project.org mailing list > > https://stat.ethz.ch/mailman/listinfo/r-help > > PLEASE do read the posting guide > > http://www.R-project.org/posting-guide.html and provide commented, > > minimal, self-contained, reproducible code. > > Thomas Lumley > Professor of Biostatistics > University of Washington, Seattle ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.