0

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。しかし、これを解決する方法がわからないので、ここで誰かがこれについて提案を持っているかどうか疑問に思いましたか?

4

1 に答える 1

2

あなたの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)

ここに画像の説明を入力

于 2015-01-27T09:39:52.753 に答える