베타 배포판을 사용하는 그룹 평균을 산정하기위한 JAGS 코드

JAGS를 사용하여 13 개 장소 (새들은 9 개, 잠재적 서식지는 4 개)에 대한 캐노피 비율의 평균과 SD를 추정하고 싶습니다. 데이터가 0과 1로 묶여 있다는 사실을 설명하기 위해 베타 배포판을 사용하고 있습니다.베타 배포판을 사용하는 그룹 평균을 산정하기위한 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

분명한 것은 R/JAGS가 모델을 초기화 한 다음 모델을 실행하는 대신 결과를 제공한다는 것 (그리고 한 시간 정도 소요)이었습니다. 그렇다면 나는 결과가 0과 1에 의해 구속 될 것을 기대했지만 그들은 벗어났다. 나는 모델 진술이 정확하지 않다는 것을 안다. 도와 주셔서 감사합니다. –

, 당신은 변수의 하나로서 cover20p이 있지만 파편 data.frame의 캐노피에 대한 데이터를 가지고있다. 모델 사양에 canopy[i] ~ dbeta(1,0.5)을 입력하고 r 코드에 canopy <- frag$canopy 및 jags.params = "median"을 입력하고 싶습니다.

출처

2014-08-29 21:28:07

감사합니다. 변수 문제는 새 데이터 프레임을 만들고 변수의 이름을 변경하지 않았기 때문에 내 감독이었습니다 (원래 데이터 프레임에는 28 개의 열이 있음). 모델 문이 완전히 잘못되었다고 가정하십시오. 내가 옳다고 생각하는 유일한 사실은 내가 우도에 베타 버전을 사용한다는 사실입니다. –

나는 당신이 당신의 가능성을 위해 로짓 모델을 사용할 수 있다고 생각한다. 어쩌면 다음과 같은 것일 수 있습니다.

먼저 캐노피 관측을 처음 생각한 형식, 즉 각 사이트에서 20 개의 샘플 중에서 캐노피 조회수로 변환합니다. 나는 1 0 0.999 0.0001를 설정하고, 정확히이 코드로 잘못 무엇 (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

베타 배포판을 사용하는 그룹 평균을 산정하기위한 JAGS 코드

답변

관련 문제