10
$\begingroup$

Is there a function to test the hypothesis that the correlation of two vectors is equal to a given number, say 0.75? Using cor.test I can test cor=0 and I can see whether 0.75 is inside the confidence interval. But is there a function to compute the p-value for cor=0.75?

x <- rnorm(10)
y <- x+rnorm(10)
cor.test(x, y)
$\endgroup$

migrated from stackoverflow.com Aug 13 '11 at 16:56

This question came from our site for professional and enthusiast programmers.

  • 2
    $\begingroup$ This question is better suited for crossvalidated.com $\endgroup$ – Sacha Epskamp Aug 13 '11 at 10:34
  • 1
    $\begingroup$ @sacha - please check a site's FAQ first, the stats.se site faq recommends that programming questions using R are posted on SO. $\endgroup$ – Kev Aug 13 '11 at 12:03
  • $\begingroup$ The question "is there a function to compute the p-value for cor=0.75?" has nothing to do with programming. It is a statistical question. $\endgroup$ – Sacha Epskamp Aug 13 '11 at 12:52
  • $\begingroup$ I will consult the stats folks and see what they think. $\endgroup$ – Kev Aug 13 '11 at 12:56
  • 1
    $\begingroup$ @mosaic Please, register your account here. This way, you will be able to associate your SO account with the present one. $\endgroup$ – chl Aug 13 '11 at 18:10
12
$\begingroup$

Using the variance stabilizing Fisher's atan transformation, you can get the p-value as

pnorm( 0.5 * log( (1+r)/(1-r) ), mean = 0.5 * log( (1+0.75)/(1-0.75) ), sd = 1/sqrt(n-3) )

or whatever version of one-sided/two-sided p-value you are interested in. Obviuosly, you need the sample size n and the sample correlation coefficient r as inputs to this.

$\endgroup$
  • $\begingroup$ +1 Thanks for your Answer - It wasn't clear to me that the Fisher transform was appropriate or not in this case, but your answer helps clear that up. $\endgroup$ – Reinstate Monica - G. Simpson Aug 13 '11 at 17:04
  • $\begingroup$ @Gavin, you tried to clarify what the OP's intention was. I just assumed the modal situation in which a question like that would arise, and it looks as if it worked out :). $\endgroup$ – StasK Aug 13 '11 at 17:27
4
$\begingroup$

The distribution of r_hat around rho is given by this R function adapted from Matlab code at the webpage of Xu Cui. It's not that difficult to turn this into an estimate for the probability that an observed value "r" is improbable given a sample size of "n" and a hypothetical true value of "ro".

corrdist <- function (r, ro, n) {
        y = (n-2) * gamma(n-1) * (1-ro^2)^((n-1)/2) * (1-r^2)^((n-4)/2)
        y = y/ (sqrt(2*pi) * gamma(n-1/2) * (1-ro*r)^(n-3/2))
        y = y* (1+ 1/4*(ro*r+1)/(2*n-1) + 9/16*(ro*r+1)^2 / (2*n-1)/(2*n+1)) }

Then with that function you can plot the distribution of a null rho of 0.75, calculate the probability that r_hat will be less than 0.6 and shade in that area on the plot:

 plot(seq(-1,1,.01), corrdist( seq(-1,1,.01), 0.75, 10) ,type="l")
 integrate(corrdist, lower=-1, upper=0.6, ro=0.75, n=10)
# 0.1819533 with absolute error < 2e-09
 polygon(x=c(seq(-1,0.6, length=100), 0.6, 0), 
         y=c(sapply(seq(-1,0.6, length=100), 
         corrdist, ro=0.75, n=10), 0,0), col="grey")

enter image description here

$\endgroup$
4
$\begingroup$

Another approach that may be less exact than Fisher's tranformation, but I think could be more intuitive (and could give ideas about practical significance in addition to statistical significance) is the visual test:

 Buja, A., Cook, D. Hofmann, H., Lawrence, M. Lee, E.-K., Swayne,
 D.F and Wickham, H. (2009) Statistical Inference for exploratory
 data analysis and model diagnostics Phil. Trans. R. Soc. A 2009
 367, 4361-4383 doi: 10.1098/rsta.2009.0120

There is an implementation of this in the vis.test function in the TeachingDemos package for R. One possibly way to run it for your example is:

vt.scattercor <- function(x,y,r,...,orig=TRUE)
{
    require('MASS')
    par(mar=c(2.5,2.5,1,1)+0.1)
    if(orig) {
        plot(x,y, xlab="", ylab="", ...)
    } else {
        mu <- c(mean(x), mean(y))
        var <- var( cbind(x,y) )
        var[ rbind( 1:2, 2:1 ) ] <- r * sqrt(var[1,1]*var[2,2])
        tmp <- mvrnorm( length(x), mu, var )
        plot( tmp[,1], tmp[,2], xlab="", ylab="", ...)
    }
}

test1 <- mvrnorm(100, c(0,0), rbind( c(1,.75), c(.75,1) ) )
test2 <- mvrnorm(100, c(0,0), rbind( c(1,.5), c(.5,1) ) )

vis.test( test1[,1], test1[,2], r=0.75, FUN=vt.scattercor )
vis.test( test2[,1], test2[,2], r=0.75, FUN=vt.scattercor )

Of course if your real data is not normal or the relationship is not linear then that will be easily picked up with the above code. If you want to simultaniously test for those, then the above code would do that, or the above code could be adapted to better represent the nature of the data.

$\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.