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

10

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)

r correlation

— mozaika
źródło

2

To pytanie jest bardziej odpowiednie dla

— crossvalidated.com

1

@ sacha - najpierw sprawdź często zadawane pytania na stronie, często zadawane pytania dotyczące statystyki w witrynie stats.se zalecają pisanie na SO.

— Kev,

Pytanie „czy istnieje funkcja obliczająca wartość p dla cor = 0,75?” nie ma nic wspólnego z programowaniem. To pytanie statystyczne.

— Sacha Epskamp

Skonsuluję się ze statystykami i zobaczę, co myślą.

— Kev,

1

@mosaic Zarejestruj swoje konto tutaj. W ten sposób będziesz mógł powiązać swoje konto SO z obecnym.

— chl

12

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.

— StasK
źródło

+1 Dzięki za odpowiedź - Nie było dla mnie jasne, czy transformacja Fishera była odpowiednia, czy nie w tym przypadku, ale twoja odpowiedź pomaga to wyjaśnić.

— Gavin Simpson,

@Gavin, próbowałeś wyjaśnić, jaki był zamiar PO. Właśnie założyłem modalną sytuację, w której pojawi się takie pytanie, i wygląda na to, że się udało :).

— StasK,

4

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

wprowadź opis zdjęcia tutaj

— DWin
źródło

4

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.

— Greg Snow
źródło