@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
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