JAGS:特定于单位的时间趋势 [英] JAGS: unit-specific time trends

问题描述

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

Using JAGS I am trying to estimate a model including a unit-specific time trend. However, the problem is that I don't know how to model this and so far I have been unable to find a solution.

作为示例，请考虑我们具有以下数据:

As an example, consider we have the following data:

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 

使用正态线性回归，我们可以将模型估算为:

Using normal linear regression, we would estimate the model as:

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

其中factor(ccode)*year是单位特定的时间趋势.现在，我想使用JAGS估算模型.因此，我为索引创建了参数:

Where factor(ccode)*year is the unit-specific time trend. Now I want to estimate the model using JAGS. So I create parameters for the indexing:

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))

具有以下模型文件，这是当前错误所在:

with the following model file, and this is where the error is at the moment:

# 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和建立索引的方式是不正确的. 但是，我不知道如何解决这个问题，所以我想知道这里有人对此有建议吗?

I am fairly sure that the way in which I define b0 and the indexing is incorrect given the error I get when using the code: Compilation error on line 7. Dimension mismatch taking subset of b0. However, I don't know how to solve this so I wondered whether someone here has suggestions about this?

推荐答案

您的lm示例也可以编写为:

Your lm example can also be written:

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

等效的JAGS模型为:

The equivalent JAGS model would be:

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:特定于单位的时间趋势的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！