The regression function:

Consider a stylized situation in which each of a series of individuals possesses a measurement for some quantitative trait denoted as y, as well as information on a possibly massive number of genomic variables, such as SNP “genotypes,” represented by a vector x. In the main, x is treated as a continuously valued vector of covariates, even though SNP genotypes are discrete (coding is done via dummy variates). Also, x could represent gene expression measurements from microarray experiments; here, it would be legitimate to regard this vector as continuous. Although gene expression measurements are typically regarded as response variables, there are contexts in which this type of information could be used in an explanatory role (Mallick et al. 2005).

Let the relationship between y and x be represented as(1)where y_i is a measurement, such as plant height or body weight, taken on individual i, x_i is a p × 1 vector of dummy SNP or microsatellite covariates observed on i, and g(.) is some unknown function relating these genotypes to phenotypes. Define g(x_i) = E(y_i|x_i) as the conditional expectation function, that is, the mean phenotypic value of an infinite number of individuals, all possessing the p-dimensional genotype x_i. e_i ∼ (0, σ²) is a random residual, distributed independently of x_i and with variance σ².

The conditional expectation function is(2)Following Silverman (1986), consider a nonparametric kernel estimator of the p-dimensional density of the covariates,(3)where is a kernel function and h is a window width or smoothing parameter. In (3), x is the value (“focal point”) at which the density is evaluated and x_i (i = 1, 2,…,n) is the observed p-dimensional SNP genotype of individual i in the sample. Hence, (3) estimates population densities (or frequencies). If is to behave as a multivariate probability density function, then it must be true that the kernel function is positive and that the conditionis satisfied. This implies thatSimilarly, and assuming that a single h parameter suffices, one can estimate the joint density of phenotype and genotypes at point (y, x) aswhere is also a kernel function; again, y_i is the observed sample value of variable y in individual i. The numerator of (2) can be estimated as(4)In (4), let , so that dy = hdz andThe kernel function is typically a proper probability density function chosen such that and . If so, the preceding expression is equal to y_i, so that (4) becomes(5)

Returning to (2), one can form the nonparametric estimatorwhich, upon replacing the numerator and denominator by (5) and (3), respectively, takes the form(6)whereis a weight that depends on the kernel function and window width chosen and on the x_i (i.e., genotypes) observed in the sample. The linear combination of the observations (6) is called the Nadaraya–Watson estimator of the regression function (Nadaraya 1964; Watson 1964). As seen in (6), the fitted value at coordinate x is a weighted average of all data points, with the value of the weight depending on the “proximity” of x_i to x and on the value of the smoothing parameter h. For instance, if the kernel function has the Gaussian formthis has a maximum value of when x = x_i and tails off to 0 as the distance between x and x_i increases. The values of w_i(x) decrease more abruptly as . This is the reason why this type of estimator is called “local,” in the sense that observations with x_i coordinates closer to the focal point x are weighted more strongly in the computation of the fitted value .

A specification that is more restrictive than (1) is the additive regression model(7)(Hastie and Tibshirani 1990; Fox 2005), where x_ij is the genotype for SNP j in individual i. In this model each of the “partial regression” functions is two-dimensional, thus allowing exploration of effects of individual SNPs on phenotypes, albeit at the expense of ignoring, e.g., epistatic interactions between genotypes. A preliminary examination of relationships can be done via calculation of Nadaraya–Watson estimators of each of the g_j(.) aswhere the kernels are unidimensional. Naturally, one may wish to account for interactions between SNPs, so this type of analysis would be merely exploratory. A specification that is intermediate between (7) and (1) could include sums of single, pairwise, tripletwise, etc., SNP regression functions.

Impact of window width:

The scatter plot shown in Figure 1, from Chu and Marron (1991), consists of data on log-income and age of 205 people. The solid line in Figure 1A is a moving weighted average of the points, with weights proportional to the curve at the bottom. Chu and Marron (1991) regard the dip in average income in the middle ages as “unexpected.” Figure 1B gives results from three different smooths at window widths of 1, 3, and 9. When h = 1, the curve displays considerable sampling variation or roughness. On the other hand, when h = 9, local features disappear, because points that are far apart receive considerable weight in the fitting procedure. If the dip is not an artifact, the oversmoothing results in a “bias.” Hence, h must be gauged carefully.

• Download figure
• Open in new tab
• Download powerpoint

