JAGS：单位特定的时间趋势

使用JAGS我试图估计一个包含单位特定时间趋势的模型。但是，问题是我不知道如何建模，到目前为止我一直无法找到解决方案。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) 
}

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

来源

2015-01-15 BlankUsername

你在'yhat'行定义'b0'是一个简单的向量（只有一个维度）。在模型的后面，'b0'似乎是一个矩阵，有两个维度。这会导致错误。 – nicola 2015-01-15 17:12:14

这是否解决了您的问题？ – jbaums 2015-01-22 10:03:51

听起来不错，但到目前为止，我一直无法正确调整代码。正如我所说，我的造型知识不是很好，所以我仍然有点卡住。 – BlankUsername 2015-01-22 10:30:15

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

enter image description here

来源

2015-01-27 09:39:52 jbaums

这是非常有帮助和信息。谢谢。 – BlankUsername 2015-01-27 15:09:32

JAGS：单位特定的时间趋势

回答

相关问题