I'm trying to write an AR1 linear model in JAGS with balanced, nested groups, that forecasts outcomes in the dependent variable. As an example, suppose we have two groups with three observations per group where the third observation in each group is missing. The objective is to write a model that forecasts the third observation in each group using previous observations in the outcome variable. I'm having trouble figuring out the correct way to write the nested model. The code below gives an error, "Index out of range taking subset of group." Any help is greatly appreciated.
data<-data.frame(group=c(1, 1, 1, 2, 2, 2), y=c(3, 5, NA, 10, 2, NA))
model<-function(){
for(i in 1:6){
y[group[i]]~dnorm(mu[group[i]], tau)
mu[group[i]]<-beta[1] + beta[2]*y[group[i-1]]
}
for(i in 1:2){
beta[i]~dnorm(0, .01)
}
tau~dgamma(.01, .01)
}
model.data<-list("group", "y")
model.data<-list(group=data$group, y=data$y)
model.params<-c("y")
model.fit<-jags(data=model.data, inits=NULL, model.params, n.chains=2, n.iter=10000, n.burnin=1000, model.file=model)
How about this:
data<-data.frame(group=c(1, 1, 1, 1, 1, 2, 2, 2, 2, 2), y=c(3, 6, 5, 7,NA, 2,5,3,4, NA))
library(R2jags)
model<-function(){
for(i in 2:5){
for(j in 1:2){
y[i,j] ~ dnorm(mu[i,j], tau)
mu[i,j] <- beta[1] + beta[2]*y[(i-1),j]
}
}
for(i in 1:2){
beta[i]~dnorm(0, .01)
}
tau~dgamma(.01, .01)
}
model.data <- list(y = matrix(data$y, ncol=2))
model.params<-c("y", "beta")
model.fit<-jags(data=model.data, inits=NULL, model.params, n.chains=2, n.iter=10000, n.burnin=1000, model.file=model)
#> module glm loaded
#> Compiling model graph
#> Resolving undeclared variables
#> Allocating nodes
#> Graph information:
#> Observed stochastic nodes: 6
#> Unobserved stochastic nodes: 5
#> Total graph size: 27
#>
#> Initializing model
up <- update(model.fit)
up$BUGSoutput
#> Inference for Bugs model at "/var/folders/qy/y5n2dh2x24d_p19jtdv2_d5w0000gn/T//Rtmpa3onfy/modele7895ec19988.txt", fit using jags,
#> 2 chains, each with 1000 iterations (first 0 discarded)
#> n.sims = 2000 iterations saved
#> mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff
#> beta[1] 4.7 2.5 -0.4 3.3 4.8 6.2 9.6 1 2000
#> beta[2] 0.1 0.6 -1.1 -0.3 0.0 0.4 1.2 1 2000
#> deviance 24.1 3.2 20.4 21.8 23.3 25.5 32.6 1 2000
#> y[1,1] 3.0 0.0 3.0 3.0 3.0 3.0 3.0 1 1
#> y[2,1] 6.0 0.0 6.0 6.0 6.0 6.0 6.0 1 1
#> y[3,1] 5.0 0.0 5.0 5.0 5.0 5.0 5.0 1 1
#> y[4,1] 7.0 0.0 7.0 7.0 7.0 7.0 7.0 1 1
#> y[5,1] 5.1 2.9 -0.7 3.5 5.1 6.8 10.9 1 2000
#> y[1,2] 2.0 0.0 2.0 2.0 2.0 2.0 2.0 1 1
#> y[2,2] 5.0 0.0 5.0 5.0 5.0 5.0 5.0 1 1
#> y[3,2] 3.0 0.0 3.0 3.0 3.0 3.0 3.0 1 1
#> y[4,2] 4.0 0.0 4.0 4.0 4.0 4.0 4.0 1 1
#> y[5,2] 5.1 2.3 0.6 3.8 5.1 6.4 10.5 1 2000
#>
#> For each parameter, n.eff is a crude measure of effective sample size,
#> and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
#>
#> DIC info (using the rule, pD = var(deviance)/2)
#> pD = 5.1 and DIC = 29.2
#> DIC is an estimate of expected predictive error (lower deviance is better).
Created on 2022-07-07 by the reprex package (v2.0.1)
If you've got balanced data, the easiest thing to do is to put y in a n_obs x n_groups matrix. Then, the AR(1) part becomes easy.