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

----------

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?

> wilcox.test(passem\$temp ~ died30in)