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：單位特定的時間趨勢

回答

相關問題