0
$\begingroup$

I'm hoping someone can shed some light on the discrepancy between the p-values obtained from two different functions in R for conducting a spearman's test. The two tests produce the same p-values in some circumstances but not others. I think I've identified what those circumstances are. When there are repeated values in the vector (i.e. ties), the two procedures produce the same results but when repeated values are replaced the two results of the two procedures diverge. See code below for example:

x2 <- c(1,2,3,4,5,6,7) 
y2 <- c(1,3,6,2,7,4,6) 
cor.test(x2, y2, method="spearman")

Spearman's rank correlation rho

data:  x2 and y2
S = 21.692, p-value = 0.1436
alternative hypothesis: true rho is not equal to 0
sample estimates:
      rho 
0.6126375 

Note the p-value of 0.1436. Using a different function for spearman's from the Hmisc library:

rcorr(x2,y2,type="spearman")
     x    y
x 1.00 0.61
y 0.61 1.00

n= 7 


P
  x      y     
x        0.1436
y 0.1436  

Again, note the p-value of 0.1436. These two procedures result in the same values -as expected of course. However, if I alter the one of the orignal vectors by replacing at least one value that is repeated, the results of the two procedures diverge:

x2 <- c(1,2,3,4,5,6,7) 
y2 <- c(1,3,6,2,7,4,9) #changed last element from 6 to 9 

rcorr(x2,y2,type="spearman")
     x    y
x 1.00 0.75
y 0.75 1.00

n= 7 


P
  x      y     
x        0.0522
y 0.0522  

Note the p-value using cor.test() from above is 0.0522 but using rcorr() I get a p-value of 0.06627 below.

cor.test(x2, y2, method="spearman")

    Spearman's rank correlation rho

data:  x2 and y2
S = 14, p-value = 0.06627
alternative hypothesis: true rho is not equal to 0
sample estimates:
 rho 
0.75 

Now say I replace that last value in y2 with 7 instead of 9. Keep in mind, that the 7 in y2 is repeated twice (in contrast to 9, which is not found in any of the original vectors and thus its addition does not produce a tie).

x2 <- c(1,2,3,4,5,6,7) 
y2 <- c(1,3,6,2,7,4,7) #This time replaced 3 with 7 instead of 9 -note 7 repeats but 9 did not
cor.test(x2, y2, method="spearman")

    Spearman's rank correlation rho

data:  x2 and y2
S = 15.638, p-value = 0.06763
alternative hypothesis: true rho is not equal to 0
sample estimates:
    rho 
0.72075 

Note the p-value of 0.06763. In this case, cor.test() and rcorr() produce the same p-value:

rcorr(x2,y2,type="spearman")
     x    y
x 1.00 0.72
y 0.72 1.00

n= 7 


P
  x      y     
x        0.0676
y 0.0676  

Does anyone have an idea of what's going on here? Why is it that when tied/repeated values are present, the two methods produce the same results yet diverge when additional non tied values are added.

$\endgroup$
2
$\begingroup$

? cor.test() says:

For Spearman's test, p-values are computed using algorithm AS 89 for n < 1290 and exact = TRUE, otherwise via the asymptotic t approximation.

(and default value is exact = TRUE). And its Warning message:

In cor.test.default(x2, y2, method = "spearman") : Cannot compute exact p-value with ties

They mean that if you don't give an argument exact, cor.test() uses algorithm AS 89 when n < 1290 and no tie. rcorr() always calculates p-values via the asymptotic t approximation, so cor.test() and rcorr() return the same p-value when any ties are present.

getS3method("cor.test", "default") gives you the code. D. J. Best and D. E. Roberts, Algorithm AS 89: The Upper Tail Probabilities of Spearman's Rho, Applied Statistics, Vol.24, No.3 1975, pp.377-379 gives the details about AS 89.

[Edited]

Algorithm AS 89 has two methods, the "exact" method and the "semi-exact" method. cor.test() uses the former when n < 10 and the latter when 9 < n < 1290. In your case, p.value was calculated by the "exact" method.

If you do it manually;

library(e1071)           # to get permutations

x2 <- c(1,2,3,4,5,6,7) 
y2 <- c(1,3,6,2,7,4,9)   # second y2
permutation <- permutations(7)
n <- length(x2)

rho <- cor.test(x2, y2, method="spearman")$estimate  # 0.75

# the function to calculate rho between an argument and y2
f.rho <- function(a) 1 - 6 * sum((rank(a) - rank(y2))^2) / (n^3 - n)

x2.all.permutation <- matrix(x2[permutation], ncol=7)
x2.all.permutation.rho <- apply(x2.all.permutation, 1, function(a) f.rho(a))

sum(x2.all.permutation.rho > rho) / factorial(7) * 2   # [1] 0.06626984

cor.test(x2, y2, method="spearman")$p.value   # [1] 0.06626984  # the same
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.