-1

我想估计 13 个组的树冠覆盖的中位数以及 5、25、75 和 95 个百分位数（记录为 20 个中有或没有植被的点数）和总共 223 个样本。我之前发布了这个假设是 beta 分布，但这是不正确的。这是一份过期的手稿，这是遗失的最后一篇。如果有人可以帮助我完成（直到代码工作），我将不胜感激。我想我已经接近了，它只需要一些调整——我想。

（我进行了编辑以修复两个反对票，但我不确定有什么不清楚的地方）。

非常感谢！

下面是我的模型声明，R 代码和数据。照原样，我得到的错误是

model(model.file, data = data, inits = init.values, n.chains = n.chains,  : 
  RUNTIME ERROR:
Compilation error on line 16.
Subset out of range: re[14]

但请注意，我删除了下面的空格，错误是指可能性陈述。

model{
# priors
for (i in 1:13){
 alpha[i] ~ dunif(0, 1)
 re[i] ~ dnorm(0, 0.001)
}
#likelihood
for (i in 1:223) {
   canopy[i] ~ dbin(p[i], 20)
   logit(p[i]) <- alpha[site[i]] + re[i] 
}   
 median <- 1/(1+exp(alpha[site[i]]))
t4est1_100  <- step(median[1]-median[4])
t5est1_10 <- step(median[3]-median[4])
t6est10_100 <- step(median[2]-median[3])

}

R代码：

cover <- read.csv("f:\\brazil\\canopy2.csv", header=T)
library(R2jags)
library(rjags)
setwd("f://brazil")
site <- frag$site
canopy <- frag$canopy*20
N <- length(frag$site)

jags.data <- list("site", "canopy")
jags.params <- c("median", "test100MF","test100MT","test100fc","test100fa", 
"test100gv","test100hm","test100mc", "test100ca","test100ct", "test10MF",
"test10MT", "test10fc","test10fa", "test10gv", "test10hm", "test10mc", "test10ca", 
"test10ct", "test1MF", "test1MT", "test1fc",  "test1fa",  "test1gv", "test1hm", 
"test1mc", "test1ca", "test1ct", "t1est1_con","t2est10_con","t3est100_con",
"t4est1_100","t5est1_10","t6est10_100")
#inits1 <- list(a=0, sd=0)
#inits2 <- list(a=100, sd=50)
#jags.inits <- list(inits1, inits2)

jags.inits <- function() {
  list(alpha = 0, re=0)}

jagsfit2 <- jags(data=jags.data, inits=jags.inits, jags.params,
n.iter=1000000, n.burnin=20000, model.file="fragmodelbinom.txt")
my.coda <- as.mcmc(jagsfit2)
summary(my.coda, quantiles=c(0.05, 0.25,0.5,0.75, 0.95))
print(jagsfit2, digits=3)

数据：

structure(list(site = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 
6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 
7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 
8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 9L, 9L, 9L, 9L, 10L, 10L, 
10L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 12L, 12L, 12L, 12L, 
13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L
), canopy = c(0, 0.05, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 
0.6, 0.6, 0.65, 0.65, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0.75, 0.75, 
0.8, 0.8, 0.85, 0.85, 0.85, 0.85, 0.85, 0.85, 0.9, 0.9, 0.9, 
0.9, 0.95, 0.95, 0.95, 0.95, 0.95, 1, 1, 1, 1, 1, 1, 0.05, 0.2, 
0.25, 0.4, 0.4, 0.5, 0.6, 0.6, 0.65, 0.65, 0.75, 0.75, 0.75, 
0.8, 0.8, 0.8, 0.8, 0.85, 0.85, 0.85, 0.9, 0.9, 0.95, 0.95, 0.95, 
0.95, 1, 1, 1, 1, 1, 0, 0.2, 0.25, 0.3, 0.35, 0.4, 0.4, 0.45, 
0.45, 0.5, 0.5, 0.55, 0.6, 0.7, 0.7, 0.7, 0.8, 0.85, 0.9, 0.9, 
0.9, 0.95, 0.95, 0.95, 0.95, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 
0.1, 0.4, 0.4, 0.45, 0.5, 0.55, 0.55, 0.7, 0.7, 0.75, 0.8, 0.8, 
0.8, 0.9, 1, 1, 0.15, 0.2, 0.25, 0.25, 0.35, 0.5, 0.5, 0.55, 
0.65, 0.7, 0.7, 0.75, 0.8, 0.85, 0.85, 0.9, 0.9, 0.95, 0.95, 
1, 1, 1, 1, 0.05, 0.4, 0.6, 0.65, 0.65, 0.65, 0.7, 0.85, 0.95, 
1, 1, 1, 0.35, 0.4, 0.4, 0.5, 0.5, 0.55, 0.65, 0.65, 0.75, 0.75, 
0.8, 0.85, 0.9, 0.9, 1, 1, 1, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 
0.8, 0.8, 0.8, 0.85, 0.95, 0.95, 1, 1, 1, 1, 0.8, 0.85, 1, 1, 
1, 1, 1, 0, 0.05, 0.1, 0.15, 0.5, 0.6, 0.6, 0.75, 0.1, 0.35, 
0.6, 1, 0.4, 0.5, 0.55, 0.65, 0.65, 0.8, 0.9, 0.9, 0.9, 0.9, 
0.95, 0.95, 1)), .Names = c("site", "canopy"), class = "data.frame", 
row.names = c(NA, -227L))
4

