标签: jags

问题

我正在尝试在JAGS中建立一个生存模型,允许时变协变量.我希望它是一个参数模型 - 例如,假设生存遵循威布尔分布(但我想让危险变化,所以指数太简单了).因此,这基本上是可以在flexsurv包中完成的贝叶斯版本,它允许参数模型中的时变协变量.

因此,我希望能够以"计数过程"形式输入数据,其中每个主题有多行,每行对应于其协变量保持不变的时间间隔(如本pdf或此处所述).这是包或包允许的(start, stop]配方.survivalflexurv

不幸的是,关于如何在JAGS中进行生存分析的每一个解释似乎都假设每个主题一行.

我试图采用这种更简单的方法并将其扩展到计数过程格式,但模型没有正确估计分布.

失败的尝试:

这是一个例子.首先,我们生成一些数据:

library('dplyr')
library('survival')

## Make the Data: -----
set.seed(3)
n_sub <- 1000
current_date <- 365*2

true_shape <- 2
true_scale <- 365

dat <- data_frame(person = 1:n_sub,
                  true_duration = rweibull(n = n_sub, shape = true_shape, scale = true_scale),
                  person_start_time = runif(n_sub, min= 0, max= true_scale*2),
                  person_censored = (person_start_time + true_duration) > current_date,
                  person_duration = ifelse(person_censored, current_date - person_start_time, true_duration)
)

  person person_start_time …

我正在使用rjagsR库.该函数coda.samples产生一个mcmc.list,例如(from example(coda.samples)):

library(rjags)
data(LINE)
LINE$recompile()
LINE.out <- coda.samples(LINE, c("alpha","beta","sigma"), n.iter=1000)
class(LINE.out)
[1] "mcmc.list"

但是,我想使用该plot.bugs函数,它需要一个bugs对象作为输入.

是否可以将对象从对象转换mcmc.list为bugs对象,以便plot.bugs(LINE.out)？

请注意,stats.SE上有一个类似的问题,一个多月没有得到答复.这个问题的结果是在2012年8月29日结束.

更多提示:

我发现R2WinBUGS包有一个函数"as.bugs.array"函数 - 但是不清楚该函数如何应用于mcmc.list.

我试图在JAGS中使用逻辑回归模型,但我有(#success y,#attempts n)形式的数据,而不是二进制变量.在R中,可以通过使用glm(y/n~)和"weights"参数将模型拟合到这些数据,但我不确定如何在JAGS中使用它.

这是一个简单的例子,我希望解决我想要问的问题.请注意,我使用的是rjags包.谢谢你的帮助!

y <- rbinom(10, 500, 0.2)
n <- sample(500:600, 10)
p <- y/n
x <- sample(0:100, 10) # some covariate

data <- data.frame(y, n, p, x)

model <- "model{
# Specify likelihood
for(i in 1:10){
    y[i] ~ dbin(p[i], n[i])
    logit(p[i]) <- b0 + b1*x
}

# Specify priors
b0 ~ dnorm(0, 0.0001)
b1 ~ dnorm(0, 0.0001)
}"

我正在尝试在R中运行rjags(通过Rstudio)来估计模型的参数alpha&beta和超参数tau.nu:

y_i|x_i~pois(eta_i),
eta_i=exp(alpha + beta*x_i + nu_i),
nu_i~N(0,tau.nu)

有我的代码:

#generating data
N = 1000
x = rnorm(N, mean=3,sd=1) 
nu = rnorm(N,0,0.01)
eta = exp(1 + 2*x + nu)
y = rpois(N,eta) 
data=data.frame(y=y,x=x)
###MCMC
library(rjags)
library(coda)
mod_string= "model {  
  for(i in 1:1000) {
    y[i]~dpois(eta[i])
    eta[i]=exp(alpha+beta*x[i]+nu[i])
    nu[i]~dnorm(0,tau.nu)
  }
  alpha  ~ dnorm(0,0.001)
  beta  ~ dnorm(0,0.001) 
  tau.nu ~ dgamma(0.01,0.01) 
}"

params = c("alpha","beta","tau.nu")

