Fast empirical Bayesian LASSO for multiple quantitative trait locus mapping

Linear model of multiple QTLs

Let y _i be the phenotypic value of a quantitative trait of the i th individual in a mapping population. Suppose we observe y _i , i = 1, ⋯, n of n individuals and collect them into a vector y = [y₁, y₂, ⋯, y _n ] ^T . Considering environmental effects, main and epistatic effects of all markers and gene-environment (G × E) interactions, we have the following linear regression model for y:

(1)

where μ is the population mean, vectors β _E and β _G represent the environmental effects and the main effects of all markers, respectively, vectors β _GG and β _GE capture the epistatic effects and the G × E interactions, respectively, X E, X G, X _GG and X _GE are the corresponding design matrices of different effects, and e is the residual error that follows a normal distribution with zero-mean and covariance . Throughout the paper we use I to denote an identity matrix whose size can be clearly identified from the context.

The design matrix X _G depends on a specific genetic model. We adopt the widely used Cockerham genetic model as also used by [18] in their generalized linear model for multiple QTL mapping. For a back-cross design, the Cockerham model defines the values of the main effect of a marker as -0.5 and 0.5 for two genotypes at the marker. For an intercross (F₂) design, there are two possible main effects named additive and dominance effects. The Cockerham model defines the values of the additive effect as -1, 0 and 1 for the three genotypes and the values of the dominance effect as -0.5 and 0.5 for homozygotes and heterozygotes, respectively. The columns of the design matrix X _GG are obtained as the element-wise product of any two different columns of X _G , and similarly the columns of X _GE are constructed as the element-wise product of any pair of columns from X _E and X _G .

Defining , X = [X _E , X _G , X _GG , X _GE ], we can write (1) in a more compact form:

(2)

Suppose that there are p environmental covariates and q markers whose main effects are additive, then the size of matrix X is n × k where k = p + q(q + 1)/2 + pq. Typically, we have k ≫ n. If dominance effects of the markers are considered, k is even larger. Our goal is to estimate all possible environmental and genetical effects on y manifested in the regression coefficients β, which is a challenging problem because k ≫ n. However, we would expect that most elements of β are zeros and thus we have a sparse linear model. Taking into account this sparsity, we will adopt the Bayesian LASSO model [10] where appropriate prior distributions are assigned to the elements of β as described in the next section.

Prior and posterior distributions

The unknown parameters in model (2) are, μ and β. While our main concern is β, parameters μ and need to be estimated so that we can infer β. To this end, we assign a noninformative uniform prior μ to and , i.e., p(μ) ∝ 1 and p( ) ∝ 1. Following the Bayesian LASSO model [10], we assume a three-level hierarchical model for β. Let us denote the elements of β as β _i , i = 1, 2, ⋯, k. At the first level, β _i , i = 1, 2, ⋯, k follow independent normal distributions with mean zero and unknown variance . At the second level, , follow independent exponential distribution with a common parameter . For a given λ, the distribution of β _i is found to be the Laplace distribution: , which is known to encourage the shrinkage of β _i toward zero [4]. However, the degree of shrinkage strongly depends on the value of λ. To alleviate the problem of choosing an inappropriate value for λ, we add another level to the hierarchical model at which we assign a conjugate Gamma prior Gamma(a, b) with a shape parameter a > 0 and an inverse scale parameter b > 0 to the parameter λ. As discussed in [10], we can pre-specify appropriate values for a and b so that the Gamma prior for λ is essentially noninformative.

Let us define vector . The joint posterior distribution of all the parameters (μ, , β, σ², λ) can be easily found [10]. In principle, MCMC simulation can be employed to draw samples from the posterior distribution for each parameter. However, since the number of parameters 2k + 3 in our model can be very large, the fully Bayesian approach based on MCMC sampling requires a prohibitive computational cost. To avoid this problem, Xu developed an empirical Bayes method for inferring β[12]. Our goal here is to develop a much faster and more accurate empirical Bayes method that can easily handle tens of thousands of variables.