2 回答 2

0

通过阅读您的代码，我可以看到两个明显的错误。在模型定义中，您有：

    for (i in 1:223) {
     canopy[i] ~ dbin(p[i], 20)
     logit(p[i]) <- alpha[site[i]] + re[i] 
    }

以上不可能是正确的，因为re参数只有 13 个组件，正如您定义的那样。我想应该是re[site[i]]。这是您遇到错误的原因。

第二个是行median <- 1/(1+exp(alpha[site[i]]))：它在任何 for 循环之外，而我想它应该在里面，因为它取决于i.

于 2014-09-21T20:35:59.070 回答
0

在 JAGS 中，dbin采用成功概率和试验次数，并返回一个表示成功次数的整数。在您的模型中，您指定的结果dbin是您的canopy值，它们都介于 0 和 1 之间 - 这是不正确的。

为什么要指定 20 次试验？您是否有 20 个点来评估该物种的存在/不存在？如果是这样，也许您应该将数据模型指定为：

y[i] ~ dbin(cover[i], 20)

其中y[i]是一个整数，给出了现场物种的存在数量i（这些是您的数据），并且cover是我们试图估计的东西，即记录为观察的物种的比例覆盖率i）。（Fukaya et al. (2010)提供了类似问题的代码，他们估计了调查地块上藤壶发生的概率——值得一读。）

这是一个示例，使用代表您的问题（或至少我对它的理解）的模拟数据，它将编译并且 JAGS 估计合理准确的覆盖值。

M <- function() {
  for (j in 1:13){
    alpha[j] ~ dnorm(0, 0.0001)
  }

  for (i in 1:length(y)) {
    y[i] ~ dbin(p[i], 20)
    logit(p[i]) <- alpha[site[i]] + eps[i]
    eps[i] ~ dnorm(0, sd^-2)
  }   
  sd ~ dunif(0, 100)
}

请注意，我使用了一个模糊的正常先验alpha，而不是您的 Uniform(0, 1) 先验，因为alpha它在 logit 标度上，因此不受 0 和 1 的限制。为了清楚起见，我已重命名re为eps并移动了它到观察水平回路。

site <- c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
          1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
          1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 
          2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 
          3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 
          3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 
          4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 
          5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 
          7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 8, 8, 8, 
          8, 8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 10, 10, 10, 11, 11, 11, 11, 
          11, 11, 11, 11, 12, 12, 12, 12, 13, 13, 13, 13, 13, 13, 13, 13, 
          13, 13, 13, 13, 13)

set.seed(1)
alpha <- rnorm(13) # 13 independent random site intercepts
sd <- 0.123        # a constant variance for obs around alphas
logit.p <- rnorm(length(site), alpha[site], sd) # logit probs
y <- sapply(plogis(logit.p), rbinom, n=1, size=20)

几个图表来显示我们模拟的数据：