Figure 1.—

Impact of window width on the regression function of log-income on age (Chu and Marron 1991). (A) Scatter plot and (B) smooths for earning power data. Kernel is N(0, 1). Window widths are represented by curves: solid curves, h = 3; dotted curves, h =1; dashed curves, h = 9.

Marron (1988) and Chu and Marron (1991) discuss data-driven procedures for assessing h. Silverman (1986) gives a discussion in the context of density estimation, whereas Mallick et al. (2005) consider Hilbert spaces kernel regression, with h treated as an unknown parameter. A conceptually simple and intuitively appealing procedure is cross-validation (e.g., Schucany 2004). For instance, in the “leave-one-out” procedure, the datum for case i, that is (y_i, x_i), is deleted, and a fit is carried out on the basis of the other n − 1 cases. Then, the prediction of y_i is formed, where the notation −i indicates that all data other than that for case i are used for estimating the regression function. This process is repeated for all n data points. Subsequently, the cross-validation criterionis formed. A cross-validation estimate of h is the minimizer of CV(h); this is found by carrying out the computations over a grid of h-values. Hart and Lee (2005), however, present evidence of large variability of the leave-one-out estimates of h. Hastie et al. (2001) discuss alternatives based on leaving out 10–20% of the sample values. Ruppert et al. (2003) present procedures for simple calculation of cross-validation statistics.

D. Gulisija, D. Gianola and K. A. Weigel (unpublished results) used another nonparametric procedure, LOESS, to study the relationship between performance and inbreeding in Jersey cows. There is some relationship between LOESS and kernel regression. In LOESS, the number of points contributing to a focal fitted value is fixed (contrary to kernel regression, where the actual number depends on the gentleness of the kernel chosen) and governed by a spanning parameter. This parameter (ranging between 0 and 1) is equivalent to h and dictates the fraction of all data points that contribute toward a fitted value. Figure 2 gives a LOESS fit for protein yield (actually, residuals from a parametric model) and 100 bootstrap replicates, illustrating uncertainty about the regression surface. Without getting into details, note that yield decreases gently at low values of inbreeding, followed by a faster linear decline, and then by an apparent increase. Irrespective of the variability (due to that few animals were either noninbred or highly inbred), neither the change of rate at low inbreeding nor the increase in yield at high consanguinity would be predicted by standard quantitative genetics theory. This is another illustration of how “irregularities” can be discovered nonparametrically, which would remain hidden otherwise.

• Download figure
• Open in new tab
• Download powerpoint

Figure 2.—

LOESS curves of protein yield deviations (EBLUP) in Jersey cows against their inbreeding coefficients (F, %). The thick curve is the fitted regression (second-degree local polynomial, spanning parameter = 0.9); the other curves are 100 bootstrap replicates.

Estimation of linear approximation:

If the kernel function is a probability density function, the nonparametric density estimator (3) will be a density function as well, retaining differentiability properties of the kernel (Silverman 1986). Consider E (y | x) = g(x) and suppose that one is interested in inferring a linear approximation to g(x) near some fixed point x* such as the mean value of the covariates; this leads to a nonparametric counterpart of additive genetic value. From a plant and animal breeding point of view, the concept of breeding value is essential in parametric models, so it seems important to develop a nonparametric counterpart as well. The linear function is(8)whereis the gradient of g(x) at x = x*. This suggests the pseudomodel(9)and the approximate variance decomposition(10)The first term in the expression above can be interpreted as the variance contributed by the linear effect of x in the neighborhood of x*. Further, if x is a vector of genotypes for molecular markers, this would be, roughly, the additive variance “due to” markers.

A nonparametric estimator of (8) is given by(11)where(12)and is an estimate of the vector of first derivatives of g(x) with respect to x, evaluated at x*. The estimate could be the gradient of (6),(13)whereandare vectors whose form depends on the kernel function used. Hence, (11) becomes(14)Likewise, an estimator of the variance contributed to (10) by the linear effect of x (nonparametric additive genetic variance from the SNPs) is given by , where is some estimate of the covariance matrix of x. Finally, the relative contribution to variance made by the linear effect of x(τ²) can be assessed as(15)where is an estimate of σ², such as(16)

Gaussian kernel:

