Bayesian inference of mixed models in quantitative genetics of crop species

Abstract

The objectives of this study were to implement a Bayesian framework for mixed models analysis in crop species breeding and to exploit alternatives for informative prior elicitation. Bayesian inference for genetic evaluation in annual crop breeding was illustrated with the first two half-sib selection cycles in a popcorn population. The Bayesian framework was based on the Just Another Gibbs Sampler software and the R2jags package. For the first cycle, a non-informative prior for the inverse of the variance components and an informative prior based on meta-analysis were used. For the second cycle, a non-informative prior and an informative prior defined as the posterior from the non-informative and informative analyses of the first cycle were used. Regarding the first cycle, the use of an informative prior from the meta-analysis provided clearly distinct results relative to the analysis with a non-informative prior only for the grain yield. Regarding the second cycle, the results for the expansion volume and grain yield showed differences among the three analyses. The differences between the non-informative and informative prior analyses were restricted to variance components and heritability. The correlations between the predicted breeding values from these analyses were almost perfect.

Appendix

(1) Code for JAGS

$$\begin{gathered} {\text{model }} \hfill \\ {\text{\{ }} \hfill \\ {\bf{\# likelihood function }} \hfill \\ \, \underbrace {\text{Y}} \sim {\text{dmnorm(mu[1:N,1],I[1:N,1:N]*tau}}\_{\text{e) }}_{{}} \hfill \\ {\text{ y}}|\beta , {\text{u}}_{ 1} , {\text{ u}}_{ 2} , { }\sigma_{{{\text{u}}}_{ 1} }^{ 2} , { }\sigma_{{{\text{u}}_{ 2} }}^{ 2} , { }\sigma_{\text{e}}^{ 2} \sim {\text{N}}(\mu ,{\text{I}}\sigma_{\text{e}}^{ 2} ) {\text{ Eq}} . { (1} . 1 )\hfill \\\# \underbrace {\text{Y}} \sim {\text{dt(mu[1:N,1],I[1:N,1:N]*tau}}\_{\text{e,v)}}_{{}} \hfill \\ {\text{ y}}|\beta , {\text{u}}_{ 1} , {\text{ u}}_{ 2} , { }\sigma_{{{\text{u}}_{ 1} }}^{ 2} , { }\sigma_{{{\text{u}}_{ 2} }}^{ 2} , { }\sigma_{\text{e}}^{ 2} , {\text{v}}\sim {\text{Student - t}}(\mu ,{\text{I}}\sigma_{\text{e}}^{ 2} , {\text{v) (assuming the Student - t distribution for the data)}} \hfill \\\# \underbrace {\text{Y}} \sim {\text{dmnorm(mu[1:N,1],I[1:N,1:N]*(1 - h}}^{ 2} ) * {\text{tau}}\_{\text{p) }}_{{}} \hfill \\ {\text{ y}}|\beta ,{\text{u}}_{ 1} , {\text{ u}}_{ 2} , {{\upsigma}}_{{{\text{u}}_{ 1} }}^{ 2} , {{\upsigma}}_{\text{p}}^{ 2} , {\text{ h}}_{{}}^{ 2} \sim{\text{N}}(\mu ,{\text{I(1 - h}}^{ 2} ) {{\upsigma}}_{\text{p}}^{ 2} ) {\text{ (assuming a reparameterization in terms of h}}^{ 2} {\text{ and}}\; \sigma _{\text{p}}^{ 2} )\hfill \\ \hfill \\ \, \underbrace {{{\text{mu[1:N,1] < - X[1:N,1:nbeta]\% *\% beta[1:nbeta,1] + Z1[1:N,1:nu1]\% *\% u1[1:nu1,1] + Z2[1:N,1:nu2]\% *\% u2[1:nu2,1]}}}}_\mu ={\text{X}}\beta + {\text{Z}}_{ 1} {\text{u}}_{ 1}+{\text{Z}}_{ 2} {\text{u}}_{ 2} \begin{array}{*{20}c} { \, } \\ \end{array}\hfill \\ \hfill \\ \end{gathered}$$

$$ \begin{gathered} \underbrace {{\text{tau}}}\_{\text{e}} \sim {\text{dgamma(ve/2, Se/2) }}\hfill \\ \, {{\upsigma}}_{\text{e}}^{ 2} | {{\upnu}}_{\text{e}} , {\text{S}}_{\text{e}} {{\sim \nu }}_{\text{e}} {\text{S}}_{\text{e}} {{\upchi}}_{{{{\upnu}}_{\text{e}} }}^{ - 2} {\text{ Eq}} . { (1} . 7 ) { } \hfill \\ \, \# {\text{ sd}}\_{\text{p}} \sim {\text{dunif (a, b)}} \hfill \\ \, {\underbrace {\#{\text{tau}}}}\_{\text{p}} < - {\text{sd}}\_p*{\text{sd}}\_{\text{p}} \hfill \\ {{\upsigma}}_{\text{p}} | {\text{a,b}} \sim {\text{U[a,b] (assuming a uniform prior for the phenotypic standard deviation)}} \hfill \\ \underbrace \# {\text{ h2}} \sim{\text{ dbeta(c,d)}} \hfill \\ {\text{ h}}^{ 2} | {\text{c,d}} \sim {\text{Beta[c,d] (assuming a Beta prior for h}}^{ 2} )\hfill \\ \underbrace \# {\text{ v}} \sim {\text{dunif(e,f)}} \hfill \\ {\text{ v}}|{\text{e,f }}\sim {\text{U[e,f] (assuming a uniform prior for degrees of freedom when using a Student - t for the data)}} \hfill \\ \hfill \\ {\bf {\# definition of variance components}} \hfill \\ {\text{ sigma2}}\_{\text{e}} < - 1/{\text{tau}}\_{\text{e }} \hfill \\ {\text{ sigma2}}\_{\text{u1}} < - 1/{\text{tau}}\_{\text{u1}} \hfill \\ {\text{ sigma2}}\_{\text{u2}} < - 1/{\text{tau}}\_{\text{u2}} \hfill \\ {\text{ sigma2}}\_{\text{a}} < - 4*{\text{sigma2}}\_{\text{u1}} \hfill \\ {\text{ h2}} < - {\text{sigma2}}\_{\text{u1}} / {\text{(sigma2}}\_{\text{u1}} + {\text{sigma2}}\_{\text{e/2 )}} \hfill \\ {\bf{\# assuming a reparametrization in terms of h}}^{ 2} {\text{ and}} \;\sigma _{\text{p}}^{ 2} \hfill \\ \#{\text{sigma2}}\_{\text{p}} < - 1/{\text{tau}}\_{\text{p }} \hfill \\ \#{\text{sigma2}}\_{\text{u1}} < -{\text{ sigma2}}\_{\text{p}}*{\text{h2 }} \hfill \\ {\text{\} }} \hfill \\ \end{gathered} $$

$$\begin{gathered} {\bf{\# prior distributions for the inverse of variance components}} \hfill \\ \, \underbrace {\text{ tau}}\_{\text{u1}} \sim {\text{dgamma(vu1/2, Su1/2)}} \hfill \\ {{\upsigma}}_{\text{u1}}^{ 2} | {{\upnu}}_{\text{u1}} , {\text{S}}_{\text{u1}} {{\sim \nu }}_{\text{u1}} {\text{S}}_{\text{u1}} {{\upchi}}_{{{{\upnu}}_{\text{u1}} }}^{ - 2} {\text{ Eq}} . { (1} . 5 )\hfill \\ \hfill \\ \, \underbrace {\text{tau}}\_{\text{u2}} \sim {\text{dgamma(vu2/2, Su2/2)}} \hfill \\ {{\upsigma}}_{\text{u2}}^{ 2} | {{\upnu}}_{\text{u2}} , {\text{S}}_{\text{u2}} {{\sim \nu }}_{\text{u2}} {\text{S}}_{\text{u2}} {{\upchi}}_{{{{\upnu}}_{\text{u2}} }}^{ - 2} {\text{ Eq}} . { (1} . 6 )\hfill \\ \end{gathered}$$

$$ \begin{gathered} {\begin{gathered} \, \underbrace {\text{tau}}\_{\text{e}} \sim {\text{dgamma(ve/2, Se/2) }}\hfill \\ {{\upsigma}}_{\text{e}}^{ 2} | {{\upnu}}_{\text{e}} , {\text{S}}_{\text{e}} {{\sim \nu }}_{\text{e}} {\text{S}}_{\text{e}} {{\upchi}}_{{{{\upnu}}_{\text{e}} }}^{ - 2} {\text{ Eq}} . { (1} . 7 ) { } \hfill \\ \# {\text{ sd}} \_{\text{p}} \sim {\text{dunif (a, b)}} \hfill \\ \underbrace {\#{\text{tau}}} \_{\text{p}} < - {\text{sd}} \_p*{\text{sd}} \_{\text{p}} \hfill \\ \end{gathered} } \hfill \\ {{\upsigma}}_{\text{p}} | {\text{a,b}} \sim {\text{U[a,b] (assuming a uniform prior for the phenotypic standard deviation)}} \hfill \\ \underbrace \# {\text{ h2}} \sim{\text{ dbeta(c,d)}} \hfill \\ {\text{ h}}^{ 2} | {\text{c,d}} \sim {\text{Beta[c,d] (assuming a Beta prior for h}}^{ 2} )\hfill \\ \underbrace \# {\text{ v}} \sim {\text{dunif(e,f)}} \hfill \\ {\text{ v}}|{\text{e,f }}\sim {\text{U[e,f] (assuming a uniform prior for degrees of freedom when using a Student - t for the data)}} \hfill \\ \hfill \\ {\bf {\# definition of variance components}} \hfill \\ {\text{ sigma2}}\_{\text{e}} < - 1/{\text{tau}}\_{\text{e }} \hfill \\ {\text{ sigma2}}\_{\text{u1}} < - 1/{\text{tau}}\_{\text{u1}} \hfill \\ {\text{ sigma2}}\_{\text{u2}} < - 1/{\text{tau}}\_{\text{u2}} \hfill \\ {\text{ sigma2}}\_{\text{a}} < - 4*{\text{sigma2}}\_{\text{u1}} \hfill \\ {\text{ h2}} < - {\text{sigma2}}\_{\text{u1}} / {\text{(sigma2}}\_{\text{u1}} + {\text{sigma2}}\_{\text{e/2 )}} \hfill \\ {\bf{\# assuming a reparametrization in terms of h}}^{ 2} {\text{ and}} \;\sigma _{\text{p}}^{ 2} \hfill \\ \#{\text{sigma2}}\_{\text{p}} < - 1/{\text{tau}}\_{\text{p }} \hfill \\ \#{\text{sigma2}}\_{\text{u1}} < -{\text{ sigma2}}\_{\text{p}}*{\text{h2 }} \hfill \\ {\text{\} }} \hfill \\ \end{gathered} $$

where Y is the phenotypic values vector; X, Z1, and Z2 are, respectively, incidence matrices for β, u ₁, and u ₂; N, nbeta, nu1, and nu2 are, respectively, the numbers of observations, of fixed effects, of families and of blocks; mean_beta, mean_u1, and mean_u2 are, respectively, the mean vectors of prior distributions for β, u ₁, and u ₂; I, Ibeta, A, Ib are, respectively, matrices related with covariance of prior distributions for e, β, u ₁, and u ₂; and ν. and S. are the hyperparameters for the inverse of the variance components.

(2) Code for R2jags

Y = as.matrix(read.table(“Yp1.txt”)) #reading phenotypic observations (Y)

X = as.matrix(read.table(“X.txt”)) #reading incidence matrix of β

Z1 = as.matrix(read.table(“Z.txt”)) #reading incidence matrix of u ₁

Z2 = as.matrix(read.table(“Jp.txt”)) #reading incidence matrix of u ₂

#specifying dimensions

N = nrow(Y) # number of observations in Y

nbeta = ncol(X) # number of fixed effects (β)

nu1 = ncol(Z1) # number of families (u ₁)

nu2 = ncol(Z2) # number of blocks (u ₂)

#mean vectors of prior distributions for location parameters

mean_beta = matrix(100,nbeta,1) # μ _β in Eq. 1.2

mean_u1 = matrix(0,nu1,1) # 0 in Eq. 1.3

mean_u2 = matrix(0,nu2,1) # 0 in Eq. 1.4

#matrices related with covariance of prior distributions for location parameters

I = diag(N) # I in Eq. 1.1

Ibeta = diag(nbeta) # $I_{{{\upbeta}}}$in Eq. 1.2

A = as.matrix(read.table(“A.txt”)) # A in Eq. 1.3

Ib = diag(nu2) # I _b in Eq. 1.4

#specifying hyperparameters for the inverse of variance components (non-informative prior)

v1 = 0.001; v2 = 0.001; ve = 0.001; S1 = v1*1; S2 = v2*1; Se = ve*1;

library(R2jags) #loading R2jags package

#listing JAGS input

jags.data = list(“Y”,“X”,“Z1”,“Z2”,“N”,“nbeta”,“nu1”, “nu2”,“mean_beta”, “mean_u1”,“mean_u2”, “Ibeta”, “A”, “Ib”, “I”,“v1”,“v2”,“ve”,“S1”,“S2”,“Se”)

#listing JAGS output

jags.params = c(“beta”,“u1”,“u2”,“sigma2_u1”,“sigma2_a”,“sigma2_u2”, “sigma2_e”,“h2”)

#listing initial values for MCMC simulation

jags.inits = function() {

list(“beta” = structure(.Data = c(4500,100,300),.Dim = c(nbeta, 1)), “u1” = structure(.Data = mean_u1,.Dim = c(nu1, 1)),

“u2” = structure(.Data = mean_u2,.Dim = c(nu2, 1)), “tau_u1” = c(0.0001), “tau_u2” = c(0.001), “tau_e” = c(0.00001))}

#calling jags function of R2jags package

bayes = jags(data = jags.data, jags.params, inits = jags.inits, n.chains = 1, n.iter = 70000, n.burnin = 20000, n.thin = 5, model.file = ”bayes_model.txt”) # “bayes_model.txt” is txt file containing model specified in Code for JAGS

#saving MCMC output

write.table(as.mcmc(bayes), “prod_noninf.txt”,row.names = FALSE,quote = FALSE)