How to add a spline to rjags model

Question

I am having difficulty finding information in fitting splines using rjags (my motivation is to try to recreate a glm in jags to impute missing dependent values). Anyhow I can find very little info on this finding only this answer on Cross-validated: https://stats.stackexchange.com/questions/79973/how-to-analyze-this-data-using-rjags-or-any-other-way/80650#80650. However, I cannot understand the splines code there (and do not have the reputation there to ask a question!). For one thing I don't understand why that code loops over both S & G.

Therefore, I have made a toy linear model in jags:

library(datasets)
library(rjags)
library(ggplot2)


# Specify a JAGS linear model

mk_jags_lin_mod <- function(prior.a, prior.b){
    sink(paste("lin_reg_jags.mod.txt", sep=""))
    cat(paste0("
               model {
               for (i in 1:N){
               y[i] ~ dnorm(y.hat[i], tau)
               y.hat[i] <- a + b * x[i]
               }
               a ~ ",prior.a,
               "\tb ~ ",prior.b,
               "\ttau <- pow(sigma, -2)
               sigma ~ dunif(0, 100)
               }
               "))
    sink()
}

# Define a default vague prior
default <- "dnorm(0, .0001)\n"

mk_jags_lin_mod(default, default)

# Initialise
jags.cars <- jags.model('lin_reg_jags.mod.txt',
                              data = list('x' = mtcars$hp,
                                          'y' = mtcars$mpg,
                                          'N' = nrow(mtcars)),
                              n.chains = 2,
                              n.adapt = 1000)

# Burn-in
update(jags.cars, 5000) 

# Sample
coda.cars <- coda.samples(jags.cars, variable.names = c('a', 'b', 'y.hat','tau'), n.iter = 1000)

# Extract posterior estimates
coda.sum <- summary(coda.cars)
q <- coda.sum$quantiles

mtcars$fit <- q[4:35 , 3]

ggplot(data=mtcars, aes(x=hp, y=mpg)) + geom_point() +
    geom_line(aes(y=fit, col="red"))

My question is: how to I add a restricted cubic spline to the 'b' estimate in the linear model ?

Answer 1

Ok so it took me most of the day but I finally figured it out... I think.... I modified the model as follows:

 model {
   for (i in 1:N){
     y[i] ~ dnorm(y.hat[i], tau)
     y.hat[i] <- a + beta[1] + beta[2]*x[i] + beta[3]*pow(x[i], 2) + beta[4]*pow(x[i], 3)
   }

   a ~ dnorm(0, .0001)

   # Specify priors for spline terms
   for (k in 1:4) {
     beta.mu[k] ~ dnorm(0, 100)
     beta.tau[k] ~ dgamma(0.01, 10)
     beta[k] ~ dnorm(beta.mu[k], beta.tau[k])
   }
     tau <- pow(sigma, -2)
     sigma ~ dunif(0, 100)
   }

The value of k must be at least 4 but thereafter increasing k increases the smooth. Otherwise I think my math is correct (though I am open to corrections!).