Maximum a posteriori estimation of variance components

Similar to the EB method of [12], our EBLASSO first estimates σ², and μ, and then finds the posterior distribution of β based on the estimated parameters. Since λ is a parameter that we do not want to estimate, we can find the prior distribution of independent of λ as follows

(3)

The posterior distribution of μ, β, σ² and is given by

(4)

The marginal posterior distribution of μ, σ² and can then be written as

(5)

Let us define the precision of β _i as , i = 1, 2, ..., k and let α= [α₁, α₂, ⋯, α _k ]. It turns out to be more convenient to estimate α rather than σ² as will be shown shortly. Let us collect all parameters that need to be estimated as θ= (μ, σ₀, α). The log marginal posterior distribution of θ can be found from (5) as follows

(6)

where is the covariance matrix of y with a given α.

Similar to the EB method [12] and the RVM [15, 16], we will iteratively estimate each parameter by maximizing the log marginal posterior distribution L(θ) with the other parameters being fixed. Since it is not difficult to find the optimal μ and in each iteration as shown in [12, 15], we will give the expressions for μ and later but focus on the estimation of α now. Let us define and . Then following the derivations in [16], we can write L(θ) in (6) as L(θ) = L(θ_-i) + L(α _i ), where L(θ_-i) does not depend on α _i and L(α _i ) is given by

(7)

with and . If we assume a > -1.5 and b > 0, we prove in Additional file 1 that L(α _i ) has a unique global maximum and that the optimal α _i maximizing L(α _i ) is given by

(8)

where we have defined , , , and . Note that the Gamma distribution requires that a > 0 and b > 0 as we mentioned earlier. When -1.5 <a ≤ 0 as assumed here, we essentially use an improper distribution. In Additional file 1, we show that this improper distribution appears appropriate for getting a point estimation of α _i given in (8). It turns out that negative values of a give one more degree of freedom to adjust the shrinkage as will be demonstrated later in the Results section. Moreover, if a = -1, the last term in the right hand side of (7) disappears. In this case, we essentially use a noninformative uniform prior. Then it is not difficult to verify that (8) gives the optimal α _i derived in [16]:

(9)

Note that α* in (8) and (9) depends on other unknown parameters through s _i and q _i , and thus, α _i will be estimated iteratively as detailed in the EBLASSO algorithm described in the next section. Comparing with the EB method [12], our method finds each α _i (and equivalently ) in a closed form, whereas the EB method generally needs to employ a numerical optimization algorithm to find each . Therefore, our method not only saves much computation but also gives more accurate estimate of each α _i . Moreover, exploiting the similar techniques used in the RVM [16], we can further increase computational efficiency as described in the ensuing section.

Fast Empirical Bayesian LASSO algorithm

Note that when α _i = ∞, we have β _i = 0. Therefore, in each iteration, we can construct a reduced model of (2) that includes only nonzero β _i s and the corresponding columns of X. Let x _i be the i th column of X. If α _i = ∞ in the previous iteration but is finite in the current iteration, then we add x _i to the model and set ; if α _i is finite in the previous iteration but in the current iteration, we delete x _i from the model and set α _i = ∞; if both α _i and are finite, we retain x _i in the model and update α _i as . This can be done for all i s in a pre-specified order in each iteration. Alternatively, we can employ a greedy and potentially more efficient method to construct the model as described in the following. We define two iteration loops: an outer iteration loop and an inner iteration loop. In each outer iteration, we estimate μ and . In the inner iterations, assuming μ and are known and fixed, we estimate each α _i and construct the model. Specifically, in each inner iteration, we first calculate each from (8) and hen find , where α _i stands for the value of α _i obtained in the previous inner iteration. This step basically identifies the x _j that gives the largest increase in the log posterior distribution. Then we add, delete or retain x _j as described early. The inner iterations can run until a local convergence criterion is satisfied. Let vector contain all finite α _i s, vector consist of the corresponding β _i s and matrix contain the corresponding columns of X. Then we essentially construct the following reduced model:

