Inferring genetic values for quantitative traits non-parametrically

Summary

Inferences about genetic values and prediction of phenotypes for a quantitative trait in the presence of complex forms of gene action, issues of importance in animal and plant breeding, and in evolutionary quantitative genetics, are discussed. Current methods for dealing with epistatic variability via variance component models are reviewed. Problems posed by cryptic, non-linear, forms of epistasis are identified and discussed. Alternative statistical procedures are suggested. Non-parametric definitions of additive effects (breeding values), with and without employing molecular information, are proposed, and it is shown how these can be inferred using reproducing kernel Hilbert spaces regression. Two stylized examples are presented to demonstrate the methods numerically. The first example falls in the domain of the infinitesimal model of quantitative genetics, with additive and dominance effects inferred both parametrically and non-parametrically. The second example tackles a non-linear genetic system with two loci, and the predictive ability of several models is evaluated.

1. Introduction

The problem considered here is that of inferring genetic values and of predicting phenotypes for a quantitative trait under complex forms of gene action, an issue of importance in animal and plant breeding, and in evolutionary quantitative genetics (Lynch & Walsh, 1998). The discussion is streamlined as follows. Current methods for dealing with epistatic variability via variance component models are discussed in section 2. Problems posed by cryptic, non-linear, forms of epistasis are identified in section 3. Section 4 proposes non-parametric definitions of additive effects (breeding values), with and without employing molecular information, and shows how these could be inferred using reproducing kernel Hilbert spaces (RKHS) regression models. Sections 5 and 6 present stylized examples to demonstrate the methods. The first example uses the infinitesimal model of quantitative genetics, and the second one tackles a non-linear genetic system. The paper ends with concluding comments.

2. Extant theory

A standard decomposition of phenotypic value in quantitative genetics (Falconer & Mackay, 1996) is

where a, d and i are additive, dominance and epistatic effects, respectively, and e is a residual, reflecting environmental (residual) variability. This linear decomposition can also be used to describe variability of latent variables, especially if assumed Gaussian (e.g. Dempster & Lerner, 1950; Gianola, 1982). The i effect can be decomposed into additive×additive, additive×dominance, dominance×dominance, etc., deviates. In what has been termed ‘statistical epistasis’ (Cheverud & Routman, 1995), these deviates are assumed to be random draws from some distributions representing ‘interactions’ between loci. Under certain assumptions (Cockerham, 1954; Kempthorne, 1954), the deviates are uncorrelated, leading to the standard variance decomposition σ²=σ_a²+σ_d²+σ_aa²+σ_ad²+σ_dd²+…+σ_e², where the σ²s are variance components, e.g., σ_ad² represents the contribution of additive by dominance effects to variance. The sum of all variance components other than σ_a², σ_d² and σ_e² is interpreted as variance ‘due to epistasis’. For the purposes of discussing problems this representation has, it suffices to assume that the only relevant epistatic effects are those of an additive×additive, additive×dominance and dominance×dominance nature.

(i) Henderson's methods for predicting epistatic effects

Using the Cockerham–Kempthorne (hereinafter, CK) assumptions, Henderson (1985) extended their model to the infinitesimal domain, and vectorially, by writing

(1)

where β is some nuisance location vector (equal to μ if it contains a single element); X is a known incidence matrix; a and d are vectors of additive and dominance effects, respectively; i_aa, i_ad and i_dd are epistatic effects, and g=a+d+i_aa+i_ad+i_dd is the ‘total’ genetic value. Assuming that g and e are uncorrelated, the variance–covariance decomposition is

(2)

where V_y,V_g and V_e are the phenotypic, genetic and residual variance–covariance matrices, respectively. Further,

(3)

Here, A is the numerator relationship matrix; D is a matrix due to dominance relationships which can be computed from entries in A, and the remaining matrices involve Hadamard (element by element) products of matrices A or D. Thus, under CK, all matrices can be computed from elements of A, as noted by Henderson (1985), because absence of inbreeding and of linkage disequilibrium are assumed. Cockerham (1956) and Schnell (1963) gave formulae for covariances between relatives in the presence of linkage, and difficulties posed by linkage disequilibrium are discussed by Gallais (1974).

In CK–Henderson, with the extra assumption that all genetic effects (having null means) and the data (with mean vector Xβ) follow a multivariate normal distribution with known dispersion components, one has

(4)

The best predictor (Henderson, 1973; Bulmer, 1980) of additive genetic value in the mean-squared error sense is

(5)

where E _g|y (g)=V_gV_y⁻¹ (y−Xβ) is the best predictor of g, and the matrix of regression coefficients V_gV_y⁻¹ is the multidimensional counterpart of heritability in the broad sense. With known variance components, β is typically estimated by generalized least-squares (equivalently, by maximum likelihood under normality) as . Then, ‘the empirical best predictor’ of the vector of additive effects is taken to be

(6)

