svss <- "model{

for(i in 1:n){
y[i] ~ dnorm(yhat[i], shat[i])
yhat[i] <-  (u[ID[i], 1]) + (X[i] * u[ID[i], 2])
shat[i] <-  exp(u[ID[i], 3] + (X[i] * u[ID[i], 4]))
}

for(i in 1:J){
u[i,1:4] ~ dmnorm.vcov(bhat[i,], Sigma[1:4,1:4])
bhat[i,1] <- beta[1]
bhat[i,2] <- beta[2]
bhat[i,3] <- beta[3]
bhat[i,4] <- beta[4]
}


beta[1] ~ dnorm(0, 1)
beta[2] ~ dnorm(0, pow(0.25, - 2))
beta[3] ~ dnorm(0, 1)
beta[4] ~ dnorm(0, 1)

sigma_11 ~  dgamma(1, 1/5)
sigma_22 ~  dgamma(1, 1/5)
sigma_33 ~  dgamma(1, 1/5)
sigma_44 ~  dgamma(1, 1/5)

Sigma[1,1] <- pow(sigma_11, 2)
Sigma[2,2] <- pow(sigma_22, 2)
Sigma[3,3] <- pow(sigma_33, 2)
Sigma[4,4] <- pow(sigma_44, 2)

Sigma[1,2] <- (rhoNS[1] * sigma_11 * sigma_22)
Sigma[2,1] <- Sigma[1,2]

Sigma[1,3] <- (rho[2] * sigma_11 * sigma_33)
Sigma[3,1] <- Sigma[1,3]

Sigma[1,4] <- (rho[3] * sigma_11 * sigma_44)
Sigma[4,1] <- Sigma[1,4]


Sigma[2,3] <- (rho[4] * sigma_22 * sigma_33)
Sigma[3,2] <- Sigma[2,3]


Sigma[2,4] <- rho[5] * sigma_22 * sigma_33
Sigma[4,2] <- Sigma[2,4]


Sigma[3,4] <- (rhoNS[2] * sigma_33 * sigma_44)
Sigma[4,3] <- Sigma[3,4]



for(i in 2:5){
IndA[i] ~ dcat(PInd[])
Ind[i] <- IndA[i] - 1
z[i] ~ dnorm(0, TauM[IndA[i]])
rho[i] <- tanh(z[i])
}

zNS[1] ~ dnorm(0, pow(0.5, -2))
zNS[2] ~ dnorm(0, pow(0.5, -2))
rhoNS[1] <- tanh(zNS[1])
rhoNS[2] <- tanh(zNS[2])


TauM[1] <- pow(0.01, -2);
TauM[2] <- pow(0.5, - 2);
PInd[1] <- 0.5;
PInd[2] <- 0.5;

}"