A p-dimensional Gaussian kernel with a single band width parameter h has the formso that the estimator of the density at x is the finite mixture of Gaussians,where, for x_i = [x_i1 x_i2 . x_ip]′, at focal x = [x₁ x₂ . x_p]′

The weights entering into estimator (6) are thenThe fitted surface is given by(17)where y is the n × 1 vector of data points, 1 is an n × 1 vector of ones, and v′(x) is an n × 1 vector with typical element

Consider next the linear approximation in (14). Using (17), writeNow, for a Gaussian kernel,Further,andHenceand

Kernels for discrete covariates:

For a biallelic SNP, there are three possible genotypes at each “locus,” as in stylized Mendelian situations. In a standard (parametric) analysis of variance representation, incidence situations (or additive and dominance effects at each of the loci) are described via two dummy binary variables per locus, and all corresponding epistatic interactions can be assessed from effects of cross-products of these variables. This leads to a highly parameterized structure and to formidable model selection problems.

Consider now the nonparametric approach. For an x vector with p coordinates, its statistical distribution is given by the probabilities of each of the 3^p combinations of binary outcomes. With SNPs, p can be very large (possibly much larger than n), so it is hopeless to estimate the probability distribution of genotypes accurately from observed relative frequencies, and smoothing is required (Silverman 1986). Kernel estimation extends as follows: for binary covariates the number of disagreements between a focal x and the observed x_i in subject i is given bywhere d(.) takes values between 0 and p. For illustration, suppose that one has genotype AABbccDd as focal point and that individual i in the sample has been genotyped as aaBbccDD. The vectors of binary covariates for the focal and observed cases (save for an intercept) are given in the matrix where N (Y) stands for disagreement (agreement). Then, d(x, x_i) = 4, which is twice the number of disagreements in genotypes because there are only 2 “d.f.” per locus. In practice, one should work with a representation of incidences that is free of redundancies. For binary covariates, Silverman (1986) suggests the “binomial” kernelwith ; alternative forms of the kernel function are discussed by Aitchison and Aitken (1976) and Racine and Li (2004). It follows that the kernel estimate of the probability of observing the focal value x is(18)If h = 1, the estimate is just the proportion of cases for which x_i = x; if , every focal point gets an estimate equal to , irrespective of the observed values x₁, x₂,…, x_n.

The nonparametric estimator of the regression function is

Since a discrete distribution does not possess derivatives, the additive genetic value must be defined in the classical sense (e.g., Falconer 1960), that is, by defining contrasts between expected values of individuals having appropriate genotypes. Here, one can use either the vectorial representation or the concept of the additive regression model in (7). Hereinafter, it is assumed that the distribution of x is continuous, and a continuous kernel function is employed throughout, as an approximation.

General considerations:

Consider now a situation for which there might be an operational or mechanistic basis for specifying at least part of a model. For instance, suppose that y is a measure on some quantitative trait, such as milk production of a cow. Animal breeders have exploited to advantage the infinitesimal model of quantitative genetics (Fisher 1918). Vectorial representations of this model are given by Quaas and Pollak (1980), and applications to natural populations are discussed by Kruuk (2004). In this section, we combine features of the infinitesimal model with a nonparametric treatment of genomic data and present semiparametric implementations.

Statistical specification:

Model (1) is expanded as(19)where is a vector of nuisance location effects and u is a q × 1 vector containing additive genetic effects of q individuals (these effects are assumed to be independent of those of the markers), some of which may lack a phenotypic record; w′_i and z′_i are known nonstochastic incidence vectors. As before, g(x_i) is some unknown function of the SNP data. It is assumed that , where is the “unmarked” additive genetic variance and A is the additive relationship matrix, whose entries are twice the coefficients of coancestry between individuals. Let be the n × 1 vector of residuals, and assume that , where is the residual variance. Note that the model implies that(20)and(21)The preceding means that: (1) the offset y_i − g(x_i) follows a standard mixed-effects model, and (2) if and u were known, one could use (6) to estimate g(x_i) employing y_i − w′_iβ − z′_iu as “observations.” This suggests two strategies for analysis of data, discussed below.

Strategy 1—mixed model analysis:

This follows from representation (20). First, estimate g(x_i), for i = 1, 2, … , n, via , as in (6) or, if a Gaussian kernel is adopted, as in (17). Then, carry out a mixed model analysis using the “corrected” data vector and pseudomodelwhere and are incidence matrices of appropriate order. The pseudomodel ignores uncertainty about , since is treated as if it were the true regression (on SNPs) surface.

