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

Question

我是一个完整的统计菜鸟和R的新手,因此这个问题.我试图找到特定情况的Rao分数的实现,当一个data是二元的并且每个观察都有bernoulli分布.我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)); HIST(p.values)