plot(plogis(logit.p) ~ site, pch=20, las=1, ylab='p', ylim=c(0, 1))

在此处输入图像描述

plot(y ~ jitter(site), pch=20, las=1, xlab='site')

在此处输入图像描述

fit <- jags(list(y=y, site=site), NULL, c('sd', 'alpha', 'p'), M, 3, 1000)

将我们的真实值p与 JAGS 估计值进行比较：

plot(plogis(logit.p), fit$BUGSoutput$summary[
  grep('^p\\[', row.names(fit$BUGSoutput$summary)), '50%'], pch=20,
  xlab='True p', ylab='Estimated p (median and 95% CrI)', las=1, ylim=c(0, 1),
  panel.first=abline(0, 1, lwd=3, col='light gray'))
segments(plogis(logit.p), 
         fit$BUGSoutput$summary[
           grep('^p\\[', row.names(fit$BUGSoutput$summary)), '2.5%'], 
         y1=fit$BUGSoutput$summary[
           grep('^p\\[', row.names(fit$BUGSoutput$summary)), '97.5%'])

在此处输入图像描述

和阿尔法：

plot(alpha, fit$BUGSoutput$summary[
  grep('^alpha\\[', row.names(fit$BUGSoutput$summary)), '50%'], pch=20,
  xlab='True alpha', ylab='Estimated alpha (median +95% CrI)', las=1, 
  ylim=c(-1.5, 2.5), panel.first=abline(0, 1, lwd=3, col='light gray'))
segments(alpha, 
         fit$BUGSoutput$summary[
           grep('^alpha\\[', row.names(fit$BUGSoutput$summary)), '2.5%'], 
         y1=fit$BUGSoutput$summary[
           grep('^alpha\\[', row.names(fit$BUGSoutput$summary)), '97.5%'])

在此处输入图像描述

摘要（例如分位数）很容易从存储的样本中提取出来，其中一些已经出现在 JAGS 摘要中（即fit$BUGSoutput$summary）。在您的情况下，您可能对估计的分位数alpha（即平均站点覆盖率）感兴趣。访问这些的另一种方法，当您想要尚未计算的分位数时很有用R2jags，如下所示：

apply(fit$BUGSoutput$sims.list$alpha, 2, quantile, c(0.2, 0.5, 0.8))

#           [,1]       [,2]       [,3]     [,4]      [,5]       [,6]      [,7]
# 20% -0.6753085 0.08214699 -0.8994924 1.213293 0.2411191 -0.9441959 0.4038927
# 50% -0.6151128 0.15716925 -0.8299434 1.310542 0.3214097 -0.8205142 0.4909790
# 80% -0.5461859 0.23436643 -0.7632734 1.405498 0.3994087 -0.6982591 0.5887423
#          [,8]      [,9]       [,10]    [,11]     [,12]      [,13]
# 20% 0.6052447 0.6722630 -0.46050504 1.513086 0.2352861 -0.6516170
# 50% 0.6971180 0.8914973 -0.24458871 1.677613 0.4217723 -0.5443512
# 80% 0.7858756 1.1142837 -0.01018129 1.839783 0.6183198 -0.4309478

不过，您可能需要考虑一种分层方法，即alpha从属于共同分布的一组 alpha 中抽取。这将涉及以下内容：

M <- function() {
  for (j in 1:13){
    alpha[j] ~ dnorm(mean.alpha, sd.alpha^-2)
  }

  for (i in 1:length(y)) {
    y[i] ~ dbin(p[i], 20)
    logit(p[i]) <- alpha[site[i]] + eps[i]
    eps[i] ~ dnorm(0, sd^-2)
  }   
  sd ~ dunif(0, 100)
  mean.alpha ~ dnorm(0, 0.0001)
  sd.alpha ~ dunif(0, 100)
}

然后，这将告诉您有关组均值中站点之间的差异的信息（即通过检查sd.alpha）。您也可以考虑允许sd（即站点内观察值之间的差异）在各个站点之间变化，尽管您可能没有太多的数据可以这样做（并且sd无论如何拥有一个单一的常量可能是完全明智的......取决于系统） .

于 2014-09-22T10:30:33.437 回答