Under the standard multivariate normality assumptions of the infinitesimal model, one can estimate the variance components and from y* via restricted maximum likelihood (REML) (Patterson and Thompson 1971) and form empirical best linear unbiased estimators and predictors of β and u, respectively, by solving the Henderson mixed model equations(22)(Henderson 1973). The ratio is evaluated at REML estimates of the variance components. Solving system (22) is a standard problem in animal breeding even for very large q, since A⁻¹ is easy to compute. The two-stage procedure could be iterated several times, i.e., use the solutions to (22), to obtain a new estimate of using as “data,” and then update the pseudodata y*, etc.

The “total” additive genetic effect of an individual possessing a vector of SNP covariates with focal value x can be defined as the sum of the additive genetic effect of the SNPs, in the sense of (8), plus the polygenic effect u_i, that is,An empirical predictor of T_i can be formed by adding (14) to the ith component of in the mixed model equations. It is not obvious how a measure of uncertainty about can be constructed using this procedure.

A Bayesian approach can be used instead, using the corrected data as observations. Under standard assumptions made for the prior and the likelihood (e.g., Wang et al. 1993, 1994; Sorensen and Gianola 2002), one can draw samples j = 1, 2, … , m from the pseudoposterior distribution via Gibbs sampling and then form “semiparametric” draws of the total genetic value asfor j = 1, 2, … , m. These can be construed as draws from a semiparametric pseudoposterior distribution; one can readily calculate the means, median, percentiles, and variance of this distribution and produce posterior summaries, leading to an approximation to uncertainty about total genetic value.

Irrespective of whether classical or Bayesian viewpoints are adopted, this approach ignores the error of estimation of g(x), as noted earlier.

Strategy 2—random g(.)function:

The estimate of can be improved as follows. Consider (21) and suppose, temporarily, that and u are known. Then, form the offset and the alternative Nadaraya–Watson estimator of g(x_i):(23)The problem is that the location vectors and u are not observable, so they must be inferred somehow.

Regard now as unknown quantities possessing some prior distribution, which is modified via Bayesian learning into the pseudoposterior process , after observation of the pseudodata y* and of the n × p matrix of observed SNP covariates X. Then, for some band width h, would also be a random variable possessing some pseudoposterior distribution. Let be a draw from the pseudoposterior process , as in the preceding section. Further, regard as a draw from the pseudoposterior distribution of , given y*. Subsequently, estimate features of the pseudoposterior distribution of by ergodic averaging. For example, an estimate of its posterior expectation would bewhere the samples are obtained via a Markov chain Monte Carlo (MCMC) procedure from the pseudoposterior distribution .

The MCMC procedure would consist of sequential, iterative, sampling from all conditional pseudoposterior distributions, as follows:

• Sample from , where ELSE denotes all other parameters, y*, X, and h. Using standard results, this conditional distribution has the formwhereand

• Sample u from . Using similar standard results, the distribution to sample from iswhereand

• Sample the two variance components from the scaled inverse chi-square distributionsandAbove, ν_u, and ν_e, are known hyperparameters of independent scaled inverse chi-square priors assigned to and , respectively.

• Form draws from the pseudoposterior distribution of as

As usual, early draws in the MCMC procedure are discarded as burn-in and, subsequently, m samples are collected to infer features of interest. Note that the procedure termed as strategy 1 gives as an estimate of the relationship between phenotypes and SNP data. In strategy 2, one can use instead, where and are means of the pseudoposterior distributions and , respectively. Under strategy 2, a point predictor of total additive genetic value could bewhere is the posterior mean of u_i.

General:

What follows is motivated by developments in Mallick et al. (2005) for classification of tumors using microarray data. The underlying theory is outside the scope of this article. Only essentials are given here, and foundations are in Wahba (1990, 1999).

Using the structure of (19), consider the penalized sum of squares(24)where, as before, h is a smoothing parameter (possibly unknown) and is some norm or “stabilizer.” For instance, in smoothing splines, is a function of the second derivatives of g(x) integrated between end points that compose the data. The second term in (24) acts as a penalty: if the unknown function g(x) is rough, in the sense of having slopes that change rapidly, the penalty increases. The main problem here is that of finding the function g(x) that minimizes (24). Since is a functional on g(x), this is a variational or calculus of variations problem over a space of smooth curves. The solution was given by Kimeldorf and Wahba (1971) and Wahba (1999), and the minimizer admits the representationwhere is called a reproducing kernel. A possible choice for the kernel (Mallick et al. 2005) is the single smoothing parameter Gaussian function

