首页 文章

R:使用模拟计算p值

提问于
浏览
0

我编写了这段代码来对两个随机分布的观察x和y运行一个测试统计

mean.test <- function(x, y, B=10000,
alternative=c("two.sided","less","greater"))
{
p.value <- 0
alternative <- match.arg(alternative)
s <- replicate(B, (mean(sample(c(x,y), B, replace=TRUE))-mean(sample(c(x,y), B, replace=TRUE))))
t <- mean(x) - mean(y) 
p.value <- 2*(1- pnorm(abs(quantile(T,0.01)), mean = 0, sd = 1, lower.tail = 
TRUE, log.p = FALSE))   #try to calculate p value 
data.name <- deparse(substitute(c(x,y)))
names(t) <- "difference in means"
zero <- 0
names(zero) <- "difference in means"
return(structure(list(statistic = t, p.value = p.value,
method = "mean test", data.name = data.name,
observed = c(x,y), alternative = alternative,
null.value = zero),
class = "htest"))
}

代码使用蒙特卡罗模拟来生成检验统计量均值(x) - 均值(y)的分布函数,然后计算p值,但显然我错过了定义这个p值,因为:

> set.seed(0)
> mean.test(rnorm(1000,3,2),rnorm(2000,4,3))

输出应该如下:

mean test
data: c(rnorm(1000, 3, 2), rnorm(2000, 4, 3))
difference in means = -1.0967, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0

但我得到了这个:

mean test
data:  c(rnorm(1000, 3, 2), rnorm(2000, 4, 3))
difference in means = -1.0967, p-value = 0.8087
alternative hypothesis: true difference in means is not equal to 0

有人可以向我解释这个错误吗?

1 回答

  • 2

    据我所知,你的代码中有很多错误和错误:

    • quantile(T, 0.01) - 这里 T == TRUE ,所以你要计算1的分位数 .

    • 永远不会使用对象 s .

    • mean(sample(c(x,y), B, replace=TRUE)) 你想在这做什么? c() 函数结合了 xy . 因为你不知道他们来自哪个人口,因此抽样是没有意义的

    • 当您计算检验统计量 t 时,它应该取决于方差(和样本量) .

相关问题