This is precisely the best linear unbiased predictor (BLUP) of additive merit when model (1) holds. Likewise, the BLUP of additive×dominance deviations is

and so on.

The BLUPs of a, d, i_aa, i_ad and i_dd in (1) can be computed simultaneously using Henderson's mixed model equations, but this requires forming the inverse matrices A⁻¹, D⁻¹, (A#A)⁻¹, (A#D)⁻¹ and (D#D)⁻¹. In a general setting, most of these inverses are impossible to obtain, contrary to that of A, which can be written directly from a genealogy.

(ii) Reformulation of Henderson's approach

Equivalently, as in de los Campos et al. (2008), one may rewrite (1) as

(7)

where a*=A⁻¹a~(0, A⁻¹σ_a²),…, i_dd*=(D#D)⁻¹i_dd~(0, (D#D)⁻¹σ_dd²). Then, the BLUP of any of the transformed genetic effects can be found by solving a system of mixed linear model equations that does not involve inverses of any of the genetic variance–covariance matrices. For example, the β-equation is

and the i_aa-equation is

Once the BLUPs of the (∗) genetic effects are obtained, linear invariance leads, for example, to , and to . A computational difficulty here is that A is typically not sparse, so all A#A, A#D, etc., are not sparse either. Note that the equations for the genetic effects in the ‘reparameterized’ model are similar to those in Henderson (1984). For example, if in an animal model the a-equation

is premultiplied by A, one obtains

(iii) Bayesian implementation of the reformulated model

In the standard Bayesian linear model (Wang et al., 1993, 1994; Sorensen & Gianola, 2002), , etc., are mean vectors of corresponding conditional posterior distributions, for which sampling procedures are well known. Further, with scaled inverse chi-square priors assigned to variance components, Gibbs sampling is straightforward, and does not require forming inverses either. For example, a draw from the conditional posterior distribution of the dominance×dominance variance component is obtained under the alternative parameterization as

where i_dd* is a draw from its conditional posterior distribution, given everything else (Sorensen & Gianola, 2002) and ν_dd and S _dd² are hyper-parameters.

While computations may still be formidable, inversion of the genetic covariance matrices is circumvented. Apart from computational issues, an important question, however, is whether or not the CK–Henderson construct can cope with complex genetic systems effectively.

3. Confronting complexity

Dealing with non-additive genetic variability may be much more difficult than what equations (1) and (3) suggest. Theoretically, at least in CK, epistatic variance can be partitioned into orthogonal additive×additive, additive×dominance, dominance×dominance, etc., variance components, only under idealized conditions. These include linkage equilibrium, absence of mutation and of selection, and no inbreeding or assortative mating (Cockerham, 1954; Kempthorne, 1954). All these assumptions are violated in mature and in breeding programmes. Also, estimation of non-additive components of variance is very difficult, even under standard assumptions (Chang, 1988), leading to imprecise inference. Linkage disequilibrium induces covariances between different types of effects; algebraically heroic attempts to deal with this problem are in Weir & Cockerham (1977) and Wang & Zeng (2006). Gallais (1974) derived expressions aimed to describe the impact of linkage disequilibrium on partition of genetic variance. The number of parameters is large, which makes estimation of the needed dispersion components intractable.

The question of whether or not standard random effects models for quantitative traits can account accurately for non-linear (non-additive) relationships between infinitesimal genotypes and phenotypes remains open. It is argued subsequently that parametric models cannot handle well complexity resulting from interactions between the hundreds or even thousands of genes expected to affect multifactorial traits, such as liability or resistance to disease. This was discussed by Templeton (2000), Gianola et al. (2006) and Gianola and van Kaam (2008).

In the standard theory for non-additive gene action, e.g. epistasis, interaction effects enter linearly into the phenotype, that is, the partial derivative of the model with respect to any effect, be it additive or of an interactive nature, is a constant that does not involve any genetic effect. This theory is unrealistic in non-linear systems, such as those used in metabolic control theory (Bost et al., 1999). Arguing from the perspective of ‘perturbation in analysis’, Feldman & Lewontin (1975) advanced the argument that linear (i.e. analysis of variance) models are ‘local’, presumably in the following sense. Suppose that the expected value of phenotype y, given some genetic values G (genotypes are discrete but their effects are continous), is some unknown function of effects of genotypes at L loci represented as f(G ₁, G ₂, …, G _L). A second-order approximation to the surface of means yields

where is the mean value of the effect of locus i (typically taken to be 0); G=(G ₁, G ₂, …, G _L)′ is a column vector and Ḡ=E(G); f′() is the row vector of first derivatives of f(.) with respect to G and

is the matrix of second derivatives, both evaluated at Ḡ in the approximation.

A linear-on-variance-components decomposition of variability may not be enlightening at all when a phenotype results, say, from a sum of sine and cosine waves. To illustrate, suppose that non-linearity (non-additivity) enters as E (y|G ₁, G ₂, G ₃)=G ₁ exp (G ₂G ₃). A log-transformation of the expected phenotypic value of individuals with genetic value (G ₁, G ₂, G ₃) would suggest linearity (in a log-scale) with respect to locus 1, and interaction between loci 2 and 3. For this hypothetical model, the first derivatives are

and the matrix of second derivatives is

In the neighbourhood of 0 for all effects, the local approximation yields

(8)

On the other hand, near 1 the local approximation produces

(9)

A local approximation near 0 suggests that phenotypes are linear on the effect of locus 1, while an approximation in the neighbourhood of 1 points towards ‘dominance’ at loci 2 and 3, and at 2 – factor epistasis involving loci (1, 2), (1, 3) and (2, 3). This relates to work by Kojima (1959), who studied conditions for equilibria when a fitness surface was affected by epistasis. He assumed free recombination, random mating and constancy of genotypic values. The arguments leading to (8) and (9) indicate that constancy of genotypic values is not tenable when local approximations are used to study the surface.

To the extent that a linear model provides a ‘local’ approximation only, it may not be surprising why attaining an understanding of epistasis within the classical paradigm has been elusive. Kempthorne (1978) disagrees with this view, however, although his argument seems more a defence of the technique of the analysis of variance (ANOVA) per se than of the lack of ability of a linear model to describe complex, interacting, systems. However, even if a linear model holds at least locally, a standard fixed effects ANOVA of a highly dimensional, multifactorial system is not feasible, because one ‘runs out’ of degrees of freedom (df) (Gianola et al., 2006). Nevertheless, methodologies for dealing with complex epistatic systems are becoming available and these include, for example, machine learning, regularized neural networks (Lee, 2004), neural networks optimized with grammatical evolution computations (Motsinger-Reif et al., 2008) and non-parametric regression (e.g. Gianola et al., 2006; Gianola & van Kaam, 2008). The latter is discussed in what follows.

4. Non-parametric breeding value

There is a wide collection of procedures for non-parametric regression available. In particular, statistical models based on RKHS have been useful, inter alia, for regression (Wahba, 1990) and classification (Vapnik, 1998). In many respects, RKHS regression is the ‘mother’ of non-parametric functional data analysis, since it includes splines, hard and soft classification and even best linear unbiased prediction (de los Campos et al., 2008) as special cases. Contrary to many ad-hoc forms of non-parametric modelling, RKHS is a variational method based on maximizing penalized likelihoods over a rich space of functions defined on a Hilbert space. Also, it uses flexible kernels, which can be adapted to many different circumstances, and allows for varying classes of information inputs, e.g. pedigrees, continuous valued covariates and molecular markers of any type. The Bayesian view of RKHS regression has been used to motivate the methodology using Gaussian processes (Rasmussen & Williams, 2006). In quantitative genetics, Gianola et al. (2006) and Gianola & van Kaam (2008) suggested this approach for incorporating information on dense whole-genome markers into models for prediction of genetic value of animals or plants for quantitative traits. The dense molecular data (e.g. SNPs) enter as covariates into a kernel (incidence) square matrix whose dimension is equal to the number of individuals with genotype information available. Thus, the dimensionality of the problem is reduced drastically, from that given by the number of SNPs to the number of individuals genotyped, typically much lower. González-Recio et al. (2008a, b) present applications to mortality rate and feed conversion efficiency in broilers. In the absence of kernels based on substantive theory, genetic interactions are dealt with in RKHS in some form of ‘black box’. However, the focus of non-parametric regression is prediction rather than inference, an issue that will be emphasized in this paper. Relationships between RKHS regression and classical models of quantitative genetics are discussed by de los Campos et al. (2008). A question that arises naturally in plant and animal breeding is whether or not measures of breeding value can be derived from an RKHS regression model.

Briefly, a RKHS regression can be represented in terms of the linear model

(10)

where Xβ is as before; K_h={k(x_i, x_j, h)} is an n×n symmetric positive-definite matrix of kernel entries dependent on covariates (e.g. SNPs) x_i and x_j (subscripts i and j define the individual whose phenotype is considered and some other individual in the sample, respectively), and possibly on a set of bandwith parameters h which must be tuned in some manner, as noted below. González-Recio et al. (2008a) discuss kernels that do not involve any h. Further, α in (10) is a set of non-parametric regression coefficients, one for each individual in the sample of data. For instance, if SNP genotypes enter into k(x_i, x_j, h), then K_hα can be construed as a vector of genetic effects marked by SNPs.

The choice of kernel k (x_i, x_j, h) is absolutely critical for attaining good predictions in RKHS regression, and it can be addressed via some suitable model comparison, e.g. by cross-validation. Each kernel is associated with a space of functions, and de los Campos et al. (2008) describe conditions under which a kernel may be expected to work well. For example, if a Hilbert space of functions associated with a given kernel spans functions of additive and non-additive genetic effects, the model would be expected to capture such effects. Otherwise, predictions may be very poor, even worse than those attained with a linear model (which, in some cases will also be a RKHS regression).

It can be shown that the solution to the RKHS regression problem is obtained by assuming that α~N(0, K_h⁻¹σ_α²), where σ_α² is a variance component entering into the smoothing parameter , and by solving

(11)

The method can be implemented in either a BLUP – residual maximum likelihood (REML) context, given h assessed previously by cross-validation or generalized cross-validation (Craven & Wahba, 1979), or in a fully Bayesian manner, with all parameters assigned prior distributions.

The predicted genetic or genomic value is then K_hᾶ, as illustrated by González-Recio et al. (2008a) in an analysis of broiler mortality in which sires had been genotyped for thousands of SNPs, and by González-Recio et al. (2008b) in a similar study of food conversion ratio conducted in the same population. Their results suggest that RKHS regression using SNP information can produce more reliable prediction of current and future (offspring) phenotypes than the standard parametric additive model used in animal breeding. Also, albeit not unambiguously, some RKHS specifications with filtered SNPs tended to outperform parametric Bayesian regression models in which additive effects of all SNPs had been fitted.

In what follows, the RKHS machinery is used to develop non-parametric measures of breeding value.

(i) Infinitesimal non-parametric breeding value

Assume that some kernel matrix K_h has been found satisfactory, in the sense mentioned above. Given the context, vector K_hα in (10) is a molecularly marked counterpart of g=a+d+i_aa+i_ad+i_dd in (1). Consider now the positive-definite kernel matrix K=A+D+(A#A)+(A#D)+(D#D). This is positive-definite, because it consists of the sum of positive-definite matrices, and leads to the RKHS model

(12)

with α~N(0, K⁻¹ σ_α²). Here, K does not depend on any bandwidth parameters, which simplifies matters greatly, relative to the parametric specification (7), in which 6 variance components enter into the problem. To obtain the RKHS predictor of genetic value (and with only 2 variance components intervening), note that solving (11) is equivalent to solving

and the genetic value is predicted as

(13)

Here, the ‘non-parametric’ breeding value would be Aα, and its predictor is Aᾶ. The posterior distribution of Aα can be arrived at by drawing samples from the posterior distribution of α.

(ii) SNP-based non-parametric breeding value

Let now K_h be a kernel matrix whose entries depend on a string x of SNP genotypes or haplotypes available for a set of individuals having phenotypic records. For example, Gianola & van Kaam (2008) consider the Gaussian kernel

(14)

Alternatively, one could use as kernel

where x_jk is the SNP string for chromosome k in individual j, and h ₁, h ₂, … , h _k are chromosome-specific positive bandwidth parameters.

Irrespective of the kernel adopted, let now

where a_K is defined as non-parametric breeding value, and γ is independent of a_K. With A being the numerator relationship matrix between the individuals having SNP information, use of the reasoning leading to (4) produces

(15)

where σ_a² is some estimate of additive genetic variance for the trait in question. Thus, the non-parametric breeding value of individual i in the sample has the form

so that, if individuals are unrelated, , with variance , where kⁱⁱ is the ith diagonal element of K⁻¹.

Note that

(16)

This is equal to Aσ_a ² only if K=A and σ_α²=σ_a²; in this case, no molecular information is used at all. Suppose now that A=I, so that individuals are genetically unrelated. However, non-parametric breeding values turn out to be correlated, since

(17)

This can be interpreted in the following manner, using a Gaussian kernel to illustrate. Here, the entries of the kernel matrix have the form

and is maximum when x_ik=x_jk, that is, when individuals i and j have the same genotype at marker k. This means that the elements of the kernel matrix (taking values between 0 and 1) will be larger for pairs of individuals that are more ‘molecularly alike’, even if unrelated by line of descent. This molecular similarity is then propagated into (16) and (17).

Given a point estimate of α(ᾶ), e.g. the posterior mean, the non-parametric estimate of breeding value is

(18)

If computations are carried out in a fully Bayesian Monte Carlo context (with m samples drawn from the posterior distribution), the posterior mean estimate is

and the uncertainty measure (akin to prediction error variance–covariance in BLUP) is

Non-parametric breeding values of individuals that are not genotyped can be inferred as follows. Let A_(−,+) be the additive relationship matrix between individuals that are not genotyped (−) and those which are genotyped (+). The non-parametric breeding value of non-genotyped individuals would be

with its point estimate being

Likewise, non-parametric breeding values of yet-to-be genotyped and phenotyped progeny would be

where A(_prog,+) is the additive relationship matrix between progenies and individuals with data.

(iii) Semi-parametric breeding value

Another possibility consists of treating additive genetic effects in a parametric manner and non-additive effects in the RKHS framework, as suggested by Gianola et al. (2006). Let now the model be

(19)

where u~(0, Aσ_a²) and Kα as before, but with

where K_h is a positive-definite matrix of Gaussian kernels with entries as in (14), so that SNP information is used. Since K is the sum of Hadamard products of positive-definite matrices, it is positive-definite as well and, hence, it is a valid kernel for RKHS regression. The estimating equations (can be rendered symmetric by premultiplying the α-set of equations by K) take now the form

where û is the predicted breeding value, and

is the predicted (SNP-based) non-additive genetic value. For instance, (D#K_h)ᾶ is interpretable as a predicted dominance value, and so on. A natural extension of this consists of removing u from (19) and then taking as kernel matrix

so that the predicted additive breeding value would now be (A#K_h)ᾶ, much along the lines of (18). This would give a completely non-parametric treatment of prediction of additive and non-additive genetic effects using dense molecular information.

5. Illustration: infinitesimal model

The toy example in Henderson (1985) is considered. The problem is to infer additive and dominance genetic effects of five individuals using data consisting of phenotypic records for only four (subject 1 lacks a record). The data and model are:

Although 1 does not have a record, all additive and dominance effects are included in the model, zeroing out the incidence of a ₁ and d ₁ via a column of 0s in the Z incidence matrix. Henderson (1985) assumed that the additive, dominance and residual variance components were σ_a²=5, σ_d²=4 and σ_e²=20, so that the phenotypic variance was 29, and used the additive and dominance relationship matrices

Application of BLUP leads to the unique solutions (rounded at the third decimal)

The predicted total genetic value is ĝ=â+, yielding

Alternatively, a RKHS representation as in (12) is used now, with K=A+D as kernel matrix, which does not involve any bandwidth parameters. One has

The RKHS model (using the part of K pertaining to individual with records) is

The value adopted for the sole smoothing parameter is σ_α²=σ_a²+σ_d²=9, and the estimating equation (11) and the solution to the RKHS mixed model equations is

The non-parametric additive genetic effects are inferred using formula (18)

and the non-parametric estimates of dominance effects are

The total genetic value is estimated at

Suppose now that one wishes to predict future performance of these five individuals, and assume that the record of 1 will be made under the same conditions as those for the record of 2, with all conditions remaining the same. The model for future records f is then

for the parametric model, whereas that for the RKHS treatment is

Above M_. and θ_. denote incidence matrices and regression coefficients, respectively; P and K indicate parametric and RKHS treatments, respectively. Using a standard Bayesian argument (Sorensen & Gianola, 2002), the mean vector and covariance matrix of the predictive distributions are

(20)

where C is the corresponding coefficient matrix for each of the procedures; the 5×5 matrix I_fσ_e² conveys uncertainty stemming from the fact that the future records (f) have not been realized yet. For the two procedures, means and standard deviations of the predictive distributions are

6. Illustration: non-linear system

(i) Two-locus model

A hypothetical system with two biallelic loci was simulated. It was assumed that phenotypes were generated according to the rule

(21)

where α_i (β_i) and α_j (β_j) are effects of alleles i and j at the α (β) locus. The system is nonlinear on allelic effects, as indicated by the first derivatives of the conditional expectation function with respect to the αs or βs. For instance,

The α-effects were two random draws from an exponential distribution with mean value equal to 2; the first draw was assigned to allele A and the second to allele a. The βs were drawn from a Weibull (2,1) distribution, having median, mean and mode equal to 0·347, 0·887 and 0·25; the first (second) deviate was the effect of allele B (b). The non-linearity of the system is illustrated in Fig. 1, where the derivative of the model for the expected phenotype with respect to the effect of Weibull allele b is plotted, with all other alleles evaluated at the values drawn. Variation in values of b produces drastic modifications near the origin, but phenotypes are essentially insensitive for b>1·5, even though these values are plausible in the Weibull process hypothesized.

Fig. 1. Rate of change of the expected phenotypic value with respect to the Weibull variable β_j.

Residuals were drawn from the normal distribution N (0, 20), and added to (21) to form phenotypes. The resulting phenotypic distribution is unknown, because y is a non-linear function of exponential and Weibull variates, plus an additive normally distributed residual. There were 5 individuals with records for each of the AABB, AABb, AAbb genotypes; 20 for each of AaBB, AaBb and Aabb and 5 for each of aaBB, aaBb and aabb. Thus, there were 90 individuals with phenotypic records, in total.

Since there are nine distinct mean genotypic values, variation among their average values can be explained completely with a linear model on 9 df; in the standard treatment, these df correspond to an overall mean, additive effects (2 df), dominance (2 df) and epistasis (4 df). Due to non-linearity, it is not straightforward to assess the proportion of variance ‘due to genetic effects’ using a random effects model based on the Weibull and exponential distributions, although an approximation can be arrived at. A linear approximation of the phenotype of an individual yields

where ᾱ and are the means of the exponential and Weibull processes, respectively, and e~N(0, σ²). An approximation to heritability is

(22)

For the situation simulated, ྱ = 2, , Var(α)=4, Var(β)=0·785, so that for σ²=1000, 100, 20, 10 one obtains h ²≈0·13, 0·60, 0·88 and 0·94, respectively. The simulation produced a high penetrance trait, that is, one with heritability near 0·88, which provides a challenge to any non-parametric treatment.

(ii) RKHS modelling

Among the many possible candidate kernels available, an arbitrarily chosen Gaussian kernel was adopted for the RKHS regression implementations, using as covariate a 2×1 vector consisting of the number of alleles at each of the two loci, e.g. x _AA=2, x _Aa=1 and x _aa=0. For example, the kernel entry for genotypes AABB and AAbb is

and the 9×9 kernel matrix for all possible genotypes (labels for genotypes are included, to facilitate understanding of entries) is

There are only five distinct entries, a result of the measure of ‘allelic disimilarity’ adopted, e.g. e ^−8/h stems from disimilarity in 2 alleles at each of the A and B loci (AABB versus aabb). Figures 2 and 3 depict values of the kernel as a function of the bandwidth parameter h at different levels of S, which enters into exp(−S/h). Values of h larger than 10 (Fig. 2) produce strong ‘prior correlations’ between genotypes; also, the kernel matrix becomes more poorly conditioned as h increases. After evaluating the kernel matrix at h=6, 4, 2 and 1·75, it was decided to adopt h=1·75 as the bandwidth parameter, producing 6 unique entries in the K matrix: 1·0 (diagonal elements, the two individuals have identical genotypes); 0·565 (3 alleles in common in a pair of individuals); 0·319 (2 alleles in common, 1 per locus) or 0·102 (2 alleles in common at only one locus); 0·06 (1 allele in common) and 0·01 (no alleles shared).

Fig. 2. Kernel value against the bandwidth parameter h. Curves, from top to bottom, correspond to S=1, 2, 4, 5, 8.

Fig. 3. Kernel value against the bandwidth parameter h. Curves, from top to bottom, correspond to S=1, 2, 4, 5, 8.

(iii) Fitting the RKHS model to the means

A ‘means’ RKHS regression was fitted, including an intercept (β). The model for the 9×1 vector of averages was

where the αs are the non-parametric regression coefficients; K_1·75 is the 9×9 kernel matrix with h=1·75 and is the mean of the residuals pertaining to observations of individuals with genotype i. As before, the assumptions were , and, for the 9×1 vector of average residuals, . Here, N=Diag {5, 5, 5, 20, 20, 20, 5, 5, 5} is a 9×9 diagonal matrix.

The solution to the RKHS regression problem was obtained as

The model was fitted for each of the following values of the shrinkage ratio : 100 (strong shrinkage towards 0), 15, 1, and ; with λ=0 there is no shrinkage of solutions at all. For each of these values, the residual sum of squares

and the weighted residual sum of squares

were computed, where 1 is a 9×1 vector of ones. Also, the effective number of parameters, or model df (Ruppert et al., 2003), was assessed as

where

and is a diagonal matrix containing the square roots of the entries of N. Note that, when , the estimating equations have an infinite number of solutions, as only 9 parameters are estimable. In this case, a solution is obtained by using a generalized inverse of the coefficient matrix, yielding

Here, the fitted values are so the model ‘copies’ the data, and the fit is perfect. Note that is the least-squares estimate of μ_ij=γ_i+δ_j+(γδ)_ij where γ_i (i=AA, Aa, aa) and δ_j (j=BB, Bb, bb) are main effects of genotypes at loci A and B, respectively; this model accounts for 8 df ‘due to’ additive and dominance effects at each of the two loci, and additive×additive, additive×dominance, dominance×additive and dominance×dominance interactions.

For the sake of comparison, the following fixed effects model was fitted as well

where, as before, γ_i (i=AA, Aa, aa) are main effects of genotypes at the A-locus, and δ_j (j=BB, Bb, bb) are the counterparts at the locus B. The linear model was parameterized as

(23)

This model has 5 free parameters, interpretable as an ‘overall mean’ plus additive and dominance effects at each of the 2 loci. Estimates of estimable functions of fixed effects were obtained from a weighted least-squares approach with solution vector where (X′NX)⁻ is a generalized inverse of X′NX and X is the 9×7 incidence matrix given above. Weighted residual sums of squares were computed as . Mean values of A-locus and B-locus genotypes were estimated as

and dominance effects were inferred as

This analysis would create the illusion of non-additivity and overdominance at the A and B loci, respectively, but without bringing light with respect to the non-linearities of (21). While this model does not have mechanistic relevance, it has predictive value, an issue which is illustrated below.

Table 1 gives estimates of the intercept β and of the RKHS regression coefficients α for each of the shrinkage ratio values employed, using h=1·75 as bandwidth parameter in all cases. The residual sum of squares (weighted and unweighted) and the effective number of parameters, or model df, are presented as well. As λ decreased from 10² to 10⁻², model fit improved, but the efective number of parameters increased from 1·46 to 8·95, near the maximum of 9. The implementations with λ=15⁻¹ and λ=10⁻² essentially produced a ‘saturated’ model, that is, one that fits to the means perfectly (as it is the case when λ=0). Large values of the variance ratio (15, 100) produced excessive shrinkage, as indicated by the small values of the RKHS regression coefficients α, and small values of the variance ratio tended to overfit, as noted. On the other hand, the linear 2-locus model on additive effects, with 5 parameters, had residual sum of squares, SSR=19·98, and weighted residual sum of squares, SSR=99·88. These values are more than 5 times those attained with the RKHS regression implementation with a variance ratio of 1 (SSR=3·55, WSSR=19·12), which are shown in Table 1.

Table 1. Estimates of the intercept (β) and of non-parametric regression coefficients (α_i) for each of the values of the variance ratio (λ=σ_e²/σ_α²) employed. RSS and WRSS are the residual and weighted residual sums of squares, respectively; df gives the model df, or effective number of parameters fitted

A more important issue, at least from the perspective taken in this paper, is ‘out of sample’ predictive ability. To examine this, 3 new (independent) samples of phenotypes were generated, assuming the residual distribution N (0, 20), as before, and with 5 individuals per genotype, i.e. there were 45 subjects in each sample. The predictive residual sums of squares (unweighted and weighted, using 5 as weight) were calculated, using the fitted values from the training sample employed to compute statistics in Table 1, and the ‘new sample’ phenotypes. The 2-locus additive model was also involved in the comparison. Results are shown in Table 2, where entries are the average, minimum and maximum values of the predictive sum of squares over the 3 new samples. As expected, predictive residual sum of squares were much larger than those observed in the training sample, notably for implementations in which overfitting to the training data was obvious see Table 1). The specifications producing strong shrinkage towards 0 in the training sample had the worse predictive performance, whereas that with had the best performance, on average, albeit close to the 2-locus additive model. It is not surprising that small variance ratios led to reasonable predictive performance, be-cause the simulation mimicked high penetrance, i.e. phenotypes are very informative about genotypes, this being so because the approximate heritability (22) for a residual variance of 20 is about 0·88, as already noted.

Table 2. Predictive (weighted and unweighted by the number of individuals per geno-type) residual sums of squares for each of the variance ratios (λ=σ_e²/σ_α²) employed in the non-parametric regression implementation, and for the two-locus model with main effects of genotypes at each of the loci. Entries are average (boldface) from three predictive samples, with minimum and maximum values over samples in parentheses

For this example, the 2-locus additive model is untenable, at least mechanistically, yet it had a reasonable predictive performance. Apart from genetic considerations discussed later, this is because any of the models considered here can be associated with a linear smoother, with predictions having the form

for some matrix L and where l_i′ is its ith row. For the 2-locus model,

whereas for any of the RKHS specifications

Hence, all predictors can be viewed as consisting of different forms of averaging observations, with the RKHS-based averages having optimality properties, in some well defined sense (Kimeldorf & Wahba, 1971; Wahba, 1990; Gianola & van Kaam, 2008; de los Campos et al., 2008). Since the data consisted of phenotypic averages, i.e. a non-parametric averaging method where a ‘bin’ is a given genotype, this represented a challenge for any of the smoothers, including the 2-locus additive model. However, some of the smoothers (e.g. RKHS with λ=1 and the 2-locus additive model) met the challenge satisfactorily.

(iv) Fitting the RKHS model to individual observations

The RKHS regression model (using the Gaussian kernel with h=1·75) was fitted again to the 90 data points in a newly simulated sample used to estimate (train) the intercept and the nine non-parametric coefficients. In addition, the following parametric specifications were fitted: additive model (3 location parameters: intercept and additive effects of the A and B loci), additive+dominance model (2 extra parameters corresponding to dominance effects at the two loci), and additive+dominance+epistasis (4 additional df pertaining to additive×additive, additive×dominance, dominance×additive and dominance×dominance interactions). The corresponding regression coefficients for the parametric models were ‘shrunken’ using a common variance ratio, corresponding to each of the 15 λ values employed in the RKHS fitting. Subsequently, 100 predictive samples of size 45 each were simulated, and the realized values were compared against the predictions obtained from the sample used to train either the RKHS regression or the three parametric models.

Figure 4 displays the average squared residual (over the 90 data points in the training sample) for each of the four models fitted. As expected, the RHKS and the additive+dominance+epistatic models fitted the data best at λ=0 (no shrinkage), because of having a larger effective number of parameters (9) than the additive (3) and additive+dominance (5) models. When effects were gradually shrunken (λ increased from 0 to 20), the parametric models maintained their relative standings, whereas RKHS voyaged through all three models, eventually giving a very poor fit, due to oversmoothing. The trajectory of the effective df is shown in Fig. 5, where the ability of RKHS to explore models of different degrees of complexity becomes clear. The out-of-sample predictive performance of the four models is shown in Fig. 6. The simplest model, that is, one with additive effects at each of the two loci fitted, had the best predictive performance, even better than the two additional parametric specifications although the simulated gene action was not additive at all! RKHS regression was competitive, but its predictive ability deteriorated markedly when λ was greater than 5 in the training sample.

Fig. 4. Average (over 90 data points) squared residual for four models plotted to the training sample (RKHS=RKHS regression with Gaussian kernel and h=1·75) for each value of the smoothing parameter λ.

Fig. 5. Effective df for four models plotted to the training sample (RKHS=RKHS regression with Gaussian kernel and h=1·75) at each value of the smoothing parameter λ.

Fig. 6. Average (over 100 samples with 45 realized observations in each) squared prediction error for four models plotted to the predictive sample (RKHS=RKHS regression with Gaussian kernel and h=1·75) for each value of the smoothing parameter λ.

How does one explain the paradox that a simple additive model had better predictive performance when gene action was non-linear, as simulated here? In order to address this question, consider the ‘true’ mean value of the 9 genotypes simulated:

• 

The ‘corrected’ sum of squares among these means is 125·23. A fixed effects ANOVA of these ‘true’ values (assuming genotypes were equally frequent) gives the following partition of sequential sum of squares, apart from rounding errors: (i) additive effect of locus A: 82·8%; (ii) additive effect of locus B after accounting for A: 7·06%; (iii) dominance effects of loci A and B: 4·2%, and (iv) epistasis: 6·2%. Thus, even though the genetic system was non-linear, most of the variation among genotypic means can be accounted for with a linear model on additive effects. The additive model had the worst fit to the data (even worse than the models that assume dominance and epistasis) and, yet, it had the best predictive ability, followed by RKHS for (roughly) 0·5<λ<3.

The simulation was also carried out at larger values of the residual variance, σ_e²=100 and 500. Again, the purely additive model had the best predictive performance, but the difference between models essentially disappeared for λ>50. In particular, the average squared prediction error of RKHS was larger than that for the additive model by 5, 3 and 1% for σ_e²=20, 100 and 500, respectively, when evaluated at the ‘optimal’ λ in each case.

It should be noted that the RKHS implementation used here was completely arbitrary. For example, the kernel chosen was not the result of any formal model comparison, so predictive performance could be enhanced by a more judicious choice of kernel. As noted earlier, the choice of a good kernel is critical in this form of non-parametric modelling.

7. Conclusion

Inference about genotypes and future phenotypes for a complex quantitative trait was discussed in this paper. In particular, it was argued that the Kempthorne–Cockerham theory for partioning variance into additive, dominance and epistatic components has doubtful usefulness, because practically all assumptions required are violated in artificial and natural populations. This theory is probably illusory when genetic systems are complex and non-linear, in agreement with views in Feldman & Lewontin (1975) and Karlin et al. (1983). Further, an ANOVA-type decomposition is inadequate for a non-linear system (because the ANOVA model is linear in the parameters), and unfeasible in a highly dimensional and interactive genetic system involving hundreds or thousands of genes. As a minimum, the ANOVA treatment would encounter a huge number of empty cells, extreme lack of orthogonality, and high-order interactions would be extremely difficult to interpret, in the usual sense. Last but not the least, the ANOVA model would require more df than the number of data points available for analysis.

For these reasons, a predictive approach was advanced in this paper, focused on the use of non-parametric methods, especifically RKHS regression. Use of this methodogy in conjunction with standard theory of quantitative genetics leads to non-parametric estimates of additive, dominance and epistatic effects. These ideas were illustrated using a stylyzed example in Henderson (1985), and it was shown how additive, dominance and total genetic values can be predicted using a single smoothing parameter (in addition to the residual variance) coupled with kernels based on substantive theory, a point that is also made in de los Campos et al. (2008). A non-linear 2-locus system was simulated as well, to illustrate the RKHS approach, which was found to have a better out-of-sample predictive performance of means than the standard 2-locus fixed effects model (with and without epistasis).

On the other hand, a 2-locus model with estimates of additive effects shrunken to different degrees had the best performance when predicting future individual observations. This was explained by the observation that most of the variation among genotypic means could be accounted for by ‘additive effects’. This is consistent with theoretical and empirical results presented by Hill et al. (2008), who gave evidence that, even in the presence of non-additive genetic action, most variance is of an additive type. Even though molecular geneticists view the additive model as irritatingly reductive, our results give reassurance to a common practice in animal breeding, i.e. predict genetic values using additive theory only. It is unknown, however, to what extent these results (from a predictive perspective) carry to more complex systems, difficult to be described suitably with naive linear models. In such situations, RKHS regression could be valuable, because of its ability to navigate through models of different degrees of complexity and, as shown here, it can be very competitive when the smoothing parameters are tuned properly.

In conclusion, it is felt that the non-parametric methods discussed here coupled with machine learning procedures, such as in Long et al. (2007) and Long et al. (2008), could play an important role in quantitative genetics in the post-genomic era.