Mixed model representation:

We embed these results into (19), leading to the specification(25)with the intercept parameter α₀ included as part of . Note that there are as many regressions α_j as there are data points. However, the roughness penalty in the variational problem leads to a reduction in the effective number of parameters in reproducible kernel Hilbert spaces (RKHS) regression, as it occurs in smoothing splines (Fox 2002).

Define the 1 × n row vectorthe n × 1 column vector , j = 1, 2,…, n; and the n × n matrixThen (25) can be written in matrix form as

Suppose, further, that the α_j coefficients are exchangeable according to the distribution . Hence, for a given smoothing parameter h, we are in the setting of a mixed-effects linear model.

Given h, , , and (at a given h, the three variance components may be estimated by, e.g., REML) one can obtain predictions of the polygenic breeding values u and of the coefficients from the solutions to the system(26)The total additive genetic value of an individual possessing a vector of SNP covariates x_i can be defined as(27)whereSince T_i is a linear function of the unknown u_i and α_j effects, its empirical best linear unbiased predictor, assuming known h, can be obtained by replacing these effects by the corresponding solutions from (26).

Incomplete genotyping:

At least in animal breeding, it is not feasible to have all individuals genotyped for SNPs. On the other hand, the number of animals with phenotypic information available is typically in the order of hundreds of thousands, and genotyping is selective, e.g., young bulls that are candidates for progeny testing in dairy cattle production. Animals lacking molecular data are not a random sample from the population, and ignoring this issue may lead to biased inferences. Unless missingness of molecular data is ignorable, in the sense of, e.g., Henderson (1975), Rubin (1976), Gianola and Fernando (1986), or Im et al. (1989), the procedures given below require modeling of the missing data process, which is difficult and may lack robustness. Here, it is assumed that missingness is ignorable, enabling use of likelihood-based or Bayesian procedures as if selection had not taken place (Sorensen et al. 2001). Two ad hoc procedures are discussed, and an alternative approach, suitable for kernel regression, is presented in the conclusion.

Let the vector of phenotypic data be partitioned as , where y₁ (n₁ × 1) consists of records of individuals lacking SNP data, whereas y₂ (n₂ × 1) includes phenotypic data of genotyped individuals. Often, it will be the case that n₁ > p ≫ n₂. We adopt the model(28)For the sake of flexibility, assume that and are mutually independent but heteroscedastic vectors. In short, the key assumption made here is that the random effect α affects y₂ but not y₁ or, equivalently, that it gets absorbed into e₁. With this representation, the mixed model equations take the form(29)If SNP data are missing completely at random and h, , , and are treated as known, then is an unbiased estimator of , and and are unbiased predictors of u and α, respectively. They are not “best,” in the sense of having minimum variance or minimum prediction error variance, because the smooth function of the SNP markers is missing in the model for individuals that are not genotyped (Henderson 1974).

An alternative consists of writing the bivariate modeland then assigning to the polygenic component, the multivariate normal distributionHere, and are additive genetic variances in individuals without and with molecular information, respectively, and is their additive genetic covariance. Computations would be those appropriate for a two-trait linear model analysis (Henderson 1984; Sorensen and Gianola 2002).

Bayesian analysis:

To illustrate, consider the first of the two options in the preceding section. Suppose a kernel has been chosen but that the value of h is uncertain, so that the model unknowns areLet the prior density have the form(30)where H denotes the set of all known hyperparameters (whose values are fixed a priori) and indicates a multivariate normal distribution with appropriate mean vector and covariance matrix. The four variance components are assigned independent scaled inverse chi-square prior distributions with degrees of freedom ν and scale parameters S², with appropriate subscripts. Assign an improper prior distribution to each of the elements of and, as in Mallick et al. (2005), adopt a uniform prior for h, with lower and upper boundaries h_min and h_max, respectively.

Given the parameters, observations are assumed to be conditionally independent, and the distribution adopted for the sampling model is

