3

我正在学习如何在 R 中拟合贝叶斯多项逻辑模型。这是我第一次尝试通过rjags. 该代码用 MWE 说明了我正在尝试做的事情：

## simulate data
set.seed(123)

n=2000

rr<-rmultinom(n, 3, c(.1,.3,.2))

r2=as.numeric(rr==1)
r3=as.numeric(rr==2)
r4=as.numeric(rr==3)

abt=rbinom(n,1,.1);smk=rbinom(n,1,.3) 
age=rnorm(n);bmi=rnorm(n)

## load programs
library("rjags")

## model
NMMmodel.string <- "
  model{
      for (i in 1:N){ 

        ## outcome levels 2, 3, and 4
        r2[i] ~ dbern(pi2[i]) 
        r3[i] ~ dbern(pi3[i]) 
        r4[i] ~ dbern(pi4[i]) 

        ## linear predictors
        logit(pi2[i]) <- g[1]+g[2]*age[i]+g[3]*abt[i]+g[4]*smk[i]
        logit(pi3[i]) <- g[5]+g[6]*bmi[i]+g[7]*age[i]+g[8]*smk[i]
        logit(pi4[i]) <- g[9]+g[10]*age[i]+g[11]*smk[i]+g[12]*bmi[i]

        ## probability that outcome is level 1
        pi1[i]  <- 1-pi2[i]-pi3[i]-pi4[i]
      } 
      for (j in 1:12) {
        g[j] ~ dnorm(0, 0.01)
      }
    }
"
NMMmodel.spec<-textConnection(NMMmodel.string)

## fit model w JAGS
jags <- jags.model(NMMmodel.spec,
                   data = list('r2'=r2,'r3'=r3,'r4'=r4,
                               'abt'=abt,'smk'=smk,
                               'age'=age,'bmi'=bmi,'N'=n),
                   n.chains=4,
                   n.adapt=100)

以下是两个问题，按重要性递减顺序排列：

问题 1：我想对由g[1]to索引的估计参数施加一个约束，g[12]使其pi1位于某个任意上限和下限之间：例如，a=0.25和b=0.75。一种方法是使用拒绝采样，其中 rjags 将拒绝所有返回pi1小于a或大于 的样本b。我怎样才能做到这一点？

问题 2：这个程序到底在做什么？例如，如果这个程序实现了一个 Gibbs 采样器，有没有办法在不使用 JAGS、STAN 或 BUGS 的情况下对其进行编码？类似于这个网站上的第一组代码？

4

0 回答 0