# This is code for the model in section 4.1.
# y[i,j] is the DV for the ith subject on the jth occasion.
# x1[i] is the experimental condition, with 0 = no prime and 1 = prime. 
model
	{
	for( i in 1 : N ) {
		for(j in 1:2){
			y[i,j] ~ dbeta(omega[i,j], tau[i,j])
# We reparameterize the beta distribution.
			omega[i,j] <- mu[i,j]*phi[i,j]
			tau[i,j] <- (1-mu[i,j])*phi[i,j]
			mu[i,j] <- exp(lambda[i,j])/(1+exp(lambda[i,j]))
# This is the dispersion submodel.
			phi[i,j] <- exp(kappa[T[i,j],j])
# This is the composition submodel. 
			T[i,j] ~ dcat(P[i,j,1:3])
			for (k in 1:3) {
				P[i,j,k] <- delta[i,j,k]
				}
				alpha[i,j,1] <- theta0[1,j] + theta1[1,j]*x1[i] 
				alpha[i,j,2] <- theta0[2,j] + theta1[2,j]*x1[i]
				delta[i,j,1] <- exp(alpha[i,j,1])/(1+exp(alpha[i,j,1])+exp(alpha[i,j,2]))
				delta[i,j,2] <- exp(alpha[i,j,2])/(1+exp(alpha[i,j,1])+exp(alpha[i,j,2]))
				delta[i,j,3] <- 1 - delta[i,j,1] - delta[i,j,2]
# This is the location submodel with the random effect for j = 3.
			lambda[i,j] <- beta0[T[i,j],j] + step(T[i,j]-3)*u1[i]
			}
		u1[i] ~ dnorm(0,tauu1)
		}
		for(j in 1:2){
		beta0[3,j] ~ dnorm(0.0, 1.0E-6)
		beta0[2,j] <- log(1)
		beta0[1,j] <- log(1/2)
		kappa[1,j] ~ dnorm(5.0,5.0)
		kappa[2,j] ~ dnorm(15.0,10.0)
		kappa[3,j] ~ dnorm(0.0, 1.0E-6)	
		theta0[1,j] ~ dnorm(0.0, 1.0E-6)
		theta0[2,j] ~ dnorm(0.0, 1.0E-6)
		theta1[1,j] ~ dnorm(0.0, 1.0E-6)
		theta1[2,j] ~ dnorm(0.0, 1.0E-6)
		dtheta0[j] <- theta0[j,2]-theta0[j,1]
		dtheta1[j] <- theta1[j,2]-theta1[j,1]
		}
		tauu1 ~ dgamma(0.001,0.001)
		sig1 <- 1/sqrt(tauu1)
	}

