JAGS：单位特定的时间趋势

Question

我正在尝试使用 JAGS 来估计一个包含单位特定时间趋势的模型。但是，问题是我不知道如何对此建模，到目前为止我一直无法找到解决方案。

例如，假设我们有以下数据：

rain<-rnorm(200)     # Explanatory variable
n1<-rnorm(200)       # Some noise
gdp<-rain+n1         # Outcome variable
ccode<-rep(1:10,20)  # Unit codes
year<-rep(1:20,10)   # Years

使用正态线性回归，我们将模型估计为：

m1<-lm(gdp~rain+factor(ccode)*year)

其中 factor(ccode)*year 是特定于单位的时间趋势。现在我想使用 JAGS 来估计模型。所以我为索引创建参数：

N<-200
J<-max(ccode) 
T<-max(year)

并估计模型，

library(R2jags)
library(rjags)

set.seed(42); runif(1)
dat<-list(gdp=gdp,
      rain=rain,
      ccode=ccode,
      year=year,
      N=N,J=J,T=T)

parameters<-c("b1","b0")
model.file <- "~/model.txt"
system.time(m1<-jags(data=dat,inits=NULL,parameters.to.save=parameters,
        model.file=model.file,
        n.chains=4,n.iter=500,n.burnin=125,n.thin=2))

使用以下模型文件，这是目前错误所在：

# Simple model 

model {
  # For N observations
  for(i in 1:N) {
    gdp[i] ~ dnorm(yhat[i], tau)
    yhat[i] <- b1*rain[i] + b0[ccode[i]*year[i]]   
  }

  for(t in 1:T) {
    for(j in 1:J) {
      b0[t,j] ~ dnorm(0, .01)
    }
  }
  # Priors
  b1 ~ dnorm(0, .01)   

  # Hyperpriors
  tau <- pow(sd, -2)
  sd ~ dunif(0,20)  
}

鉴于我在使用代码 Compilation error on line 7. Dimension mismatch taking subset of b0 时遇到的错误，我相当确定我定义 b0 和索引的方式是不正确的。但是，我不知道如何解决这个问题，所以我想知道这里是否有人对此有建议？

Answer 1

你的lm例子也可以这样写：

m1 <- lm(gdp ~ -1 + rain + factor(ccode) + factor(ccode):year)

等效的 JAGS 模型为：

M <- function() {
  for(i in 1:N) {
    gdp[i] ~ dnorm(yhat[i], sd^-2)
    yhat[i] <- b0[ccode[i]] + b1*rain[i] + b2[ccode[i]]*year[i]   
  }

  b1 ~ dnorm(0, 0.001)   
  for (j in 1:J) {
    b0[j] ~ dnorm(0, 0.001)
    b2[j] ~ dnorm(0, 0.001)
  }
  sd ~ dunif(0, 100)  
}

parameters<-c('b0', 'b1', 'b2')
mj <- jags(dat, NULL, parameters, M, 3)

比较系数：

par(mfrow=c(1, 2), mar=c(5, 5, 1, 1))
plot(mj$BUGSoutput$summary[grep('^b0', row.names(mj$BUGSoutput$summary)), '50%'],
     coef(m1)[grep('^factor\(ccode\)\d+$', names(coef(m1)))],
     xlab='JAGS estimate', ylab='lm estimate', pch=20, las=1,
     main='b0')
abline(0, 1)

plot(mj$BUGSoutput$summary[grep('^b2', row.names(mj$BUGSoutput$summary)), '50%'],
     coef(m1)[grep('^factor\(ccode\)\d+:', names(coef(m1)))],
     xlab='JAGS estimate', ylab='lm estimate', pch=20, las=1,
     main='b2')
abline(0, 1)

JAGS：单位特定的时间趋势

JAGS: unit-specific time trends

r

bayesian

jags