R：计算的随机分布的P值

Question

欲得到两个无规分布的观测x和y的P值，例如：

> set.seed(0) 
> x <- rnorm(1000, 3, 2) 
> y <- rnorm(2000, 4, 3) 

或：

> set.seed(0) 
> x <- rexp(50, 10) 
> y <- rexp(100, 11) 

假设T是我的测试统计量，定义为mean（x） - mean（y）= 0（这是H0），那么P值定义为：p-value = P [T> T_observed | H0成立]。
我试着这样做：

> z <- c(x,y) # if H0 holds then x and y are distributed with the same distribution 
> f <- function(x) ecdf(z) # this will get the distribution of z (x and y) 

然后计算p值我想这：

> T <- replicate(10000, mean(sample(z,1000,TRUE))-mean(sample(z,2000,TRUE))) # this is 
supposed to get the null distribution of mean(x) - mean(y) 
> f(quantile(T,0.05)) # calculating the p-value for a significance of 5% 

显然，这似乎并没有工作，我失去了什么？

Answer 1

您的意图非常好 - 通过自举采样（aka bootstrapping）来计算统计显着性。但是，平均值（样本（x，1000，TRUE）） - 平均值（样本（z，2000，TRUE））无法正常工作，因为这需要平均1000个z样本 - 平均2000个z样本。无论x和y的真实方式如何，这肯定会非常接近0。

我建议如下：x和y的

diff <- (sample(x, size = 2000, replace = TRUE) - sample(y, size = 2000, replace = TRUE)) 

2000样品（与替换）采取并计算差值。当然你也可以按照你的建议增加重复次数来增加信心。与pvalue相比，我更喜欢置信区间（confidence interval，CI），因为我认为它们更具信息性（与p值相比统计准确度相当）。使用平均值和标准误差如下顺然后可以计算：

stderror <- sd(diff)/sqrt(length(x)) 
upperCI <- mean(diff)+stderror 
lowerCI <- mean(diff)-stderror 
cat(lowerCI, upperCI) 

由于CI不包括0时，零假设被拒绝。请注意，结果将接近t检验（对于您的正常示例）CI结果R：

t <- t.test(x, y) 
cat(t$conf.int)