伯努利数据模型的 Rao 分数测试的 R 代码是否正确？

Question

我是一个完全的统计菜鸟，也是 R 的新手，因此才有这个问题。当一个人的 data 是二进制的并且每个观察值都具有伯努利分布时，我试图找到 Rao score 的特定情况的实现。我偶然发现了 R 语言中的 anova，但未能理解如何使用它。因此，我尝试自己为这个特殊案例实施 Rao 评分：

rao.score.bern <- function(data, p0) {
  # assume `data` is a list of 0s and 1s
  y <- sum(data)
  n <- length(data)
  phat <- y / n

  z <- (phat - p0) / sqrt(p0 * (1 - p0) / n)
  p.value <- 2 * (1 - pnorm(abs(z)))
}

我很确定我的代码中存在错误，因为它在以下情况下只产生两个不同的 p 值：

p0 <- 1 / 4
p <- seq(from=0.01, to=0.5, by=0.01)
n <- seq(from=5, to=70, by=1)
g <- expand.grid(n, p)

data <- apply(g, 1, function(x) rbinom(x[1], 1, x[2]))
p.values <- sapply(data, function(x) rao.score.bern(x[[1]], p0))

谁能告诉我问题出在哪里？你能告诉我 R 中的内置解决方案吗？

Answer 1

先测试，再调试。

测试

rao.score.bern 有用吗？

rao.score.bern(c(0,0,0,1,1,1), 1/6))

这个returns...没什么！通过将最终行替换为

来修复它

2 * (1 - pnorm(abs(z)))

这消除了不必要的赋值。

rao.score.bern(c(0,0,0,1,1,1), 1/6))

[1] 0.02845974

好的，现在我们开始了。

调试

不幸的是，代码仍然无效。让我们通过取消对 rao.score.bern 的调用并将其替换为向我们显示输入的内容来进行调试。不要将它应用于您创建的大输入！使用一小块：

sapply(data[1:5], function(x) x[[1]])

[1] 0 0 0 0 0

这不是您所期望的，是吗？对于 data 的每个元素，它只返回一个零。这个呢？

sapply(data[1:5], function(x) x)

[[1]]
[1] 0 0 0 0 0
[[2]]
[1] 0 0 0 0 0 0 
...
[[5]]
[1] 0 0 0 0 0 0 0 0 0

好多了！ sapply 调用中的变量 x 指的是整个向量，这是您要传递给例程的内容。从何而来

p.values <- sapply(data, function(x) rao.score.bern(x, p0)); hist(p.values)

伯努利数据模型的 Rao 分数测试的 R 代码是否正确？

Is this R code of Rao score test for the Bernoulli data model correct?

r

statistics

bernoulli-probability

测试

调试