WinBUGS code for the article Simultaneous estimation of multiple quantitative trait loci and growth curve parameters through hierarchical Bayesian modeling. Heredity 108: 134-146.

Simulated data
Phenotypic residual variance 0.1
no block effect
Marker inits in mark inits.txt file
to monitor:
- alpha1[k], alpha2[k], alpha3[k]
- x1[k], x2[k], x3[k]
- NQTL1, NQTL2, NQTL3
- her[1], her[2], her[3]
- gcov12, 13, 23
- gcor12, 13, 23
- sigma
- intercept1, intercept2, intercept3
- rcov12, rcov13, rcov23

model
{

for(j in 1 : loc) {	 
   for(k in 1 : 2 ){	 
      p[ j, k ] <- 1/2		
   }	  
}
	 
for(i in 1 : ind) {	 
   for(j in 1 : loc ) {	 
      mark[i, j ]  ~ dcat( p[ j, 1:2] )	
   }	 
}

for(i in 1 : ind) {	 
   for(j in 1 : loc ) {	  
      Z[i, j ] <- mark[i,j] - 1	
   }	 
}


for (i in 1:ind) {
		for(j in 2:years) {
			Y[i, j] ~ dnorm(mu[i,j], 10) 
			mu[i, j] <- beta[i, 1] + beta[i, 2] * age[j] + beta[i,3 ] * pow(age[j],2) + beta[i,4]* Y[i,j-1]
		 }                                  
	     gbv1[i]<- inprod(alpha1[], Z[i,])
	     gbv2[i]<- inprod(alpha2[], Z[i,])
	     gbv3[i]<- inprod(alpha3[], Z[i,])
		my[i,1]<-intercept1 + gbv1[i]
		my[i,2]<-intercept2 + gbv2[i]+roo21*beta[i,1]
		my[i,3]<-intercept3 + gbv3[i]+roo31*beta[i,1]+roo32*beta[i,2]
		my[i,4]<-0
		for(j in 1:3){
			beta[i,j]~dnorm(my[i,j], tau[j])
		}
		beta[i,4]<-0
	}


intercept1~dnorm(0,0.001)
intercept2~dnorm(0,0.001)
intercept3~dnorm(0,0.001)  

roo21~dnorm(0,0.001)
roo31~dnorm(0,0.001)
roo32~dnorm(0,0.001)

tau[1] ~dgamma(1.0E-3,1.0E-3)
sigma[1] <- 1/tau[1]
tau[2] ~dgamma(1.0E-3,1.0E-3)
sigma[2] <- 1/tau[2]
tau[3] ~dgamma(1.0E-3,1.0E-3)
sigma[3] <- 1/tau[3]
#tau[4] ~dgamma(1.0E-3,1.0E-3)
#sigma[4] <- 1/tau[4]

for(j in 1:loc){
   a1[j]~dnorm(0,prec1[j])
   lntau1[j] ~ dunif(-5,50)
   prec1[j]<-exp(lntau1[j])
   var1[j]<-1/prec1[j]
}

for(j in 1:loc){
   a2[j]~dnorm(0,prec2[j])
   lntau2[j] ~ dunif(-5,50)
   prec2[j]<-exp(lntau2[j])
   var2[j]<-1/prec2[j]
}

for(j in 1:loc){
   a3[j]~dnorm(0,prec3[j])
   lntau3[j] ~ dunif(-5,50)
   prec3[j]<-exp(lntau3[j])
   var3[j]<-1/prec3[j]
}


#lntautot~ dunif(-2,50)
#prectot<-exp(lntautot)
#sigma2tot<-1/prectot

#for(j in 1:years) {
#    tautot[j]~dgamma(1.0E-2,1.0E-2)
#    sigma2tot[j]<-1/tautot[j]
#   tau.block[j]~dgamma(1.0E-3,1.0E-3)
#   sigma2.block[j]<-1/tau.block[j]
#  }

for(k in 1:loc) {
	p1[k]<-1/loc
	x1[k]~dbern(p1[k])
	alpha1[k]<-x1[k]*a1[k]
	varalpha1[k]<-x1[k]*var1[k]
	}	
for(k in 1:loc) {
	p2[k]<-1/loc
	x2[k]~dbern(p2[k])
	alpha2[k]<-x2[k]*a2[k]
	varalpha2[k]<-x2[k]*var2[k]
	}	
for(k in 1:loc) {
	p3[k]<-1/loc
	x3[k]~dbern(p3[k])
	alpha3[k]<-x3[k]*a3[k]
	varalpha3[k]<-x3[k]*var3[k]
	}	
NQTL1<-sum(x1[1:loc])	
NQTL2<-sum(x2[1:loc])	
NQTL3<-sum(x3[1:loc])	

varb[1]<-sd(beta[1:ind,1])*sd(beta[1:ind,1])
varb[2]<-sd(beta[1:ind,2])*sd(beta[1:ind,2])
varb[3]<-sd(beta[1:ind,3])*sd(beta[1:ind,3])

vareff[1]<-sd(gbv1[1:ind])*sd(gbv1[1:ind])
vareff[2]<-sd(gbv2[1:ind])*sd(gbv2[1:ind])
vareff[3]<-sd(gbv3[1:ind])*sd(gbv3[1:ind])

her[1] <- (varb[1]-sigma[1])/varb[1]
her[2] <- (varb[2]-sigma[2])/varb[2] 
her[3] <- (varb[3]-sigma[3])/varb[3] 

for(i in 1:ind) {
	prod12[i]<-beta[i,1]*beta[i,2]
	prod13[i]<-beta[i,1]*beta[i,3]
	prod23[i]<-beta[i,2]*beta[i,3] 
	}
cov12<-(sum(prod12[1:ind])-sum(beta[1:ind,1])*sum(beta[1:ind,2])/ind)/(ind-1)
cov13<-(sum(prod13[1:ind])-sum(beta[1:ind,1])*sum(beta[1:ind,3])/ind)/(ind-1)
cov23<-(sum(prod23[1:ind])-sum(beta[1:ind,2])*sum(beta[1:ind,3])/ind)/(ind-1)

gcor12<-gcov12/(sd(beta[1:ind,1])*sd(beta[1:ind,2]))
gcor13<-gcov13/(sd(beta[1:ind,1])*sd(beta[1:ind,3]))
gcor23<-gcov23/(sd(beta[1:ind,2])*sd(beta[1:ind,3]))

gcov12 <- cov12-rcov12
gcov13 <- cov13-rcov13
gcov23 <- cov23-rcov23

rcov12<-roo21*sd(beta[1:ind,1]) 
rcov13<-roo31*sd(beta[1:ind,1])
rcov23<-roo32*sd(beta[1:ind,2]) 
}