Jak sprawdzić hipotezę, że korelacja jest równa podanej wartości przy użyciu R?

Czy istnieje funkcja do testowania hipotezy, że korelacja dwóch wektorów jest równa danej liczbie, powiedzmy 0,75? Za pomocą testu cor.test mogę przetestować cor = 0 i mogę sprawdzić, czy 0,75 znajduje się w przedziale ufności. Ale czy istnieje funkcja obliczania wartości p dla cor = 0,75?

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

Używając wariancji stabilizującej transformację Atana Fishera , możesz uzyskać wartość p jako

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

lub jakiejkolwiek wersji jednostronnej / dwustronnej wartości p, którą jesteś zainteresowany. Oczywiście potrzebujesz wielkości próbki ni współczynnika korelacji próbki rjako danych wejściowych do tego.

Rozkład r_hat wokół rho jest podany przez tę funkcję R zaadaptowaną z kodu Matlab na stronie Xu Cui . Przekształcenie tego w oszacowanie prawdopodobieństwa, że ​​zaobserwowana wartość „r” jest nieprawdopodobna, biorąc pod uwagę wielkość próby „n” i hipotetyczną wartość rzeczywistą „ro”, nie jest trudne.

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

Następnie za pomocą tej funkcji możesz wykreślić rozkład zerowego rho na 0,75, obliczyć prawdopodobieństwo, że r_hat będzie mniejszy niż 0,6 i zacienić w tym obszarze na wykresie:

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

Innym podejściem, które może być mniej dokładne niż transformacja Fishera, ale myślę, że może być bardziej intuicyjne (i może dawać pomysły na znaczenie praktyczne oprócz znaczenia statystycznego) to test wizualny:

 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

Istnieje implementacja tego w vis.testfunkcji w TeachingDemospakiecie dla R. Jednym z możliwych sposobów uruchomienia go dla twojego przykładu jest:

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 )

Oczywiście, jeśli twoje rzeczywiste dane nie są normalne lub relacja nie jest liniowa, można to łatwo odczytać za pomocą powyższego kodu. Jeśli chcesz jednocześnie je przetestować, powyższy kod by to zrobił, lub powyższy kod można dostosować, aby lepiej reprezentował charakter danych.