11
$\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$
2
  • $\begingroup$ why not just use the correlation coefficient and its SE? With 1.96*SE you can easily obtain the confidence interval for correlation coefficient! Or is this idea not correct? $\endgroup$
    – Tomas
    Aug 18, 2011 at 9:58
  • $\begingroup$ Doesn't look right since this can give a confidence interval that is not contained in (-1, 1). I guess this is why you use the Fisher transform. $\endgroup$
    – mosaic
    Aug 18, 2011 at 14:28

3 Answers 3

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$
2
  • $\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$ Aug 13, 2011 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, 2011 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$
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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