inits = function() {
  inits = list("alpha"=rnorm(1,0,100),"beta"=rnorm(1,0,80),"tau.nu"=rgamma(1,1,1))
}
mod = jags.model(textConnection(mod_string), data=data, inits=inits, n.chains =3)
update(mod,5000)
mod_sim = coda.samples(model=mod, …

我正在尝试编写一个Winbugs/Jags模型来建模多粒度主题模型(完全是本文 - > http://www.ryanmcd.com/papers/mg_lda.pdf)

在这里,我想根据特定值选择不同的分布.对于Eg:我想做点什么

`if ( X[i] > 0.5 )
{
Z[i] ~ dcat(theta-gl[D[i], 1:K-gl])
W[i] ~ dcat(phi-gl[z[i], 1:V])
}
else 
{
Z[i] ~ dcat(theta-loc[D[i], 1:K-loc])
W[i] ~ dcat(phi-loc[z[i], 1:V])
}
`

这可以在Winbugs/JAGS中完成吗？

我第一次尝试用 rjags 拟合分层模型；它在纸上看起来很简单，但我收到一个“错误解析模型文件：第 5 行靠近“[”的语法错误”，我完全无法解释。

你能帮助我并告诉我我做错了什么吗？

data = list('P.hat'=c(0.0032, 0.0045, 0.077), 'R'=c(34580, 37932, 46724), 'N'=c(10028321, 15674923, 21426662), 's.over.rootn'=c(0.02, 0.006, 0.017), 'n'=1, 'tmax'=3)

cat('model{
## likelihoods ##
for(i in 1:n){
    for(j in 1:tmax){
        P.hat[i,j]  ~  dnorm(pi[j], (1/pow(s.over.rootn,2))[j])   
        R[i,j]      ~  dbin(theta[j], N[j])
}}
## daterministic relations ##
        gam         <- m*vs+(1-m)*va
for(j in 1:tmax){   
        theta[j]    <- (pi[j]*beta*gam)/(gam*dt+(1-gam)*du)
}
## priors ##
for(j in 1:tmax){
        pi[j]       ~  dbeta(1, 1)
}
        beta        ~  dbeta(1, 1)
        m           ~  dbeta(1, 1)                   
        vs          ~  dbeta(1, 1) …

我正在使用 RJAGS 修改现有模型。我想并行运行链，偶尔检查 Gelman-Rubin 收敛诊断，看看我是否需要继续运行。问题是，如果我需要根据诊断值恢复运行，重新编译的链会从第一个初始化的先前值而不是链停止的参数空间中的位置重新启动。如果我不重新编译模型，RJAGS 就会抱怨。有没有办法在链条停止时存储它们的位置，以便我可以从停止的地方重新初始化？这里我将举一个非常简单的例子。

示例1.bug：

model {
  for (i in 1:N) {
      x[i] ~ dnorm(mu,tau)
  }
  mu ~ dnorm(0,0.0001)
  tau <- pow(sigma,-2)
  sigma ~ dunif(0,100)
}

parallel_test.R:

#Make some fake data
N <- 1000
x <- rnorm(N,0,5)
write.table(x,
        file='example1.data',
        row.names=FALSE,
        col.names=FALSE)

library('rjags')
library('doParallel')
library('random')

nchains <- 4
c1 <- makeCluster(nchains)
registerDoParallel(c1)

jags=list()
for (i in 1:getDoParWorkers()){
  jags[[i]] <- jags.model('example1.bug',
                          data=list('x'=x,'N'=N))
}

# Function to combine multiple mcmc lists into a single one
mcmc.combine <- function( ... ){
  return( …

我正在尝试研究基于"做贝叶斯数据分析:R,JAGS和斯坦(2015)的教程"一书中的贝叶斯分析.

在本书中,有一些例子.所以,我试图在R中复制这个例子.但是,在这个例子中我收到了一条错误信息.

具体而言,这是示例数据.

data
   y        s
1  1 Reginald
2  0 Reginald
3  1 Reginald
4  1 Reginald
5  1 Reginald
6  1 Reginald
7  1 Reginald
8  0 Reginald
9  0     Tony
10 0     Tony
11 1     Tony
12 0     Tony
13 0     Tony
14 1     Tony
15 0     Tony

y<-data$y
s<-as.numeric(data$s)
Ntotal=length(y)
Nsubj=length(unique(s))

dataList=list(y=y, s=s, Ntotal=Ntotal, Nsubj=Nsubj)

另外,这是我的模特.

modelString=" 
model{
  for(i in 1:Ntotal){
    y[i] ~ dbern(theta[s[i]])
  }
  for(s in 1:Nsubj){
    theta[s] ~ dbeta(2,2)
  }
}
"
writeLines(modelString, con="TEMPmodel.txt") …

在 JAGS 中，我想为参数 w[i] 定义泊松分布，如果另一个参数 e[i] 大于 0，它也会被截断（大于或等于 2）。

基本上我希望它代表：

w[i] ~ ifelse( e[i] > 0, dpois(mu) T(2,) , dpois(mu) )

我尝试通过调整响应其他人的帖子而给出的代码来使用 step 函数，该帖子要求类似的内容： 根据 WinBugs/JAGS 中的 if - else 条件选择不同的发行版

但这似乎不起作用？

谢谢

我正在编写一个JAGS脚本(分层贝叶斯模型),其中事件的时间被建模为两个进程之间的竞争.

观察: time是事件的测量时间.

模型:具有高斯率的两个过程 - 无论哪个过程首先触发事件.

目标:估算两个过程的比率.

model{
  # Priors
  mu1 ~ dnorm( 0,1 )                      # rate of one process
  mu2 ~ dnorm( 0,1 )                      # rate of other process
  sigma1 <- 1                             # variability in rate
  sigma2 <- 0.1                           # variability in rate

  # Observations
  for (i in 1:N)
    rate1[i] ~ dnorm( mu1, sigma1 )        # Sample the two
    rate2[i] ~ dnorm( mu2, sigma2 )        # racing processes.

    rmax[i] <- max( rate1[i], rate2[i] ) …