model {
  # offspring
  for (i in 1:273) {
     pheno[i] ~ dnorm(mu[i], tau.within)
     mu[i] <- nu + alpha[sire[i]] + beta[dam[i]]
  }

  # sires
  for (i in 1:9) {
     alpha[i] ~ dnorm(0.0, tau.sire)
  }
  alpha[10] <- 0

  # dams
  # sire 1
  for (i in 1:3) {
     beta[i] ~ dnorm(0.0, tau.dam)
  }
  beta[4] <- 0
  # sire 2
  for (i in 5:9) {
     beta[i] ~ dnorm(0.0, tau.dam)
  }
  beta[10] <- 0
  # sire 3
  for (i in 11:12) {
     beta[i] ~ dnorm(0.0, tau.dam)
  }
  beta[13] <- 0
  # sire 4
  for (i in 14:18) {
     beta[i] ~ dnorm(0.0, tau.dam)
  }
  beta[19] <- 0
  # sire 5
  beta[20] ~ dnorm(0.0, tau.dam)
  beta[21] <- 0
  # sire 6
  for (i in 22:25) {
     beta[i] ~ dnorm(0.0, tau.dam)
  }
  beta[26] <- 0
  # sire 7
  for (i in 27:29) {
     beta[i] ~ dnorm(0.0, tau.dam)
  }
  beta[30] <- 0
  # sire 8
  for (i in 31:33) {
     beta[i] ~ dnorm(0.0, tau.dam)
  }
  beta[34] <- 0
  # sire 9
  for (i in 35:40) {
     beta[i] ~ dnorm(0.0, tau.dam)
  }
  beta[41] <- 0
  # sire 10
  for (i in 42:44) {
     beta[i] ~ dnorm(0.0, tau.dam)
  }
  beta[45] <- 0

  # relationship between phenotype and fitness
  for (i in 1:273) {
    fit[i] ~ dnorm(mu.fit[i], tau.reg)
    mu.fit[i] <- b0 + b1*pheno[i]
  }
  b0 ~ dnorm(0.0, 0.001)
  b1 ~ dnorm(0.0, 0.001)
  w.bar <- nu*b1 + b0

  # observational components
  v.sire   <- sd.sire*sd.sire
  v.dam    <- sd.dam*sd.dam
  v.within <- sd.within*sd.within

  # priors
  nu <- mean(pheno[])
  sd.sire ~ dunif(0.01, sd.dam)
  sd.dam ~ dunif(sd.sire, min.sd.within)
  min.v.within.try <- 3*v.dam - v.sire
  min.v.within <- max(0.01, min.v.within.try)
  min.sd.within <- sqrt(min.v.within)
  sd.within ~ dunif(min.sd.within, 250)

  # precisions
  tau.sire <- 1.0/v.sire
  tau.dam <- 1.0/v.dam
  tau.within <- 1.0/v.within

  # priors for regression
  sd.reg ~ dunif(0, 10)
  v.reg <- sd.reg*sd.reg
  tau.reg <- 1.0/v.reg
}

# initial values to be loaded after compiling and before hitting
# "gen inits"
list(sd.sire=1, sd.dam=2, sd.within=10)