model;
{
   ############ 1st condition
   for( i in 1 : n ) {
	 for( j in 1 : nreps) {
       y1[i, j] ~ dnorm(x1[i, j], tau1[i])
       ynew1[i, j] ~ dnorm(x1[i, j], taunew1[i])
       ## equal variances 
       #y1[i, j] ~ dnorm(x1[i, j], tau1)
       #ynew1[i, j] ~ dnorm(x1[i, j], tau1)
       x1[i, j] <- alpha[i] - 0.5*dstar[i] + c1[i, j]
	 }
       ##### predictive p-values
	 s21[i] <- pow(sd(y1[i, ]), 2)
       s2new1[i] <- pow(sd(ynew1[i, ]), 2)
       pval1[i] <- step(s2new1[i] - s21[i])
   } 
   ###### exch variances
   for( i in 1 : n ) {
       tau1[i] <- 1.0/sig21[i]
       sig21[i] <- exp(lsig21[i])
       lsig21[i] ~ dnorm(mm1,tt1) 
       sig2new1[i] <- exp(lsig2new1[i])
       taunew1[i] <- 1.0/sig2new1[i]
       lsig2new1[i] ~ dnorm(mm1,tt1) 
   }
   mm1 ~ dnorm( 0.0,1.0E-3)
   tt1 ~ dgamma(0.01,0.01)
   ######### equal variances
   #tau1 <- 1.0/sig21
   #sig21 <- exp(lsig21)
   #lsig21 ~ dnorm( 0.0, 0.001)  

for( i in 1 : n ){
	 for( j in 1 : nreps  ){
	   cdum1[i, j] <- b1[j] * pow(alpha[i]-a1[j],2) * step(alpha[i]-a1[j])
	   c1[i, j] <- cdum1[i, j] + beta10[j] + beta11[j] * (alpha[i]-ma) + 
	                 beta12[j]  * pow(alpha[i]-ma,2) * step(alpha[i]-ma)
	 }
   }

   ############ 2nd condition
   for( i in 1 : n ) {
	 for( j in 1 : nreps) {
       y2[i, j] ~ dnorm(x2[i, j], tau2[i])
       ynew2[i, j] ~ dnorm(x2[i, j], taunew2[i])
       ## equal variances
       #y2[i, j] ~ dnorm(x2[i, j], tau2)
       #ynew2[i, j] ~ dnorm(x2[i, j], tau2)
       x2[i, j] <- alpha[i] + 0.5*dstar[i] + c2[i, j]
	 }
       ##### predictive p-values
	 s22[i] <- pow(sd(y2[i, ]), 2)
       s2new2[i] <- pow(sd(ynew2[i, ]), 2)
       pval2[i] <- step(s2new2[i] - s22[i])
   } 
   ######## exch variances
   for( i in 1 : n ) {
       tau2[i] <- 1.0/sig22[i]
       sig22[i] <- exp(lsig22[i])
       lsig22[i] ~ dnorm(mm2,tt2) 
       taunew2[i] <- 1.0/sig2new2[i]
       sig2new2[i] <- exp(lsig2new2[i])
       lsig2new2[i] ~ dnorm(mm2,tt2) 
   }
   mm2 ~ dnorm( 0.0,1.0E-3)
   tt2 ~ dgamma(0.01,0.01)
   ####### equal variances
   #tau2 <- 1.0/sig22
   #sig22 <- exp(lsig22)
   #lsig22 ~ dnorm( 0.0, 0.001) 


  for( i in 1 : n ){
	 for( j in 1 : nreps  ){
	   cdum2[i, j] <- b2[j] * pow(alpha[i]-a2[j],2) * step(alpha[i]-a2[j])
	   c2[i, j] <- cdum2[i, j] + beta20[j] + beta21[j] * (alpha[i]-ma) + 
	               beta22[j]  * pow(alpha[i]-ma,2) * step(alpha[i]-ma)
	 }
   }

 ############# both conditions together

 for( i in 1 : n ){
    alpha[i] ~ dunif(ma,mb)
	dstar[i] ~ dnorm(0.0,1.0E-4)
	dstarprime[i] <- dstar[i] - dstarmean
      #d[i] <- dstar[i]/sqrt((sig21[i] + sig22[i])/nreps)
	dcorrect[i] <- dstarprime[i]/sqrt((sig21[i] + sig22[i])/nreps)
	pd[i] <- step(dstarprime[i] - log(3))
 }
 dstarmean <- mean(dstar[])

############ priors on coefs
 for( j in 1 : nreps ) {
   beta10[j] ~ dnorm(0.0,1.0E-1)
   beta11[j] ~ dnorm(0.0,1.0E-1)
   beta12[j] ~ dnorm(0.0,1.0E-1)
   beta13[j] ~ dnorm(0.0,1.0E-1)
   beta20[j] ~ dnorm(0.0,1.0E-1)
   beta21[j] ~ dnorm(0.0,1.0E-1)
   beta22[j] ~ dnorm(0.0,1.0E-1)
   beta23[j] ~ dnorm(0.0,1.0E-1)
   b1[j] ~ dnorm(0.0,1.0E-1) 
   b2[j] ~ dnorm(0.0,1.0E-1) 
 } 

######## knots
 for( j in 1 : 3 ) {
     a1[j] ~ dunif(ma,mb)
     a2[j] ~ dunif(ma,mb)
 }

######### constraints using dummy data
for( i in 1 : n ) {
    scr1[i] <- sum(c1[i, ])
    scr2[i] <- sum(c2[i, ])
    zc1row[i] ~ dnorm( scr1[i], 1.0E+6 )
    zc2row[i] ~ dnorm( scr2[i], 1.0E+6 )
}

}