(10)

where the number of columns of , k _r , is typically much less than the number of rows, n. This property will be used to reduce computation.

To calculate in (8), we need to first calculate s _i and q _i which requires . Since C_-iis different for different i, it may need large computation to calculate all . However, it was shown in [16] that we can calculate s _i and q _i as follows

(11)

where and . This requires only one matrix inversion in each iteration for calculating all s _i and q _i , i = 1, ⋯, k. However, since C is a relatively large matrix of size n × n, direct calculation of C^-1 may still require large computation. To avoid this problem, we can use the Woodbury matrix identity to derive an expression for C^-1:

(12)

where

(13)

with . The size of matrix ∑ is k _r × k _r which is typically much smaller than the size of C. Since inverting a matrix of N × N using an efficient method such as QR decomposition needs computation of O(N³), calculating ∑ requires much less computation than directly inverting C. Using (12), we can calculate S _i and Q _i as follows:

(14)

So far we have derived the method for efficiently estimating α. Other two unknown parameters μ and can be obtained by setting and . This gives

(15)

where 1 is a vector whose elements are all 1, and [15]

(16)

where ∑ _ii is the i th diagonal element of ∑ and

(17)

After we get estimates of μ, and α, we finalize the model (10), where the prior distribution of is now . For those x _i s not in the model, we can declare that they do not affect the quantitative trait because their regression coefficient is zero. For those x _i s in the matrix , the posterior distribution of their regression coefficients is Gaussian with covariance ∑ in (13) and mean u in (17) [15]. We can then use the z-score or more conservative t-statistics to test if at certain significance level. In this paper, the posterior mean u _j of the j th effect is reported as the empirical Bayes estimate of β _j , denoted by , and the posterior standard deviation, , is used as the standard error of .

We now summarize our fast EBLASSO algorithm as follows.

Algorithm 1 (EBLASSO algorithm)

1. Initialize parameters: choose a > -1.5 and b > 0, calculate μ = 1 ^T y/n, and set to be a small number, e.g., .

2. Initialize the model: Find , and calculate α _j from (9), set all other α _i s to be ∞ and .

3. Calculate ∑ from (13), S _i and Q _i , ∀i, from (14).

4. Update the model

   while the local convergence criterion is not satisfied

      Calculate q _i and s _i from (11), ∀i.

      Calculate from (8), ∀i.

Find .

   if

      if x _j is in the model, delete it and update ∑ , S _i and Q _i , ∀i.

   else

      if x _j is in the model, set and update ∑ , S _i and Q _i , ∀i.

      if x _j is not in the model, add it, set and update ∑ , S _i

      and Q _i , ∀i.

   end if

   end while

5. Update the residual variance using (16).

6. Calculate ∑ from (13) and C^-1from (12).

7. Update the fixed effects μ using (15) and update .

8. Calculate S _i and Q _i from (14).

9. If the global convergence criterion is not satisfied, go to step 4.

10. Find the posterior mean of from (17) and the posterior variance ∑ _ii from (13).

11. Use t-statistics to test if .

The parameters a and b can be set to be a small number (e.g., a = b = 0.01) so that the Gamma prior distribution is essentially noninformative [10]. Alternatively, we can use the predicted error (PE) obtained from cross-validation [4] to choose the values of a and b. As done in [16], the initial value of is chosen in step 1 to be a small number and the initial model is selected in step 2 to have a single effect that corresponds to the maximum L(α _i ) with a = -1. The outer iteration loop consists of steps 4-9, while the inner iteration loop is step 4, where we use the greedy method described earlier to update the model. In step 4, we use the method given in the Appendix of [16] to efficiently update ∑, S _i and Q _i after a x _j is added to or deleted from the model or after α _j is updated. The local convergence criterion can be defined as the simultaneous satisfaction of the following three conditions: 1) no effect can be added to or deleted from the model, 2) the change of L(θ) between two consecutive inner iterations is smaller than a pre-specified small value, and 3) the change of between two consecutive inner iterations is less than a pre-specified value. The global convergence criterion can be defined as the simultaneous satisfaction of the following two conditions: the change of L(θ) between two consecutive outer iterations is smaller than a pre-specified small value, and 2) the total change in μ, and between two consecutive outer iterations is less than a pre-specified value. A Matlab program has been developed to implement the algorithm; and a more efficient C++ program is under development.

