How to test hypothesis that correlation is equal to given value using R?

Question

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)

Answer 1

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.

Answer 2

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")

Answer 3

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.