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

Я достаточно уверен, что способ, которым я определяю b0, и индексирование неверно, учитывая ошибку, которую я получаю при использовании кода: Compilation error on line 7. Dimension mismatch taking subset of 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)