# Data
list(N = 155, y = structure(.Data = c(0.15,0.25,0.8,0.316666667,0.1,0.093333333,0.25,0.333333333,0.999,0.999,0.5,0.583333333,0.8,0.633333333,0.333333333,0.416666667,0.45,0.45,0.516666667,0.583333333,0.5,0.7,0.33,0.33,0.5,0.5,0.5,0.5,0.33,0.356666667,0.5,0.566666667,0.33,0.33,0.33,0.366666667,0.33,0.333333333,0.5,0.666666667,0.333333333,0.3,0.5,0.7,0.1,0.1,0.33,0.333333333,0.03,0.03,0.5,0.5,0.1,0.133333333,0.3,0.366666667,0.5,0.166666667,0.32,0.32,0.433333333,0.566666667,0.4,0.4,0.6,0.466666667,0.33,0.33,0.5,0.566666667,0.333333333,0.4,0.6,0.533333333,0.25,0.333333333,0.333333333,0.333333333,0.333,0.186666667,0.5,0.566666667,0.11,0.11,0.5,0.483333333,0.333333333,0.333333333,0.33,0.333333333,0.07,0.07,0.1,0.1,0.33,0.333333333,0.5,0.5,0.33,0.333333333,0.33,0.33,0.33,0.33,0.333333333,0.5,0.001,0.233333333,0.2,0.233333333,0.533333333,0.633333333,0.4,0.473333333,0.333,0.333,0.33,0.333333333,0.5,0.633333333,0.333,0.433333333,0.633333333,0.666666667,0.5,0.5,0.3,0.466666667,0.1,0.466666667,0.716666667,0.766666667,0.65,0.616666667,0.333333333,0.333333333,0.33,0.426666667,0.333,0.333,0.33,0.33,0.773333333,0.7,0.2,0.3,0.2,0.333333333,0.6,0.6,0.5,0.5,0.7,0.5,0.33,0.5,0.9,0.9,0.7,0.666666667,0.333333333,0.333333333,0.1,0.1,0.333333333,0.553333333,0.25,0.566666667,0.33,0.4,0.5,0.383333333,0.5,0.5,0.6,0.9,0.33,0.5,0.5,0.333333333,0.5,0.5,0.583333333,0.3,0.25,0.25,0.333,0.333,0.6,0.616666667,0.5,0.633333333,0.33,0.4,0.333333333,0.333333333,0.333,0.333,0.333333333,0.333333333,0.5,0.666666667,0.5,0.5,0.5,0.6,0.266666667,0.333333333,0.5,0.533333333,0.333333333,0.333333333,0.4,0.333333333,0.5,0.333333333,0.633333333,0.6,0.5,0.733333333,0.5,0.5,0.33,0.333333333,0.33,0.333333333,0.999,0.6,0.333333333,0.366666667,0.5,0.566666667,0.5,0.5,0.333333333,0.6,0.33,0.566666667,0.666666667,0.7,0.3,0.4,0.333333333,0.333333333,0.416666667,0.633333333,0.25,0.25,0.383333333,0.483333333,0.333333333,0.333333333,0.33,0.333333333,0.5,0.616666667,0.566666667,0.533333333,0.3,0.233333333,0.6,0.75,0.3,0.566666667,0.25,0.3,0.516666667,0.733333333,0.333333333,0.333333333,0.4,0.733333333,0.1,0.2,0.5,0.5,0.5,0.6,0.7,0.766666667,0.2,0.3,0.33,0.4,0.33,0.333333333,0.333333333,0.36,0.3,0.4,0.6,0.533333333,0.5,0.8,0.683333333,0.7,0.9,0.733333333,0.3,0.5,0.3,0.333333333,0.5,0.6,0.1,0.1,0.5,0.5,0.2,0.233333333), .Dim = c(155,2)), x1 = c(1,0,0,0,0,0,1,0,0,1,1,0,0,0,1,0,1,1,1,0,1,1,0,1,0,1,1,1,0,0,1,0,0,1,0,1,0,0,1,1,1,0,0,1,1,1,1,1,1,1,0,0,1,0,0,1,1,1,1,1,1,0,1,0,1,1,0,1,0,1,1,1,0,0,0,0,0,0,1,1,1,1,0,0,0,0,0,0,1,1,0,0,0,1,1,1,1,1,1,1,0,1,0,1,1,1,1,1,1,1,0,0,0,0,1,0,0,1,1,1,1,1,1,0,1,1,0,0,0,0,0,0,1,0,0,0,0,0,0,1,0,1,1,1,1,1,0,0,0,0,0,0,0,0,0))

#Inits for one chain
list(beta0 = structure(.Data = c(NA,NA,NA,NA,0.1,-0.2), .Dim = c(3,2)), kappa = structure(.Data = c(5.0,9.0,15.0,17.0,0.7,0.5), .Dim = c(3,2)), theta0 = structure(.Data = c(-1.0,-0.5,1.0,0.5), .Dim = c(2,2)), theta1 = structure(.Data = c(1.0,-0.5,1.0,-0.5), .Dim = c(2,2)), tauu1 = 0.25, T = structure(.Data = c(3,3,3,3,3,3,3,1,3,3,2,3,3,3,1,3,3,3,3,3,2,3,1,1,2,2,2,2,1,3,2,3,1,1,1,3,1,1,2,3,1,3,2,3,3,3,1,1,3,3,2,2,3,3,3,3,2,3,3,3,3,3,3,3,3,3,1,1,2,3,1,3,3,3,3,1,1,1,1,3,2,3,3,3,2,3,1,1,1,1,3,3,3,3,1,1,2,2,1,1,1,1,1,1,1,2,3,3,3,3,3,3,3,3,1,1,1,1,2,3,1,3,3,3,2,2,3,3,3,3,3,3,3,3,1,1,1,3,1,1,1,1,3,3,3,3,3,1,3,3,2,2,3,2,1,2,3,3,3,3,1,1,3,3,1,3,3,3,1,3,2,3,2,2,3,3,1,2,2,1,2,2,3,3,3,3,1,1,3,3,2,3,1,3,1,1,1,1,1,1,2,3,2,2,2,3,3,1,2,3,1,1,3,1,2,1,3,3,2,3,2,2,1,1,1,1,3,3,1,3,2,3,2,2,1,3,1,3,3,3,3,3,1,1,3,3,3,3,3,3,1,1,1,1,2,3,3,3,3,3,3,3,3,3,3,3,3,3,1,1,3,3,3,3,2,2,2,3,3,3,3,3,1,3,1,1,1,3,3,3,3,3,2,3,3,3,3,3,3,2,3,1,2,3,3,3,2,2,3,3), .Dim = c(155,2)), u1 = c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0))