Given h, one is again in the setting of the Bayesian analysis of a mixed linear model, and Markov chain Monte Carlo procedures for this situation are well known. Under standard conjugate prior parameterizations, all conditional posterior distributions are known, save for that of h. Hence, one can construct a Gibbs–Metropolis sampling scheme in which conditional distributions are used for drawing and a Metropolis update is employed for h.

Location effects are drawn from a multivariate normal distribution with mean vector and covariance matrix(31)and(32)respectively. Likewise, the additive genetic effects can be sampled from a normal distribution centered at(33)and with covariance matrix(34)The conditional posterior distribution of the coefficients is multivariate normal as well, with mean vector(35)and variance–covariance matrix(36)

All four variance components have scaled inverse chi-square conditional posterior distributions and are conditionally independent. The conditional posterior distributions to sample from are(37)(38)(39)and(40)

The most difficult parameter to sample is h. Its conditional posterior density can be represented as(41)whereDensity (41) is not in a recognizable form. However, a Metropolis algorithm (Metropolis et al. 1953) can be tailored for obtaining samples from the distribution . Suppose that the Markov chain is at state h^[t]. A proposal value h* is drawn from some symmetric candidate-generating distribution and accepted with probabilityIf the proposal is accepted, then set h^[t+1] = h*; otherwise, stay at h^[t]. If a Hastings update is employed, instead, an adaptive proposal distribution could be, e.g., a finite mixture of the prior density and of a scaled inverse chi-square density f with mean value h^[t] and variance ,where h_largest and h_smallest are the largest and smallest value of h, respectively, accepted up to time t. The weight w in is assessed such that a reasonable acceptance rate for the Metropolis–Hastings algorithm is attained.

Illustrative example:

Phenotypic and genotypic values were simulated (single replication) for a sample of N unrelated individuals, for each of two situations. The trait was determined either by five biallelic QTL, having additive gene action, or by five pairs of biallelic QTL, having additive-by-additive gene action. Under nonadditive gene action, the additive genetic variance was null, so that all genetic variance was of the additive-by-additive type. Heritability (ratio between genetic and phenotypic variance) in both settings was set to 0.5.

Genotypes were simulated for a total of 100 biallelic markers, including the “true” QTL; all loci were simulated to be in gametic-phase equilibrium. Since all individuals were unrelated and all loci were in gametic-phase equilibrium, only the QTL genotypes and the trait phenotypes would provide information on the genotypic value. In most real applications, the location of the QTL will not be known, and so many loci that are not QTL will be included in the analysis.

A RKHS mixed model was used to predict genotypic values, given phenotypic values and genotypes at the QTL and at all other loci. The model included a fixed intercept α₀ and a random RKHS regression coefficient α_i for each subject; additive effects were omitted from the model, as individuals were unrelated, precluding separation of additive effects from residuals, in the absence of knowledge of variance components. The genetic value g_i of a subject was predicted aswhere and were obtained by solving (26) for this model, using a Gaussian kernel at varying functions of h. The mean squared error of prediction of genetic value (MSEP) was calculated as(42)and a grid search was used to determine the values of h and that minimized (42). To evaluate the performance of as a predictor, another sample of 1000 individuals (“PRED”) was simulated, including genotypes, genotypic values, and phenotypes. This was deemed preferable to doing prediction in the training sample, to reduce dependence between performance and (h, ), whose values were assessed in the training sample The genotypic values of the subjects in PRED were predicted, given their genotypes, using . Genotypic values were also predicted using a multiple linear regression (MLR) mixed model with a fixed intercept and random regression coefficients on the linear effects of genotypes.

Results for the RKHS mixed model are in Table 1, and those for the MLR mixed model are in Table 2, where “accuracy” is the correlation between true and predicted genetic values. When gene action was strictly additive, the two methods (each fitting k = 100 loci) had the same accuracy, indicating that RKHS performed well even when the parametric assumptions were valid. On the other hand, when inheritance was purely additive by additive, the parametric MLR was clearly outperformed by RKHS, irrespective of the number of loci fitted. An exception occurred at k = 100 and N = 1000; here, the two methods were completely inaccurate. However, when N was increased from 1000 to 5000, the accuracy of RKHS jumped to 0.85 (Table 1), whereas MLR remained inaccurate (Table 2). Note that the accuracy of RKHS decreased (with N held at 1000) when k increased. We attribute this behavior to the use of a Gaussian kernel when, in fact, covariates are discrete.