Bayesian Interrupted Time Series

[hlsmith]

@Dason

I have questions regarding some JAGS code I am trying to adapt. I will start off with trying to change the listed priors. Below is the code I am using. This fits a time series that is interrupted. So outputs an intercept, slope and second intercept and second slope as well as the autocorrelation and standard deviation error terms.

Toy example

library(runjags)
y <- c(85, 85, 83, 83, 82, 83, 84, 85, 84, 84, 94, 96, 96, 98, 98, 97, 97, 98, 100, 100)


str(Y)
T <- length(y)
P <- 2
Tb <- 5

beta1 <- mean(y[1:Tb])
beta1
beta2 <- mean(y[(Tb + 1):T])
beta2
#Model is defined:
BITS.model2 <- "model {
yhat[1] <- beta[1, 1]             #INTERCEPT
yhat[Tb+1] <- beta[2, 1]          #INTERCEPT

#baseline phase
for (i in 2:Tb) {

#slope beta[1, 2] is added to the equation
yhat[i] <- beta[1, 1] + beta[1, 2] * i
#MODEL: ESTIMATE + AUTOCORRELATION * ERROR AND PRECISION
y[i] ~ dnorm(yhat[i] + rho * (y[i - 1] - yhat[i - 1]), tau)
}



#treatment phase
for (i in (Tb + 2):T) {

#slope is multiplied by (T-Tb) as given in equation 7 such that
#(i-Tb)  = x represents the xth time-point in the intervention phase
yhat[i] <- beta[2, 1] + beta[2, 2] * (i - Tb)
#MODEL FOR POST PERIOD
y[i] ~ dnorm(yhat[i] + rho * (y[i - 1] - yhat[i - 1]), tau)
}
y[1] ~ dnorm(yhat[1], tau)
y[(Tb + 1)] ~ dnorm(yhat[(Tb + 1)], tau)

#PRIOR specifications, P=number of periods
for (i in 1:P){
#the slopes for baseline and treatment phases are drawn
#from distributions with means mu.s[1] and mu.s[2], respectively
# and precision 0.01
beta[i, 1] ~ dnorm(mu[i], 0.01)      #Intercept
beta[i, 2] ~ dnorm(mu.s[i], 0.01)    #Slope
mu[i] ~ dnorm(40, 0.5)               #Intercept
mu.s[i] ~ dnorm(0, 0.1)              #Slope
}
sigma ~ dunif(0.1, 5)                #Model standard deviation uniformly vary between 0.1 to 5.
tau <- pow(sigma, -2)                #Precision, T=1/variance, pow does this
rho ~ dunif(-1, 1)                   #Auto-correlation term
}"
############# end of model definition##########################

#Begin running the model with the data

results <- autorun.jags(
  model = BITS.model2,
  data = list(y = y, T = T, P = P, Tb = Tb),
  monitor = c("beta", "sigma", "rho"),
  n.chains = 2, #4 increase once model runs
  startsample = 10000, #30000 increase once model runs
  inits = function() {
    list(
      beta = rbind(c(rnorm(1, beta1, 1), 1),
                   c(rnorm(1, beta2, 1), 1)),
      sigma = runif(1, 0.1, 5),
      rho = runif(1, -1, 1)
    )
  },
  method = "rjparallel"
  #run the chains in parallel
)


results$draws <- combine.mcmc(results$mcmc)
results

 

[hlsmith]

Initial questions, I am not seeing Priors for the two intercepts and slopes. My brain sees the following and does know how to input two slopes and intercepts. Is the code just recycling this for use in both models?

beta[i, 1] ~ dnorm(mu, 0.01) #Intercept
beta[i, 2] ~ dnorm(mu.s, 0.01) #Slope
mu ~ dnorm(40, 0.5) #Intercept
mu.s ~ dnorm(0, 0.1) #Slope


In my actual project my priors would be:
Intercept #1: dnorm(40, 2)
Slope #1: dnorm(0, 0.0005)
Intercept #2: dnorm(8, 1)
Slope #2: dnorm(0, 0.0005)

 

[hlsmith]

Yes there is a loop. I need to write the code so the priors don't get feed into both models. Any suggestions?

 

[hlsmith]

I switched the priors to the following, which I think addresses the same priors getting used for both model portions. Any thoughts?

#PRIOR specifications, P=number of periods

#the slopes for baseline and treatment phases are drawn 
#from distributions with means mu.s[1] and mu.s[2], respectively 
# and precision 0.01
beta[1, 1] ~ dnorm(mu[1], 2.0)      #Intercept
beta[1, 2] ~ dnorm(mu.s[1], 2.0)    #Slope
mu[1] ~ dnorm(40, 0.5)              #Intercept
mu.s[1] ~ dnorm(0, 0.1)             #Slope

beta[2, 1] ~ dnorm(mu[2], 2.0)      #Intercept
beta[2, 2] ~ dnorm(mu.s[2], 2.0)    #Slope
mu[2] ~ dnorm(8, 0.5)               #Intercept
mu.s[2] ~ dnorm(0, 0.1)             #Slope
sigma ~ dunif(0.1, 5)               #Model standard deviation uniformly vary between 0.1 to 5.
tau <- pow(sigma, -2)               #Precision, T=1/variance, pow does this
rho ~ dunif(-1, 1)                  #Auto-correlation term
}"