r - R中的成功/失败错误估计

Question

我有成功/失败数据(在一定时期内存活/死亡的树木)，并想估计与我的每个观察(7 个地点)相关联的二项式分布的误差。到目前为止我一直在使用 glm这样做:

s <- c(1,20,0,40,2,1,0) # success
f <- c(2,0,20,4,50,0,1) # failure

#for each observation I would calculate this error: 

error <- vector ()  
z_scores <- vector ()  
p_value <- vector ()  

  for (i in 1:7) {
    models <- glm (cbind (s[i], f[i]) ~ 1, family = 'binomial')
    error [i] <- summary (models)$coefficients[2]
    z_scores [i] <- summary (models)$coefficients[3]
    p_value [i] <- summary (models)$coefficients[4]
  }

这会是最好的方法吗？

这里如何估计二项式分布的概率？

请注意，无论成功和失败的次数如何，当 s 时，我的错误都非常高。或 f是 =0

Answer 1

这里有一些代码可以在不使用 glm 的情况下重新计算大部分结果(除了由零引起的极端结果) ，我解释了它们背后的含义。

s <- c(1, 20, 0, 40, 2, 1, 0) # success
f <- c(2, 0, 20, 4, 50, 0, 1) # failure

#for each observation I would calculate this error: 

error <- vector()  
z_scores <- vector()  
p_value <- vector()  

for (i in 1:7) {
    models <- glm(cbind(s[i], f[i]) ~ 1, family = 'binomial')
    error[i] <- summary(models)$coefficients[2]
    z_scores[i] <- summary(models)$coefficients[3]
    p_value[i] <- summary(models)$coefficients[4]
}

logit <- function(x){
    log(x / (1 - x))
}

dlogit <- function(x){
    1 / x / (1 - x)
}

p_hat <- s / (s + f)
## sqrt(p_hat * (1 - p_hat) / (s + f))
## is the standard error of p_hat
## error1 is the standard error of logit(p_hat)
error1 <- dlogit(p_hat) * sqrt(p_hat * (1 - p_hat) / (s + f))
## divide the estimation by the standard error, you get z-score
z_scores1 <- logit(p_hat) / error1
p_value1 <- 2 * pnorm(-abs(z_scores1))

您需要知道的第一件事是标准误差、z 得分、p 值等背后的基本原理。在统计中，我们首先有一些模型(在这种情况下，二项式模型:s|(s+f) ~ Binomial(s + f, p))，我们想使用它以适应我们拥有的数据和

1) 获得估计值(在本例中为 p)

2) 由于数据是随机生成的，我们想知道我们的估计有多好，这里有标准误差、z 分数和 p 值来“衡量估计中的随机性”，这里有一些重要的“技巧”:因为我们不知道产生数据的真正机制，我们只能通过假设来近似计算我们估计中的随机性

a) 我们的模型是(或类似于)数据生成的真实机制，并且

b) 实际参数与我们的估计相似(这通常需要大样本量，在这种情况下，样本量只是 s + f ，因此 s + f 必须足够大才能进行推断(标准误差、z-score 和p 值)验证)。我们可以看到，在 i = 1、6 和 7 的情况下，样本量非常小，这使得相应的标准误、z 分数和 p 值令人难以置信。

然后我可以谈谈我的计算背后的技术细节以及它们的含义。在 glm , 除了 Binomial(n, p)模型，您还假设 p 的模型像这样:
logit(p) ~ N(mu, sigma^2)
而 logit 函数就像我代码中的那样。

在这个简单的例子中，二项式概率的估计 p只是 p_hat <- s / (s + f) (不管是否使用glm)，从二项式变量的方差公式，我们可以得到估计概率p的方差。是 p * (1 - p) / n , 在这里，如果我们认为 p_hat <- s / (s + f)与真实的相似p由假设b，并用它代替p ，我们可以得到估计值 p 的标准误差.遵循CLT和Delta方法，当样本量足够大时，我们可以处理s / (s + f)或 logit(s / (s + f))遵循正态分布，例如，s / (s + f)大约是 N(p, s * f / (s + f) ^ 3)和 logit(s / (s + f))大约是 N(logit(p), dlogit(s / (s + f)) ^ 2 * s * f / (s + f) ^ 3) .

简单地说，glm 的标准误、z 分数和 p 值计算只是 logit(s / (s + f)) 的标准误差、z 分数和 p 值.这些是零假设的有效结果:logit(p) = 0 ，换句话说，p = 0.5 .因此，从 glm 获得的 z 分数和 p 值是测试是否s和 f当样本大小 s + f 时发生的概率相等很大。

然后我再说说0引起的极值。当s或 f等于 0，f 的估计概率或 s发生将是1，如果这是真的，数据生成机制实际上是非随机的!!一开始我说过我们使用我们的估计来近似计算我们估计中的随机性，并且在 s 的情况下或 f等于 0，如果我们使用我们的估计作为基本事实，我们应该 100% 相信我们的估计，这有点荒谬。在这种情况下，很多方法如 glm将无效。一般来说，如果样本量s + f足够大，我们认为s的概率或 f如果 s = 0，发生的情况真的很小或 f = 0 ，但如果样本量真的很小，如案例 6 或案例 7，我们实际上无法得出任何结论。

总之，如果二项式模型为真，则来自 glm结果，我的代码和上面提供的分析，我们可以说在 i = 2, 3, 4, 5 的情况下，s的概率和 f彼此显着不同。