首页 文章

使用Bebber等人2007查找估计的总物种和该数字的置信区间[已结束]

提问于
浏览 次
1

我的问题是统计数据和R相关的混合 . 我正在尝试实施Bebber等人的模型 . 2007（使用发现曲线预测未知物种数量）估算我的数据集中未发现物种的总数（通过找到预期的物种总数），以及围绕物种总数预测的置信区间 .

物种总数的估计由E（St）= k（Ntot-Nt-1）给出，其通过拟合Nt-1上的St的glm（具有拟线性误差）来计算，其中St是物种数每年，Nt-1是前一年的累积物种数，k是发现率 . 从glm输出，我使用截距和斜率来计算估计的Ntot：Ntot = -intercept / slope我认为我正在理解如何做到这一点 . 这是一个示例df和我到目前为止使用的代码：

St <- c(3, 1, 6, 6, 2, 4, 2, 2, 3, 2, 4)
Nprev <- c(43, 46, 47, 53, 59, 61, 65, 67, 69, 72, 74)
samp_data <- data.frame(St, Nprev)
mymodel <- glm(samp_data$St ~ samp_data$Nprev, family = quasipoisson())
summary(mymodel)
Ntot <- -1.75918/-0.01019

Here is where I start having trouble understanding what to do. Bebber等人给出了我不完全理解的Ntot置信限制的说明 . 他们说首先定义一个新变量R = Ntot M - Nt-1 . 接下来是使用Poisson误差和没有截距的身份链接为St对抗R的glm .

The first thing I am confused about is where M comes from... 这是第一个glm的色散参数吗？然后他们继续说，"This forces the regression through Ntot + M. The model will give the same residual deviance as the best fit when M is zero, and a larger residual deviance as M diverges from zero. Changes in deviance scaled by the dispersion have an F distribution, allowing calculation of confidence limits."

I am not understanding how to get confidence limits for Ntot from this glm. 我是广义线性模型的新手，所以也许这就是我的问题 . 我希望这里的某些人可以为我阐明这一点 . 可以在此处找到该论文的副本（相关部分在第1652页"The model."下）：Bebber paper link

许多人提前感谢您提供的任何见解！

r glm confidence-interval

1 回答

  • 0

    以下是Ntot的自举置信区间的R解决方案：

    # exampel from question
St <- c(3, 1, 6, 6, 2, 4, 2, 2, 3, 2, 4)
Nprev <- c(43, 46, 47, 53, 59, 61, 65, 67, 69, 72, 74)
samp_data <- data.frame(St, Nprev)
mymodel <- glm(samp_data$St ~ samp_data$Nprev, family = quasipoisson())
summary(mymodel)
mymodel %>% str()
mymodel$coefficients
Ntot <- -1.75918/-0.01019

# bootstrap version
library(boot)
# put the process into a function, to be used in boot()
calcNtot <- function(data, inds, lhs, rhs){
  f <- as.formula(paste(lhs, " ~ ", rhs))
  # lhs/rhs means left/right hand side of the formula
  model <- glm(f, family = quasipoisson(), data = data, subset = inds)
  Ntot <- -model$coefficients[1]/model$coefficients[2]
  unname(Ntot)
}

calcNtot(samp_data, 1:nrow(samp_data), "St", "Nprev")
# it gives the same Ntot as the example

Ntot_boot <- boot(data = samp_data, calcNtot, 999, lhs = "St", rhs = "Nprev")
boot.ci(Ntot_boot, type = "basic")
    回复于

相关问题