Success/failure R 中的误差估计

Question

我有 success/failure 数据（在特定时期内 survived/died 的树），我想估计二项分布的误差与我的每个观察结果（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-score、p-value等背后的基本原理。在统计中，我们首先有一些模型（在这种情况下，二项式模型：s|(s+f) ~ Binomial(s + f, p)) 我们想用它来拟合我们拥有的数据

1) 获取估计值（在本例中为 p）

2) 由于数据是随机生成的，我们想知道我们的估计有多好，这里是标准误差、z 分数和 p 值 "measure the randomness in the estimation"，这里是一些重要的 "trick": 由于我们不知道产生数据的真实机制，我们只能通过假设来近似计算我们估计中的随机性

a) 我们的模型是（或类似于）真正的数据生成机制

b) 真实参数与我们的估计相似（这通常需要较大的样本量，在这种情况下，样本量刚好s + f，因此s + f必须足够大才能使推理（标准误差、z 分数和 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-scores 和 p-values 是为了检验 s 和 f 在样本量 s + f 大的情况下是否等概率发生。

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

综上所述，如果二项式模型成立，根据glm结果，我的代码和我上面提供的分析，我们可以说在i = 2, 3, 4, 5的情况下，[=的概率40=] 和 f 彼此明显不同。

Success/failure R 中的误差估计

Success/failure error estimation in R

binary

r

distribution

confidence-interval

glm