JAGS线性AR1预测模型具有嵌套组

Question

我正在尝试用平衡的嵌套组编写AR1线性模型，以预测因变量中的结果。例如，假设我们有两组每组三个观测值，其中每组中的第三个观察结果都缺少。目的是编写一个模型，该模型使用结果变量中的先前观察结果预测每个组中的第三个观察结果。我很难找出编写嵌套模型的正确方法。下面的代码给出了一个错误，“索引超出了组的子集”。任何帮助将不胜感激。

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)

Answer 1

怎么样：

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

^{在2022-07-07创建的 preprex package （v2.0.1）< /sup>}

如果您有平衡的数据，最简单的事情是将y放在n_obs中x n_groups矩阵。然后，AR（1）部分变得容易。