Simulation study

We simulated a single large chromosome of 2400 centiMorgan (cM) long covered by evenly spaced q = 481 markers with a marker interval of 5 cM. The simulated population was an F₂ family derived from the cross of two inbred lines with a sample size n = 1000. The genotype indicator variable for individual i at marker k is defines as X _ik = 1, 0, -1 for the three genotypes, A₁A₁, A₁A₂ and A₂A₂, respectively. Twenty markers are QTLs with main effects and 20 out of the marker pairs have interaction effects. The locations and effects of the markers and maker pairs are shown in Table 1. Environmental effects and G × E effects were not simulated. The true population mean is μ = 100 and the residual variance is .

Markers a	Position	True b	EBLASSO c	EB c
( i, j )	(cM, cM)	β ( h 2 )
(11,11)	(50,50)	4.47(0.0975)	4.5801(0.1612)	4.8593(0.2075)
(26,26)	(125,125)	3.16(0.0524)	3.0768(0.1576)	3.3221(0.2035)
(42,42)	(205,205)	-2.24(0.0250)	-2.3169(0.1796)	-2.2769(0.2262)
(48,48)	(235,235)	-1.58(0.0128)	-1.3171(0.1720)	-1.3634(0.2205)
(72,72)	(355,355)	2.24(0.0247)	-	1.6537(0.4277)
(73,73)	(360,360)	3.16(0.0506)	5.1247(0.1555)	3.8771(0.4219)
(123,123)	(610,610)	1.10(0.0062)	-	1.5168(0.2432)
(127,127)	(630,630)	-1.10(0.0063)	-	-1.1834(0.2460)
(161,161)	(800,800)	0.77(0.0030)	-	-
(181,181)	(900,900)	1.73(0.0152)	-	-
(182,182)	(905,905)	3.81(0.0725)	5.6744(0.2400)	5.5127(0.2894)
(185,185)	(920,920)	2.25(0.0263)	1.7123(0.2327)	1.7070(0.2858)
(221,221)	(1100,1100)	-1.30(0.0088)	-1.4276(0.1506)	-1.0867(0.1956) d
(243,243)	(1210,1210)	-1.00(0.0051)	-0.8603(0.1486)	-
(262,262)	(1305,1305)	-2.24(0.0245)	-2.2539(0.1826)	-1.6078(0.2417)
(268,268)	(1335,1335)	1.58(0.0120)	2.4264(0.2040)	2.1736(0.2509)
(270,270)	(1345,1345)	1.00(0.0049)	-	-
(274,274)	(1365,1365)	-1.73(0.0147)	-1.4114(0.1800)	-1.4935(0.2254)
(361,361)	(1800,1800)	0.71(0.0026)	0.7856(0.1457)	0.6520(0.1859) d
(461,461)	(2300,2300)	0.89(0.0040)	-	-
(5,6)	(20,25)	2.24(0.0230)	1.7839(0.1654)	1.5752(0.2886)
(6,39)	(25,190)	2.25(0.0128)	1.9691(0.2168)	-
(42,220)	(205,1095)	4.47(0.0511)	4.3836(0.2198)	4.6414(0.3394)
(75,431)	(370,2150)	0.77(0.0014)	1.1360(0.2124) d	-
(81,200)	(400,995)	-2.24(0.0128)	-2.4190(0.2460)	-
(82,193)	(405,960)	1.58(0.0063)	1.6345(0.2442)	-
(87,164)	(430,815)	3.16(0.0235)	2.9263(0.2254)	1.7059(0.3319) d
(87,322)	(430,1605)	3.81(0.0342)	4.1019(0.2274)	3.7040(0.3632)
(92,395)	(455,1970)	1.73(0.0081)	1.5714(0.2065) d	-
(104,328)	(515,1635)	1.00(0.0024)	0.8081(0.1979) d	-
(118,278)	(585,1385)	-2.24(0.0120)	-2.0796(0.2221)	-2.2590(0.3460)
(150,269)	(745,1340)	1.10(0.0028)	1.0740(0.2142)	-
(237,313)	(1180,1560)	0.71(0.0014)	-	-
(246,470)	(1225,2345)	-1.10(0.0032)	-1.2381(0.2114) d	-
(323,464)	(1610,2315)	0.89(0.0020)	-	-
(328,404)	(1635,2015)	-1.73(0.0079)	-2.3036(0.2123)	-1.9428(0.3330)
(342,420)	(1705,2095)	-1.30(0.0041)	-1.3886(0.2121) d	-
(344,407)	(1715,2030)	-1.00(0.0025)	-	-
(373,400)	(1860,1995)	-1.58(0.0070)	-1.4732(0.2028)	-
(431,439)	(2150,2190)	3.16(0.0278)	2.6700(0.2121)	2.2454(0.3366) d
μ		100	100.70	100.59
		10	11.76	0.25
