Bayesian Interrupted Time Series

hlsmith

Less is more. Stay pure. Stay poor.
#1
@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
Code:
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

Less is more. Stay pure. Stay poor.
#2
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

Less is more. Stay pure. Stay poor.
#3
Yes there is a loop. I need to write the code so the priors don't get feed into both models. Any suggestions?
 

hlsmith

Less is more. Stay pure. Stay poor.
#4
I switched the priors to the following, which I think addresses the same priors getting used for both model portions. Any thoughts?

Code:
#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
}"