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

Question

我是一个完整的统计菜鸟，并且是 R 的新手，因此提出了这个问题。当一个人是二元的并且每个观察都具有伯努利分布时，我试图为特定情况找到Rao 分数的实现。data我偶然发现anova了 R 语言，但不明白如何使用它。因此，我自己尝试为这个特殊情况实施 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))

这返回......什么都没有！通过将最终行替换为来修复它

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

这消除了不必要的分配。

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

[1] 0.02845974

好的，现在我们到了某个地方。

调试

不幸的是，代码仍然不起作用。让我们通过拉出调用rao.score.bern并将其替换为向我们显示输入的内容来进行调试。不要将其应用于您创建的大型输入！使用它的一小部分：

sapply（数据[1:5]，函数（x）x[[1]]）

[1] 0 0 0 0 0

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

sapply（数据[1:5]，函数（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

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

p.values <- sapply(data, function(x) rao.score.bern(x, p0)); 历史（p.values）