1

我想使用 JAGS 估计 13 个地点（9 个是鸟类，4 个是潜在栖息地）的冠层覆盖百分比的平均值和 sd。我正在使用 beta 分布来解释数据受 0 和 1 约束的事实。

我有适用于其他分布（泊松和对数正态）的模型语句的代码，我试图调整该代码，但我惨遭失败。

下面是 R 代码、模型语句和数据。我在 Windows Vista 中使用 R 3.1.1。如果您可以查看模型声明并将我理顺，我将非常感激。

谢谢，

杰夫

######## MODEL ##############
model{
  for (i in 1:227) {
    log(mean[i]) <- a[site[i]] 
    cover20p[i] ~ dbeta(1, 0.5)   
  }
  for (i in 1:13){
    a[i] ~ dnorm(0, tau) 
    median[i] <- exp(a[i])
  }
  sd ~ dunif(0, 10) 
  tau <- 1 / (sd*sd) # precision
} 

#########  R  code ########## 
frag <- read.csv("f:\\brazil\\TIandFRAG.csv", header=T)
library(R2jags)
library(rjags)
setwd("f://brazil")
site <- frag$site
cover20p <- frag$cover20p/100
N <- length(frag$site)

jags.data <- list("site", "cover20p")
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(a=c(0,0,0,0,0,0,0,0,0,0,0,0,0), sd=1)}

jagsfit <- jags(data=jags.data, inits=jags.inits, jags.params,
n.iter=1000000, n.burnin=20000, model.file="fragmodelbeta.txt")

my.coda <- as.mcmc(jagsfit)
summary(my.coda, quantiles=c(0.05, 0.25,0.5,0.75, 0.95))
print(jagsfit, digits=3)

##### DATA ###################    
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.95, 0.8, 0.85, 0.9, 0.35, 0.999, 0.999, 0.999, 
0.95, 0.55, 0.9, 0.85, 0.7, 0.65, 0.05, 0.6, 0.999, 0.999, 0.85, 
0.9, 1e-04, 0.45, 0.999, 0.7, 0.95, 0.5, 0.95, 0.6, 0.65, 0.7, 
0.4, 0.85, 0.6, 0.95, 0.75, 0.9, 0.85, 0.75, 0.7, 0.85, 0.3, 
0.7, 0.8, 0.7, 0.75, 0.8, 0.75, 0.95, 0.9, 0.05, 0.85, 0.6, 0.65, 
0.5, 0.85, 0.95, 0.85, 0.25, 0.75, 0.999, 0.65, 0.95, 0.8, 0.9, 
0.6, 0.8, 0.999, 0.2, 0.8, 0.4, 0.999, 0.95, 0.4, 0.999, 0.999, 
0.95, 0.45, 0.2, 0.7, 0.95, 0.7, 0.8, 0.5, 0.85, 0.55, 1e-04, 
0.25, 0.45, 0.999, 0.95, 0.999, 0.9, 0.6, 0.35, 0.95, 0.3, 0.999, 
0.999, 0.5, 0.4, 0.9, 0.999, 0.7, 0.999, 0.9, 0.999, 0.4, 0.55, 
0.8, 0.7, 0.999, 1e-04, 0.8, 1e-04, 0.7, 0.5, 0.8, 0.75, 1e-04, 
0.45, 0.1, 1e-04, 0.4, 0.55, 0.4, 0.999, 0.9, 0.9, 0.15, 0.55, 
0.35, 0.9, 0.65, 0.25, 0.999, 0.85, 0.999, 0.95, 0.7, 0.5, 0.7, 
0.2, 0.95, 0.999, 0.999, 0.25, 0.85, 0.5, 0.8, 0.75, 0.85, 0.7, 
0.95, 0.05, 0.65, 0.65, 0.999, 0.999, 0.999, 0.65, 0.4, 0.6, 
0.9, 0.85, 0.75, 0.5, 0.65, 0.999, 0.65, 0.55, 0.75, 0.4, 0.9, 
0.35, 0.999, 0.999, 0.4, 0.5, 0.8, 0.95, 0.95, 0.55, 0.7, 0.85, 
0.8, 0.8, 0.65, 0.999, 0.6, 0.5, 0.999, 0.8, 0.999, 0.45, 0.999, 
0.999, 0.8, 0.85, 0.999, 0.999, 0.999, 0.999, 0.5, 0.6, 0.15, 
0.75, 0.6, 0.1, 0.05, 1e-04, 0.999, 0.6, 0.1, 0.35, 0.9, 0.9, 
0.95, 0.95, 0.9, 0.55, 0.65, 0.9, 0.4, 0.999, 0.65, 0.5, 0.8)), .Names = c("site", 
"canopy"), class = "data.frame", row.names = c(NA, -227L))
4

2 回答 2

0

我认为您可以使用 logit 模型来计算概率。也许像下面这样。

首先，我将您的冠层观测值转换回我怀疑它们开始的格式，即每个站点的 20 个样本中冠层命中的数量。我将 0.0001 设置为 0，将 0.999 设置为 1，然后将其他canopy值乘以 20。

d$hits <- ifelse(d$canopy < 0.05, 0, ifelse(d$canopy > 0.95, 20, d$canopy * 20))

M <- function() {
  for (i in 1:n) {
    hits[i] ~ dbin(p[site[i]], 20)
  }
  for (j in 1:nsites) {
    logit.p[j] ~ dnorm(mu, sigma^-2)
    logit(p[j]) <- logit.p[j]
  }
  mu ~ dnorm(0, 0.0001) # uninformative prior for grand mean of logit(p)
  sigma ~ dunif(0, 10) # uninformative prior for sd of logit(p)
}

j <- jags(list(site=d$site, hits=d$hits, n=nrow(d), nsites=length(unique(d$site))), 
          NULL, 'p', M)

plot(j$BUGSoutput$summary[-1, '50%'], pch=20, xlab='site', xaxt='n', las=1, 
     ylim=c(0, 1), ylab=expression("p (median" %+-% "95% credible interval)"))
segments(1:13, j$BUGSoutput$summary[-1, '2.5%'], 
         y1=j$BUGSoutput$summary[-1, '97.5%'])
axis(1, 1:13, 1:13)

在此处输入图像描述

于 2014-10-29T12:45:35.697 回答
0

在您的模型中，您将 cover20p 作为变量之一，但在 frag data.frame 中有树冠数据。我怀疑你想要canopy[i] ~ dbeta(1,0.5)在你的模型规范和canopy <- frag$canopy你jags.params = "median"的 r 代码中。

于 2014-08-29T21:28:07.657 回答