Different P value of Mann-Whitney-Wilcoxon between R and Stata I want to compare the median of body temperature (passem$temp) across two groups (whether the patient died or not, died30in.numeric). I ran Mann Whitney Wilcoxon first in R using this command wilcox.test(passem$temp, died30in.numeric) and got these results:
> wilcox.test(passem$temp, died30in.numeric)

    Wilcoxon rank sum test with continuity correction

data:  passem$temp and died30in.numeric
W = 10719076, p-value < 2.2e-16
alternative hypothesis: true location shift is not equal to 0

However, when I ran the test in Stata, I got the following results which were different from that in R:
ranksum temp, by(died30in)

Two-sample Wilcoxon rank-sum (Mann–Whitney) test

    died30in |      Obs    Rank sum    Expected
-------------+---------------------------------
           0 |     2987   4924500.5   4891212.5
           1 |      287    436674.5    469962.5
-------------+---------------------------------
    Combined |     3274     5361175     5361175

Unadjusted variance   2.340e+08
Adjustment for ties  -746907.26
                     ----------
Adjusted variance     2.332e+08

H0: temp(died30in==0) = temp(died30in==1)
         z =  2.180
Prob > |z| = 0.0293

I also double check with SPSS, and the results were concordant with Stata analysis.
Of note, it might be an oversimplification (or misconception by some) to assume that Mann-Whitney-Wilcoxon is used to compare two medians; still, I want to conduct the test on our data.
What is the reason for this unexplained difference in results?

The reason for the dramatic difference is that you have run the Wilcoxon test incorrectly in R. You probably intended:
> wilcox.test(passem$temp ~ died30in)

