JAGS and R: Obtain posterior predictive distribution for specific x

I am trying to obtain a posterior predictive distribution for specified values of x from a simple linear regression in Jags. I could get the regression itself to work by adapting this example (from https://biometry.github.io/APES//LectureNotes/StatsCafe/Linear_models_jags.html) to my own data. I have supplied a smal chunk of this data here so that the code works here also.

library(rjags)
library(R2jags)

#create data
dw=c(-15.2,-13.0,-10.0,-9.8,-8.5,-8.5,-7.7,-7.5,-7.2,-6.1,-6.1,-6.1,-5.5,-5.0,-5.0,-5.0,-4.5,-4.0,-2.0,-1.0,1.3)
phos=c(11.8,12.9,15.0,14.4,17.3,16.1,20.8,16.6,16.2,18.2,18.8,19.2,15.6,17.0,18.9,22.1,18.9,22.8,21.6,20.5,21.1)

#convert to list
jagsdwphos=list(dw=dw,phos=phos,N=length(phos))

#write model function for linear regression
lm1_jags <- function(){
    # Likelihood:
    for (i in 1:N){
    phos[i] ~ dnorm(mu[i], tau) # tau is precision (1 / variance)
    mu[i] <- intercept + slope * dw[i]
  }
    # Priors:
  intercept ~ dnorm(0, 0.01)
  slope ~ dnorm(0, 0.01)
  sigma ~ dunif(0, 100) # standard deviation
  tau <- 1 / (sigma * sigma)
}

#specifiy paramters of MCMC sampler, choose posteriors to be reported and      run the jags model
#set initial values for MCMC
init_values <- function(){
  list(intercept = rnorm(1), slope = rnorm(1), sigma = runif(1))
}
#choose paramters to report on
params <- c("intercept", "slope", "sigma")
#run model in jags
lm_dwphos <- jags(data = jagsdwphos, inits = init_values, parameters.to.save = params, model.file = lm1_jags,
          n.chains = 3, n.iter = 12000, n.burnin = 2000, n.thin = 10, DIC = F)

In addition to this regression, I would like to have an output of the posterior predictive distributions of particular phos values, but I cannot get it to work with this simple example I have written. I found a tutorial here https://doingbayesiandataanalysis.blogspot.com/2015/10/posterior-predicted-distribution-for.html and tried to implement it like this:

#create data
dw=c(-15.2,-13.0,-10.0,-9.8,-8.5,-8.5,-7.7,-7.5,-7.2,-6.1,-6.1,-6.1,-5.5,-5.0,-5.0,-5.0,-4.5,-4.0,-2.0,-1.0,1.3)
    phos=c(11.8,12.9,15.0,14.4,17.3,16.1,20.8,16.6,16.2,18.2,18.8,19.2,15.6,17.0,18.9,22.1,18.9,22.8,21.6,20.5,21.1)
#specifiy phos values to use for posterior predictive distribution
phosprobe=c(14,18,22)

#convert to list
jagsdwphos=list(dw=dw,phos=phos,N=length(phos),xP=phosprobe)

#write model function for linear regression
lm1_jags <- function(){
  # Likelihood:
  for (i in 1:N){
    phos[i] ~ dnorm(mu[i], tau) # tau is precision (1 / variance)
    mu[i] <- intercept + slope * dw[i]
  }
  # Priors:
  intercept ~ dnorm(0, 0.01) # intercept
  slope ~ dnorm(0, 0.01) # slope
  sigma ~ dunif(0, 100) # standard deviation
  tau <- 1 / (sigma * sigma) # sigma^2 doesn't work in JAGS
  nu <- nuMinusOne+1
  nuMinusOne ~ dexp(1/29.0)
  #prediction
  for(i in 1:3){
    yP ~ dt(intercept+slope*xP[i],tau,nu)
  }
}

#specifiy paramters of MCMC sampler, choose posteriors to be reported and run the jags model
#set initial values for MCMC
init_values <- function(){
  list(intercept = rnorm(1), slope = rnorm(1), sigma = runif(1))
}
#choose paramters to report on
params <- c("intercept", "slope", "sigma","xP","yP")
#run model in jags
lm_dwphos <- jags(data = jagsdwphos, inits = init_values, parameters.to.save = params, model.file = lm1_jags,
          n.chains = 3, n.iter = 12000, n.burnin = 2000, n.thin = 10, DIC = F)

But I get the following error message:

Error in jags.model(model.file, data = data, inits = init.values, n.chains = n.chains, : RUNTIME ERROR: Compilation error on line 14. Attempt to redefine node yP[1]

I confess I don't quite understand exactly how the prediction was implemented in that example I used and could not find an explanation on what exactly nu is or where those numbers come from. So I presume that is where I made some mistake adapting to my example, but it was the only tutorial in Jags that I could find that gives the whole distribution of y values for the probed x instead of just the mean.

I would appreciate any help or explanation.

Thanks!

Solution

This error occurs because you are not indexing yP. You have written this loop like this:

  #prediction
  for(i in 1:3){
    yP ~ dt(intercept+slope*xP[i],tau,nu)
  }

As i moves from 1 to 3 the element yP is being written over. You need to index it like you have done with xP.

  #prediction
  for(i in 1:3){
    yP[i] ~ dt(intercept+slope*xP[i],tau,nu)
  }