是否有一个函数可以检验两个向量的相关性等于给定数字（例如 0.75）的假设？使用 cor.test 我可以测试 cor=0 并且可以查看 0.75 是否在置信区间内。但是有一个函数可以计算 cor=0.75 的 p 值吗？
x <- rnorm(10)
y <- x+rnorm(10)
cor.test(x, y)
如何使用 R 检验相关性等于给定值的假设？
机器算法验证
r
相关性
2022-03-13 15:22:21
3个回答
使用方差稳定费舍尔的 atan 变换，您可以得到 p 值
pnorm( 0.5 * log( (1+r)/(1-r) ), mean = 0.5 * log( (1+0.75)/(1-0.75) ), sd = 1/sqrt(n-3) )
或您感兴趣的任何版本的单侧/双侧 p 值。显然，您需要样本量n和样本相关系数r作为输入。
r_hat 在 rho 周围的分布由这个 R 函数给出，该函数改编自Xu Cui 网页上的Matlab 代码。在给定样本大小“n”和假设真值“ro”的情况下，将其转化为对观测值“r”不可能的概率的估计并不难。
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)) }
然后使用该函数，您可以绘制 0.75 的空 rho 分布，计算 r_hat 小于 0.6 的概率，并在绘图上的该区域着色：
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")
另一种可能不如费舍尔变换精确的方法，但我认为可能更直观（并且除了统计显着性之外，还可以提供有关实际意义的想法）是视觉测试：
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
vis.testR 包中的函数中有一个实现。TeachingDemos为您的示例运行它的一种可能方法是：
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 )
当然，如果你的真实数据不正常或者关系不是线性的，那么上面的代码很容易找到。如果你想同时测试这些，那么上面的代码会这样做，或者上面的代码可以适应更好地表示数据的性质。