用於估算具有Beta分佈的羣體平均值的JAGS代碼

我想使用JAGS估算13個地點（9個是鳥類，4個是潛在棲息地）的樹冠覆蓋百分比的均值和sd。我正在使用測試版分佈來說明數據被0和1綁定的事實。用於估算具有Beta分佈的羣體平均值的JAGS代碼

我有模型語句的代碼，可以完美地適用於其他分佈（泊松和對數正態），我正在嘗試以適應代碼，但我失敗了。

下面是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))

來源

2014-08-29 Jeff Stratford

乘其他canopy值？它是拋出一個錯誤（如果是這樣，錯誤是什麼？）或給你你沒有想到的結果（如果是這樣，你期望什麼以及結果有什麼不同？）？ – tkmckenzie 2014-08-29 19:36:53

顯而易見的是R/JAGS初始化模型，然後給出結果而不是運行模型（並花費一個小時左右的時間來完成）。然後我期望結果受到0和1的約束，但他們關閉了。我知道模型聲明是不正確的。謝謝您的幫助。 – 2014-08-31 01:15:03

在你的模型，你有cover20p作爲變量之一，但有數據雨棚在FRAG data.frame。我懷疑你的代碼中需要canopy[i] ~ dbeta(1,0.5)，canopy <- frag$canopy和jags.params = "median"。

來源

2014-08-29 21:28:07

謝謝吉姆。自從我創建了一個新的數據框並沒有更改變量的名稱（原始數據框有28列）後，變量問題就是我的疏忽。假設模型聲明是完全錯誤的。我認爲唯一正確的是我使用beta分佈來評估可能性。 – 2014-08-31 01:14:29

我想你可以用你的概率的logit模型。可能類似於以下內容。

首先，我將您的冠層觀測數據轉換回我懷疑他們開始的格式，即每個地點20個樣本中的冠層數量。我設置0.0001〜0和0.999至1，和究竟出了毛病，這個代碼由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)

enter image description here

來源

2014-10-29 12:45:35 jbaums

用於估算具有Beta分佈的羣體平均值的JAGS代碼

回答

相關問題