¿Cómo probar la hipótesis de que la correlación es igual al valor dado usando R?

10

¿Existe una función para probar la hipótesis de que la correlación de dos vectores es igual a un número dado, digamos 0.75? Usando cor.test puedo probar cor = 0 y puedo ver si 0.75 está dentro del intervalo de confianza. Pero, ¿hay una función para calcular el valor p para cor = 0.75?

x <- rnorm(10)
y <- x+rnorm(10)
cor.test(x, y)
mosaico
fuente
2
Esta pregunta es más adecuada para crossvalidated.com
Sacha Epskamp
1
@sacha: primero consulte las preguntas frecuentes de un sitio, las preguntas frecuentes del sitio stats.se recomiendan que las preguntas de programación que usan R se publiquen en SO.
Kev
La pregunta "¿hay una función para calcular el valor p para cor = 0.75?" No tiene nada que ver con la programación. Es una pregunta estadística.
Sacha Epskamp
Consultaré a la gente de estadísticas y veré qué piensan.
Kev
1
@mosaic Por favor, registre su cuenta aquí. De esta manera, podrá asociar su cuenta SO con la actual.
chl

Respuestas:

12

Usando la variación que estabiliza la transformación atan de Fisher , puede obtener el valor p como

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

o cualquier versión del valor p unilateral / bilateral que le interese. Obviamente, necesita el tamaño de nla muestra y el coeficiente de correlación de la muestra rcomo entradas para esto.

StasK
fuente
+1 Gracias por su respuesta: no estaba claro para mí que la transformación de Fisher fuera apropiada o no en este caso, pero su respuesta ayuda a aclarar eso.
Gavin Simpson
@Gavin, trataste de aclarar cuál era la intención del OP. Simplemente asumí la situación modal en la que surgiría una pregunta como esa, y parece que funcionó :).
StasK
4

La distribución de r_hat alrededor de rho está dada por esta función R adaptada del código de Matlab en la página web de Xu Cui . No es tan difícil convertir esto en una estimación de la probabilidad de que un valor observado "r" sea improbable dado un tamaño de muestra de "n" y un valor hipotético verdadero de "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)) }

Luego, con esa función, puede trazar la distribución de un nulo rho de 0.75, calcular la probabilidad de que r_ sea menor que 0.6 y sombrear esa área en el gráfico:

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

ingrese la descripción de la imagen aquí

DWin
fuente
4

Otro enfoque que puede ser menos exacto que la transformación de Fisher, pero creo que podría ser más intuitivo (y podría dar ideas sobre la significación práctica además de la significación estadística) es la prueba visual:

 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

Hay una implementación de esto en la vis.testfunción en el TeachingDemospaquete para R. Una posible forma de ejecutarlo para su ejemplo es:

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 )

Por supuesto, si sus datos reales no son normales o la relación no es lineal, eso se captará fácilmente con el código anterior. Si desea realizar una prueba simultánea de esos, entonces el código anterior haría eso, o el código anterior podría adaptarse para representar mejor la naturaleza de los datos.

Greg Snow
fuente