CPU time			3.4 mins	249 hrs

^a When i = j, the QTL is a main effect; otherwise, it is an epistatic effect.

^b The true value of a QTL effect is denoted by β and the proportion of variance contributed by the QTL is denoted by h².

^c The estimated QTL effect is denoted by and the standard error is denoted by . The EBLASSO algorithm used hyperparameters a = b = 0.1 and the EB algorithm used hyperparameters τ = -1 and ω = 0.001.

^d The estimated QTL effect was obtained from a neighboring marker (≤ 20 cM away) rather than from the maker with the true effect.

The total phenotypic variance for the trait can be written as

(18)

where cov(x _j , x_j') is the covariance between X _j and X_j'if j ≠ j' or the variance of X _j if j = j', which can be estimated from the data. The total phenotypic variance was calculated from (18) to be and the total genetic variance contributed by the main and epistatic effects of the markers was calculated from the second term of the right hand side of (18) to be 88.67. If we ignore the contributions from the covariance terms which are relatively small, the proportions of the phenotypic variance explained by a particular QTL effect j can be approximated by

(19)

where var(x _j ) is the variance of X _j . In the simulated data, the proportion of contribution from an individual QTL varied from 0.30% to 9.75%, whereas the proportion of contribution from a pair of QTLs ranged from 0.26% to 7.25%, as shown in Table 1. Some of the markers had only main or epistatic effect, while the others had both main and epistatic effects. The QTL model contained a total of possible effects, a number about 115 times of the sample size.

The data were analyzed in Matlab on a personal computer (PC) using the EBLASSO algorithm, the EB method, the RVM and the LASSO. The Matlab program SPARSEBAYES for the RVM was downloaded from http://www.miketipping.com. We translated the original SAS program for the EB method [12] into Matlab, and slightly modified the program to avoid possible memory overflow due to the large number of possible effects. We also got the program glmnet [19] that is a very efficient program implementing the LASSO and other related algorithms. The PC version of glmnet uses Matlab to initialize and call the core LASSO algorithm that is implemented efficiently with Fortran code.

We used the PE [4] obtained from ten-fold cross validation to select the values of hyperparameters a and b in our EBLASSO algorithm. Ideally, we should test a large set of values for a ≥ -1 and b > 0, but this may be time consuming. Therefore, we first ran cross-validation for the following values: a = b = 0.001, 0.01, 0.05, 0.1, 0.5, 1; the degree of shrinkage generally decreases along this path. Table 2 lists the PEs and the standard errors of the PEs for different values of a and b. It is seen that the PE for a = b = 0.1 is the smallest, although it is close to the PEs for a = b = 0.05 and 0.5 but relatively smaller than PEs for the other values of a and b. To see if a = b = 0.1 is the best set of values, we further ran cross-validation with b = 0.1 and a = 0.5 or -0.01. For a fixed b, the degree of shrinkage decreases when a decreases. It is seen from Table 2 that the PE for b = 0.1 and a = 0.5 or -0.01 is greater than that for a = b = 0.1. Therefore, cross-validation gave a = b = 0.1 as the best set of values. Table 2 also lists the number of effects detected with different values of a and b. All 30 effects detected with a = b = 0.1 are presented in Table 1 and shown in Figures 1 and 2.

Algorithm	Parameters◇	PE ± STE*	Number of effects†‡		CPU time (mins)
	(0.001, 0.001)	16.49 ± 0.8908	25/0	13.5	3.4
	(0.01, 0.01)	15.95 ± 0.7477	28/0	12.46	3.4
	(0.05, 0.05)	15.89 ± 0.7498	30/0	11.72	3.4
	(0.1, 0.1)	15.81 ± 0.8359	30/0	11.72	3.4
EBLASSO	(0.5, 0.5)	15.86 ± 0.7717	31/0	11.57	3.4
	(1, 1)	16.07 ± 0.7203	29/0	12.31	3.4
	(0.5, 0.1)	16.14 ± 0.8557	28/0	12.5	3.4
	(-0.01, 0.1)	15.92 ± 1.0161	32/1	11.31	3.4
	(-1, 0.0001)	-	14/1	21.22	2,760.0
EB	(-1, 0.0005)	-	13/1	12.15	4,140.0
	(-1, 0.001)	-	22/1	0.25	14,940.0
	(-1, 0.01)	-	8/0	0.01	2,760.0

^◇ Parameters are (a, b) for the EBLASSO, and (τ, ω) for the EB.

*The average PE and the standard error were obtained from ten-fold cross validation.

^†The number of effects and residual variance were obtained using all 1000 samples not from cross validation.

^‡The first number is the number of true positive effects; the second number is number of false positive effects. All the effects counted have a p-value ≤ 0.05.

To test the performance of the EB method, we ran the program with the following parameter values: τ = -1, ω = 0.0001, 0.0005, 0.001 and 0.01, which yielded 14, 13, 22, and 8 true simulated effects, respectively, and 1, 1, 1, and 0 false effects, respectively, as shown in Table 2. Cross-validation was not done because it was too time-consuming, and thus, the optimal values for the parameters could not be determined. Nevertheless, we listed 22 true positive effects estimated with τ = -1 and ω = 0.001 in Table 1, which reflects the best performance of the EB method with the set of parameters values tested. We also plotted these 22 effects in Figures 3 and 4. Comparing the effects detected by EBLASSO and EB methods shown in Table 1 and in Figures 1, 2, 3 and 4, the EBLASSO method detected 13 (17) true main (epistatic) effects, whereas the EB method detected 15 (7) true main (epistatic) effects. Overall, the EBLASSO detected 8 more true effects than the EB method without any false positive effects, whereas the EB method gave one false positive effect. We would like to emphasize here that the EB method may detect less number of true effects in practice because as we mentioned earlier it is too time-consuming to choose the optimal values for the parameters τ and ω. To see if EBLASSO could estimate QTL effects robustly, we simulated three replicates of the data: each replicate consists of 1000 individuals whose genotypes at 481 markers were independently generated and whose phenotypes were calculated from (2) with e independently generated from Gaussian random variables with zero mean and covariance 10I. We performed cross-validation and determined the best values of a and b for each replicate. Using these values, we ran the EBLASSO and identified 35, 38, 34 true positive effects and 4, 6, 2 false positive effects, respectively, for three replicates. These results showed that the EBLASSO could detect QTL effects robustly.

The EBLASSO took about 3.4 minutes of CPU time for each set of values of a and b listed in Table 2, whereas the EB took 249 hours of CPU time for τ = -1 and ω = 0.001 and about 46, 69, 46 hours for τ = -1 and ω = 0.0001, 0.0005, 0.01, respectively. This simulation study showed that the EBLASSO method not only can detect more effects, but also offers a huge advantage in terms of computational time. Note that all simulations were done in Matlab. It is expected that the EBLASSO algorithm will be even faster, after its implementation in C++ is completed.

We wished to test the performance of the RVM and the LASSO on the simulated data. To this end, we replaced the inner iteration in our EBLASSO algorithm with the program SPARSEBAYES that implemented the RVM. Although we carefully modified SPARSEBAYES to avoid possible memory overflow due to the large number of possible effects, the program ran out of memory after one or two outer iterations. Hence, we did not get any results from the RVM method. Considering the QTL model in (2), the LASSO tries to estimate μ and the QTL effects β as follows

(20)

where λ is a positive constant that can be determined with cross-validation [4]. We tried to run the program glmnet [19] with the simulated data. However, glmnet could not handle the big design matrix X of 1, 000 × 115, 922 in our QTL model, and we did not get any results from glmnet.

The optimal values for the parameters of the EB method were τ = -1 and ω = 0.01, since they gave the smallest PE in cross-validation as listed in Table 3. With τ = -1 and ω = 0.01, the EB method detected 19 true effects and 4 false positive effects. The RVM detected all 20 true effects as the EBLASSO did, but it also output a large number of 42 false positive effects. This result is consistent the observation [13] that the uniform prior distribution used in the RVM usually yields many false positive effects. To choose the optimal value of λ for the LASSO, we ran ten-fold cross validation starting from λ = 4.9725 (which gave only one nonzero effect) and then decreasing λ to 0.0025 with a step size of 0.0768 on the logarithmic scale (Δ ln(λ) = 0.0768). The smallest PE was achieved at λ = 0.0675. We then used this value to run glmnet on the whole data set, which yielded 97 nonzero effects. For each of these nonzero effects, we calculated their standard error using equation (7) in [4], and then calculated the p-value of each nonzero effect. This gave 20 true effects and 48 false positive effects with a p-value less than 0.05. Comparing the number of effects detected by the EBLASSO, EB, RVM and LASSO, the EBLASSO offered the best performance because it detected all true effects and a very small number of false positive effects.

It is seen from Table 3 that the EBLASSO and the LASSO took much less time than the EB method and the RVM on analyzing this data set. It is expected that the EBLASSO is much faster than the EB method because as we discussed earlier the EB needs a numerical optimization procedure. The RVM and EBLASSO generally should have similar speed because two algorithms use the similar technique to estimate effects. However, when applying to the same data set, the RVM often yields a model with much more nonzero effects than the EBLASSO as is the case here, because the RVM does not provide sufficient degree of shrinkage. Due to this reason, the RVM algorithm requires more time than the EBLASSO. The LASSO took slightly less CPU time than the EBLASSO in this example. However, we would emphasize that the LASSO was implemented with Fortran but our EBLASSO was implemented with Matlab. The speed of EBLASSO is expected to increase significantly once it is implemented in C/C++.

Real data analysis

This dataset was obtained from [20]. This dataset consists of n = 150 double haploids (DH) derived from the cross of two spring barley varieties Morex and Steptoe. The total number of markers was q = 495 distributed along seven pairs of chromosomes of the barley genome, covering 206 cM of the barley genome. The phenotype was the spot blotch resistance measured as the lesion size on the leaves of barley seedlings. Note that spot blotch is a fungus named Cochliobolus sativus. This dataset was used as an example for the application of the EBLASSO method. Genotype of the markers were encoded as +1 for genotype A (the Morex parent), -1 for genotype B (the Steptoe parent), and 0 for missing genotype. Ideally, the missing genotypes should be imputed from known genotypes of neighboring markers. For simplicity, we replaced the missing genotypes with 0 in order to use the phenotypes of the individuals with missing genotypes. The total missing genotypes only account for about 4.2% of all the genotypes. Including the population mean, the main and the pair-wise epistatic effects, the total number of model effects was , about 818 times as large as the sample size.

Table 4 gives the average PE and the standard error obtained from 5-fold cross validation, the residual variance and the number of effects detected by the EBLASSO method for different values of a and b. It is seen that the PEs for a = b = 0.001, 0.01.0.05 are almost the same but are smaller than the PEs for other larger a and b. However, when a = b = 0.001 or 0.01, only one or two effects were detected. When a = b = 0.1, 0.5 or 1, the residual variance is very small, implying that the model is likely over-fitted. Specifically, the number of columns of matrix in the model (10) is equal to the total number nonzero effects, which is more than 120 for a = b = 0.1, 0.5 or 1 as indicated in Table 4. Hence, since the number of samples (150) is relative small, y - μ can be almost completely in the column space of , which results in very small residual variance. Based on these observations, it seems that a = b = 0.05 gives reasonable results, because the PE is among the smallest and the residual variance is relatively but not unreasonably small. Nevertheless, in order to estimate effects more reliably, we searched over all effects detected with a = b = 0.05, 0.07, 0.1, 0.5 or 1, and found eight effects are detected with all these values. Markers or marker pairs of these eight effects and their values estimated with a = b = 0.05 were listed in Table 5. All 11 effects estimated with a = b = 0.05 were also plotted in Figure 5. As seen from Table 5, one single QTL at marker 446 contributes most of the phenotypic variance (about 76%), while the other QTL effects contribute from about 0.8% to 2.6% of the variance. We also analyzed the data with the EB method. One effect at marker 446 was detected with a p-value < 0.05 when τ = -1, ω = 0.0001 or 0.0005, no effect was detected with p-value ≤ 0.05 when τ = -1, ω = 0.001 or 0.01. The CPU time of the EB was about 46 minutes for each set of values of τ and ω tested; whereas the CPU time of the EBLASSO method was about 3 minutes for a = b = 0.05.

a= b	PE ± STE*	Number of effects†‡
0.001	0.70 ± 0.21	1/1/1	0.6706
0.01	0.79 ± 0.31	2/2/2	0.5996
0.05	0.70 ± 0.21	11/11/11	0.2699
0.07	0.96 ± 0.30	10/15/15	0.2104
0.1	1.20 ± 0.18	13/128/132	2.59E-06
0.5	1.21 ± 0.09	9/112/122	2.59E-06
1	1.25 ± 0.17	8/115/132	2.59E-06

*The average PE and the standard error were obtained from five-fold cross validation.

^†The number of effects and residual variance were obtained using all 150 samples not from cross validation.

^‡The first number is the number of effects with a p-value ≤ 0.05 and proportion of variance ≥ 0.5%; the second number is the number of effects with a p-value ≤ 0.05; the third number is the total number of non-zero effects reported by the program.

Markers
(446,446)	1.4173(0.0432)	0.7639
(187,187)	0.2624(0.0421)	0.0262
(77,77)	0.1881(0.0413)	0.0132
(238,238)	-0.1742(0.0427)	0.0101
(197,483)	0.1748(0.0405)	0.0117
(37,130)	0.1557(0.0422)	0.0085
(53,270)	-0.1649(0.0452)	0.0081
(149,175)	0.1697(0.0464)	0.0078
	5.7037
	0.2699
		0.8495

Eights effects were detected with all of the following parameters: a = b = 0.05, 0.07, 0.1, 0.5, 1; , and listed here were obtained with a = b = 0.05.