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分数和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).

简单来说,计算的标准误差,z分数和p值glm只是标准误差,z分数和p值logit(s / (s + f)).这些是零假设的有效结果:logit(p) = 0换句话说,p = 0.5.因此,从获得的z得分和p值glm是测试是否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彼此显着不同.