Estimation and testing for the effect of a genetic pathway on a disease outcome using logistic kernel machine regression via logistic mixed models (original) (raw)

• Research article
• Open access
• Published: 24 June 2008

BMC Bioinformatics volume 9, Article number: 292 (2008)Cite this article

• 12k Accesses
• 133 Citations
• Metrics details

Abstract

Background

Growing interest on biological pathways has called for new statistical methods for modeling and testing a genetic pathway effect on a health outcome. The fact that genes within a pathway tend to interact with each other and relate to the outcome in a complicated way makes nonparametric methods more desirable. The kernel machine method provides a convenient, powerful and unified method for multi-dimensional parametric and nonparametric modeling of the pathway effect.

Results

In this paper we propose a logistic kernel machine regression model for binary outcomes. This model relates the disease risk to covariates parametrically, and to genes within a genetic pathway parametrically or nonparametrically using kernel machines. The nonparametric genetic pathway effect allows for possible interactions among the genes within the same pathway and a complicated relationship of the genetic pathway and the outcome. We show that kernel machine estimation of the model components can be formulated using a logistic mixed model. Estimation hence can proceed within a mixed model framework using standard statistical software. A score test based on a Gaussian process approximation is developed to test for the genetic pathway effect. The methods are illustrated using a prostate cancer data set and evaluated using simulations. An extension to continuous and discrete outcomes using generalized kernel machine models and its connection with generalized linear mixed models is discussed.

Conclusion

Logistic kernel machine regression and its extension generalized kernel machine regression provide a novel and flexible statistical tool for modeling pathway effects on discrete and continuous outcomes. Their close connection to mixed models and attractive performance make them have promising wide applications in bioinformatics and other biomedical areas.

Background

The rapid progress in gene expression array technology in the past decade has greatly facilitated our understanding of the genetic aspect of various diseases. Knowledge-based approaches, such as gene set or pathway analysis, have become increasingly popular. In such gene sets/pathways, groups of genes act in concert to accomplish tasks related to a cellular process and the resulting genetic pathway effects may manifest themselves through phenotypic changes, such as occurrence of disease. Thus it is potentially more meaningful to study the overall effect of a group of genes rather than a single gene, as single-gene analysis may miss important effects on pathways and difficult to reproduce from studies to studies [1]. Researchers have made significant progress in identifying metabolic or signaling pathways based on expression array data [2, 3]. Meanwhile, new tools for identification of pathways, such as GenMAPP [4], Pathway Processor [5], MAPPFinder [6], have made pathway data more widely available. However, It is a challenging task to model the pathway data and test for a potentially complex pathway effect on a disease outcome.

One way to model pathway data is through the linear model approach, where the pathway effect is represented by a linear combination of individual gene effects. This approach has several limitations. Activities of genes within a pathway are often complicated, thus a linear model is often insufficient to capture the relationship between these genes. Furthermore, genes within a pathway tend to interact with each other. Such interactions are not taken into account by the linear model approach.

In this paper we propose a nonparametric approach, the kernel machine regression, to model a pathway effect. The kernel machine method, with the support vector machine (SVM) as a most popular example, has emerged in the last decade as a powerful machine learning technique in high-dimensional settings [7, 8]. This method provides a flexible way to model linear and nonlinear effects of variables and gene-gene interactions, unifies the model building procedure in both one- and multi-dimensional settings, and shows attractive performance compared to other nonparametric methods such as splines.

Liu et al. [9] proposed a kernel machine-based regression model for continuous outcomes. In this paper, we propose a logistic kernel machine regression model for binary outcomes, where covariate effects are modeled parametrically and the genetic pathway effect is modeled parametrically or nonparametrically using the kernel machine method. A main contribution of this paper is to establish a connection between logistic kernel machine regression and the logistic mixed model. We show that the kernel machine estimator of the genetic pathway effect can be obtained from the estimator of the random effects in the corresponding logistic mixed model. This connection provides a convenient vehicle to connect the powerful kernel machine method with the popular mixed model method in the statistical literature. This mixed model connection also provides an unified framework for statistical inference for model parameters, including the regression coefficients, the nonparametric genetic pathway function, and the regularization and kernel parameters. Based on the proposed logistic kernel machine regression model, we develop a new test for the nonlinear pathway effect on disease risk. An appealing feature of the proposed test is that it performs well without the need to correctly specify the functional form of the effects of each gene or their interactions. This feature has a significant practical implication when analyzing genetic pathway data, as the true relationship between the pathway and the disease outcome is often unknown. We extend the results to generalized kernel machine regression for a class of continuous and discrete outcomes and discuss its connection with generalized linear mixed models [10].

Recently, Wei and Li [11] proposed a nonparametric pathway-based regression (NPR) to model pathway data. NPR is a pathway-based gradient boosting procedure, where the base learner is usually a regression or classification tree. It provides a flexible approach in modeling pathways and interactions among genes within a pathway. Michalowski et al. [12] proposed a Bayesian Belief Network approach for pathway data. Neither method is likelihood-based. Thus parameter estimation and inference cannot be casted within a unified likelihood framework. It is hence difficult to estimate and quantify the overall pathway effect on disease risk and assess its statistical uncertainty. Secondly, a primary interest in this paper is to test for the statistical significance of the overall pathway effect on the risk of a disease. Both NPR and Bayesian belief network do not provide such a statistical test for the pathway effect. For example, NPR uses an importance score to rank the relative importance of each pathway. It lacks formal inferential procedure for assessing the statistical significance of a pathway. Further, when considering a single pathway, the importance score loses its meaning in assessing the importance of a pathway. Our method, on the other hand, is based on penalized likelihood, and estimation and inference can be conducted in a systematic manner within the likelihood framework. We also propose a formal statistical test for the significance of a pathway effect on the risk of a disease.

Goeman et al. [13] proposed a linear mixed model to relate the pathway effect with a continuous outcome. They modeled the pathway effect using a linear function with each gene entering into the model as a regressor. They assumed the regression coefficients of the gene as random from a common distribution with mean 0 and an unknown variance. The pathway effect can then be tested through a variance component test for random effects. Our approach is different from theirs in the following aspects. First, we model the pathway effect by allowing for a nonparametric model rather than a parametric one. As we commented earlier, the highly complicated nature of activities of genes within a pathway makes the linear model assumption untenable. Secondly, the kernel function used in kernel machine regression usually contains unknown tuning parameters. The parameter is present under the alternative hypothesis but disappears under null hypothesis. This makes tests as proposed in [13, 14] not applicable. Our proposed test, on the other hand, works quite well under this scenario. Third, Goeman et al. [14] extended their linear model results to discrete outcomes using basis functions. A key advantage of the kernel machine approach over this basis approach for modeling multi-gene effects is that one does not need to specify bases explicitly, which is often difficult for high-dimensional data especially when interactions are modeled.

Results

Analysis of prostate cancer data

In this section, we apply the proposed logistic kernel machine regression model (3) as described in the Methods section to the analysis of a prostate cancer data set. The data came from the Michigan prostate cancer study [15]. This study involved 81 patients with 22 diagnosed as non-cancerous and 59 diagnosed with local or advanced prostate cancer. Besides the clinical and demographic covariates such as age, cDNA microarray gene expressions were also available for each patient. The early results of Dhanasekaran et al. [15] indicate that certain functional genetic pathways seemed dys-regulated in prostate cancer relative to non-cancerous tissues. We are interested in studying how a genetic pathway is related to the prostate cancer risk, controlling for the covariates. We focus in this analysis on the cell growth pathway, which contains 5 genes. The pathway we describe was annotated by the investigator (A. Chinnaiyan) and is simply used to illustrate the methodology. Of course, one could take the pathways stored in commercial databases such as Ingenuity Pathway Analysis (IPA) and use the proposed methodology based on those gene sets.

The outcome was the binary prostate cancer status and the covariate includes age. Since the functional relationship between the cell growth pathway and the prostate cancer risk is unknown, the kernel machine method provides a convenient and flexible framework for the evaluation of the pathway effect on the prostate cancer risk. Specifically, we consider the following semiparametric logistic modellogit(P(y = 1)) = _β_0 + _β_1age + h(gene1, ..., gene5),

where h(·) is a nonparametric function of 5 genes within the cell growth pathway. The detail of the estimation procedure is provided in the Methods section. In summary, we fit this model using the kernel machine method via the logistic mixed model representation and using the Gaussian kernel function in estimating h(·). Under the mixed model representation, we estimated (_β_0, _β_1) and h(·) using penalized quasi-likelihood (PQL), and estimated the smoothing parameter τ and the Gaussian kernel scale parameter ρ simultaneously by treating them as variance components. The results are presented in Table 1.

Table 1 Analysis of prostate cancer data.

Full size table

The test for the cell growth pathway effect on the prostate cancer status _H_0: h(z) = 0 vs _H_1:h(z) ≠ 0 was conducted using the proposed score test as described in the Methods section. For the purpose of comparison, we also conducted the global test proposed by Goeman et al. [13] that assumed a linear pathway effect. Note that our test allows a nonlinear pathway effect and gene-gene interactions. Table 1 gives the p-values for both tests. The p-value of our test suggests that cell growth pathway has a highly significant effect on the disease status, while the test from Goeman et al. [13] indicates only marginal significance of the growth pathway effect.

Simulation Study for the Parameter Estimates

We conducted a simulation study to evaluate the performance of the parameter estimates of the proposed logistic kernel machine regression by using the logistic mixed model formulation. We considered the following modellogit(P(y i = 1)) = x i + h(z i 1, ⋯, z ip ),

where the true regression coefficient β = 1. We consider p = 5 and set h(_z_1, ..., _z_5) = 2{sin(_z_1) – z 2 2 MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemOEaO3aa0baaSqaaiabikdaYaqaaiabikdaYaaaaaa@2F63@ + _z_1 exp(-_z_3) -sin(_z_2) cos(_z_3) + z 4 2 MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemOEaO3aa0baaSqaaiabisda0aqaaiabikdaYaaaaaa@2F67@ + sin(_z_4) cos(_z_1) + z 5 2 MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemOEaO3aa0baaSqaaiabiwda1aqaaiabikdaYaaaaaa@2F69@ + z_3_z_5}. To allow x i and (z i 1, ⋯, z ip ) to be correlated, x i was generated as x i = sin(z i 1) + 2_u i , where u i and z ij (j = 1, ⋯, p) follow independent Uniform(-0.5, 0.5). The Gaussian kernel was used throughout the simulation. All simulations ran 300 times. Settings 1, 2, and 3 correspond to sample size n = 100, 200, and 300, respectively.

The simulation results are shown in Table 2. Due to the multi-dimensional nature of the variables z, it is difficult to visualize the fitted curve h ^ MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGafmiAaGMbaKaaaaa@2D3E@ (z). We hence summarized the goodness-of-fit of h ^ MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGafmiAaGMbaKaaaaa@2D3E@ (·) in the following way. For each simulated data set, we regressed the true h on the fitted value h ^ MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGafmiAaGMbaKaaaaa@2D3E@ , both evaluated at the design points. We then empirically summarized the goodness-of-fit of h ^ MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGafmiAaGMbaKaaaaa@2D3E@ (·) by calculating the average intercepts, slopes and _R_2's obtained from these regressions over the 300 simulations. If the kernel machine method fits the nonparametric function well, then we would expect the intercept to be close to 0, the slope close to 1, and _R_2also close to 1.

Table 2 Simulation results on estimation.

Full size table

Our results show that even when the sample size is as low as 100, estimation of the regression coefficient and nonparametric function only has small bias. When the kernel parameter ρ is estimated, these biases tend to be small compared with those when ρ is held fixed. With the increase of sample size the estimates of β and h become closer to the true values, especially when ρ is estimated, while there are still some bias when ρ is fixed at values farther away from the estimated one. Table 3 compares the estimated standard errors of β ^ MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGafqOSdiMbaKaaaaa@2D86@ with the empirical standard errors. Our results show that they agree to each other well when ρ is estimated.

Table 3 Simulation results on standard errors.

Full size table

Simulation Study of the Score Test for the Pathway Effect

We next conducted a simulation study to evaluate the performance of the proposed variance component score test for the pathway effect _H_0: h(·) = 0 vs _H_1: h(·) ≠ 0. In order to compare the performance of our test with the linearity-based global test proposed by Goeman et al. [13], both tests were applied to each simulated data set. Nonlinear and linear functions of h(z) were both considered. For the nonlinear pathway effect, the true model is logit(y) = x +ah(z), where h(z) = 2(_z_1 - _z_2)2 + _z_2_z_3 + 3 sin(2_z_3)_z_4 + z 5 2 MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemOEaO3aa0baaSqaaiabiwda1aqaaiabikdaYaaaaaa@2F69@ + 2 cos(_z_4)_z_5. For the linear pathway effect, the true model is logit(y) = x + ah(z), where h(z) = 2_z_1 + 3_z_2+_z_3 + 2_z_4+_z_5. All z's were generated from the standard normal distribution, and a = 0, 0.2, 0.4, 0.6, 0.8. To allow x and (z i 1, ⋯, z ip ) to be correlated, x was generated as x = _z_1+e/2 with e being independent of _z_1 and following N (0, 1). We studied the size of the test by generating data under a = 0, and studied the power by increasing a. The sample size was 100. For the size calculations, the number of simulations was 2000; whereas for the power calculations, the number of runs was 1000. Based on the discussions in Section "Test for the genetic Pathway Effect", the bound of ρ is set up by interval [ min i ≠ j ∑ l = 1 5 ( z i l − z j l ) 2 / 5 , 10 max i ≠ j ∑ l = 1 5 ( z i l − z j l ) 2 ] MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaei4waSLagiyBa0MaeiyAaKMaeiOBa42aaSbaaSqaaiabdMgaPjabgcMi5kabdQgaQbqabaGcdaaeWaqaaiabcIcaOiabdQha6naaBaaaleaacqWGPbqAcqWGSbaBaeqaaOGaeyOeI0IaemOEaO3aaSbaaSqaaiabdQgaQjabdYgaSbqabaGccqGGPaqkdaahaaWcbeqaaiabikdaYaaakiabc+caViabiwda1iabcYcaSiabigdaXiabicdaWiGbc2gaTjabcggaHjabcIha4naaBaaaleaacqWGPbqAcqGHGjsUcqWGQbGAaeqaaaqaaiabdYgaSjabg2da9iabigdaXaqaaiabiwda1aqdcqGHris5aOWaaabmaeaacqGGOaakcqWG6bGEdaWgaaWcbaGaemyAaKMaemiBaWgabeaakiabgkHiTiabdQha6naaBaaaleaacqWGQbGAcqWGSbaBaeqaaOGaeiykaKYaaWbaaSqabeaacqaIYaGmaaaabaGaemiBaWMaeyypa0JaeGymaedabaGaeGynaudaniabggHiLdGccqGGDbqxaaa@6A7F@ , and the interval is divided by 500 equally spaced grid points. All simulations were conducted using R 2.5.0, and the package "globaltest" v4.6.0 was used for the test proposed by Goeman et al. [13] as a comparison.

Table 4 reports the empirical size (a = 0) and power (a > 0) of the variance component score test for the pathway effect. When the true function h(z) is non-linear in z, the results show that the size of our test was very close to the nominal value 0.05, while the size of the global test of Goeman et al. [13] is inflated. The results also show that our test had a much higher power. This was not surprising since the test of Goeman et al. [13] was based on a linearity assumption of the pathway effect. When the true underlying model is far from linear, the linearity assumption breaks down and the test quickly loses power. The results also show that the proposed test works well for moderate sample sizes. When the pathway effect is linear, the results show that the size of both tests were very close to the nominal value 0.05 and their power were also very close. This demonstrates that our test is as powerful as the global test when the true underlying h(z) is linear. Therefore our test could be used as a universal test for testing the overall effect of a set of variables without the need to specify the true functional forms of each variable. This feature is especially desirable for genetic pathway data, because the relationship between genes and clinical outcome is often unknown.

Table 4 Simulation results on score test.

Full size table

Conclusions and Discussion

In this paper, we developed a logistic kernel machine regression model for binary outcomes, where the covariate effects are modeled parametrically and the genetic pathway effect is modeled nonparametrically using the kernel machine method. This method provides an attractive way to model the pathway effect, without the need to make strong parametric assumptions on individual gene effects or their interactions. Our model also allows for parametric pathway effects if a parametric kernel, such as the first-degree polynomial kernel, is used.

A key result of this paper is that we have established a close connection between the generalized kernel machine regression and generalized linear mixed models, and show that the kernel machine estimators of regression coefficients and the nonparametric multi-dimensional pathway effect can be easily obtained from the corresponding generalized linear mixed models using PQL. The mixed model connection provides a unified framework for estimation and inference and can be easily implemented in existing software, such as SAS PROC GLIMMIX or R GLMMPQL. The mixed model connection also makes it possible to test for the overall pathway effect through the proposed variance component test. A key advantage of the proposed score test for the pathway effect is that it does not require an explicit functional specification of individual gene effects and gene-gene interactions. This feature is of practical significance as the pathway effect is often complex. Our simulation study shows the proposed test performs well for moderate sample size. It has similar power to the linearity-based pathway test of Goeman et al. [13] when the true effect is linear, but much higher power when the true effect is nonlinear.

We have considered in this paper a single pathway. One could generalize the proposed semiparametric model to incorporate multiple pathways by fitting an additive model:logit(P(y = 1)) = xT β+ _h_1(z1) + ⋯ + h m (z m ),

where z j (j = 1, ⋯, m) denotes a p j × 1 vector of genes in the j th pathway and hj(·) denotes the nonparametric function associated with the j th genetic pathway.

Machine learning is a powerful tool in advancing bioinformatics research. Our effort helps to build a bridge between kernel machine methods and traditional statistical models. This connection will undoubtedly provide a new and convenient tool for the bioinformatics community and opens a door for future research.

Methods

The Logistic Kernel Machine Model

Throughout the paper we assume that gene expression data have been properly normalized. Suppose the data consist of n samples. For subject i (i = 1, ⋯, n), y i is a binary disease outcome taking values either 0 (non-disease) or 1 (disease), x i is a q × 1 vector of covariates, z i is a p × 1 vector of gene expression measurements in a pathway/gene set. We assume that an intercept is included in x i . The binary outcome y i depends on x i and z i through the following semiparametric logistic regression model:

logit ( μ i ) = x i T β + h ( z i ) , MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeeiBaWMaee4Ba8Maee4zaCMaeeyAaKMaeeiDaqNaeiikaGIaeqiVd02aaSbaaSqaaiabdMgaPbqabaGccqGGPaqkcqGH9aqpieWacqWF4baEdaqhaaWcbaGaemyAaKgabaGaemivaqfaaGGadOGae4NSdiMaey4kaSIaemiAaGMaeiikaGIae8NEaO3aaSbaaSqaaiabdMgaPbqabaGccqGGPaqkcqGGSaalaaa@46C6@

(3)

where μ i = P (y i = 1| x i , z i ), β is a q × 1 vector of regression coefficients, and h(z i ) is an unknown centered smooth function.

In model (3), covariate effects are modeled parametrically, while the multi-dimensional genetic pathway effect is modeled parametrically or nonparametrically. A nonparametric specification for h(·) reflects our limited knowledge of genetic functional forms. Note that h(·) = 0 means that genes in the pathway have no association with the disease risk. If h(z) = _γ_1_z_1 + ⋯ + γ p z p , model (3) becomes the linear model considered by Goeman et al. [13].

In nonparametric modeling, such as smoothing splines, the unknown function is usually assumed to lie in a certain function space. For the kernel machine method, this function space, denoted by H K MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaWenfgDOvwBHrxAJfwnHbqeg0uy0HwzTfgDPnwy1aaceaGae83cHG0aaSbaaSqaaiabdUealbqabaaaaa@37D6@ , is generated by a given positive definite kernel function K(·, ·). The mathematical properties of H K MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaWenfgDOvwBHrxAJfwnHbqeg0uy0HwzTfgDPnwy1aaceaGae83cHG0aaSbaaSqaaiabdUealbqabaaaaa@37D6@ imply that any unknown function h(z) in H K MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaWenfgDOvwBHrxAJfwnHbqeg0uy0HwzTfgDPnwy1aaceaGae83cHG0aaSbaaSqaaiabdUealbqabaaaaa@37D6@ can be written as a linear combination of the given kernel function K(·, ·) evaluated at each sample point. Two popular kernel functions are the d th polynomial kernel K ( z 1 , z 2 ) = ( z 1 T z 2 + ρ ) d MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaem4saSKaeiikaGccbmGae8NEaO3aaSbaaSqaaiabigdaXaqabaGccqGGSaalcqWF6bGEdaWgaaWcbaGaeGOmaidabeaakiabcMcaPiabg2da9iabcIcaOiab=Pha6naaDaaaleaacqaIXaqmaeaacqWGubavaaGccqWF6bGEdaWgaaWcbaGaeGOmaidabeaakiabgUcaRiabeg8aYjabcMcaPmaaCaaaleqabaGaemizaqgaaaaa@421C@ and the Gaussian Kernel K(**z**1, **z**2) = exp{-|| **z**1 – z2||2/_ρ_2}, where | | z 1 − z 2 | | 2 = ∑ k = 1 p ( z 1 k − z 2 k ) 2 MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeiiFaWNaeiiFaWhcbmGae8NEaO3aaSbaaSqaaiabigdaXaqabaGccqGHsislcqWF6bGEdaWgaaWcbaGaeGOmaidabeaakiabcYha8jabcYha8naaCaaaleqabaGaeGOmaidaaOGaeyypa0ZaaabmaeaacqGGOaakcqWG6bGEdaWgaaWcbaGaeGymaeJaem4AaSgabeaakiabgkHiTiabdQha6naaBaaaleaacqaIYaGmcqWGRbWAaeqaaOGaeiykaKYaaWbaaSqabeaacqaIYaGmaaaabaGaem4AaSMaeyypa0JaeGymaedabaGaemiCaahaniabggHiLdaaaa@4CB6@ and ρ is an unknown parameter. The first and second degree polynomial kernels (d = 1, 2) correspond to assuming h(·) to be linear and quadratic in z's, respectively. The choice of a kernel function determines which function space one would like to use to approximate h(z). The unknown parameter of a kernel function plays a critical role in function approximation. It is a challenging problem to optimally estimate it from data. In the machine learning literature, this parameter is usually pre-fixed at some values based on some ad-hoc methods. In this paper, we show that we can optimally estimate it from data based on a mixed model framework.

The Estimation Procedure

Assuming h(·) ∈ H K MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaWenfgDOvwBHrxAJfwnHbqeg0uy0HwzTfgDPnwy1aaceaGae83cHG0aaSbaaSqaaiabdUealbqabaaaaa@37D6@ , the function space generated by a kernel function K(·, ·), we can estimate β and h(·) by maximizing the penalized log-likelihood function

J ( h , β ) = ∑ i = 1 n { y i log ( μ i 1 − μ i ) + log ( 1 − μ i ) } − 1 2 λ | | h | | H K 2 = ∑ i = 1 n ( y i { x i T β + h ( z i ) } − log [ 1 + exp { x i T β + h ( z i ) } ] ) − λ 2 | | h | | H K 2 , MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaqbaeaabiWaaaqaaiabdQeakjabcIcaOiabdIgaOjabcYcaSGGadiab=j7aIjabcMcaPaqaaiabg2da9aqaamaaqahabaWaaiWaaeaacqWG5bqEdaWgaaWcbaGaemyAaKgabeaakiGbcYgaSjabc+gaVjabcEgaNnaabmaajuaGbaWaaSaaaeaacqaH8oqBdaWgaaqaaiabdMgaPbqabaaabaGaeGymaeJaeyOeI0IaeqiVd02aaSbaaeaacqWGPbqAaeqaaaaaaOGaayjkaiaawMcaaiabgUcaRiGbcYgaSjabc+gaVjabcEgaNjabcIcaOiabigdaXiabgkHiTiabeY7aTnaaBaaaleaacqWGPbqAaeqaaOGaeiykaKcacaGL7bGaayzFaaaaleaacqWGPbqAcqGH9aqpcqaIXaqmaeaacqWGUbGBa0GaeyyeIuoakiabgkHiTKqbaoaalaaabaGaeGymaedabaGaeGOmaidaaOGaeq4UdWMaeiiFaWNaeiiFaWNaemiAaGMaeiiFaWNaeiiFaW3aa0baaSqaamrtHrhAL1wy0L2yHvtyaeHbnfgDOvwBHrxAJfwnaGabaiab+TqiinaaBaaameaacqWGlbWsaeqaaaWcbaGaeGOmaidaaaGcbaaabaGaeyypa0dabaWaaabCaeaacqGGOaakcqWG5bqEdaWgaaWcbaGaemyAaKgabeaakiabcUha7Hqadiab9Hha4naaDaaaleaacqWGPbqAaeaacqWGubavaaGccqWFYoGycqGHRaWkcqWGObaAcqGGOaakcqqF6bGEdaWgaaWcbaGaemyAaKgabeaakiabcMcaPiabc2ha9jabgkHiTiGbcYgaSjabc+gaVjabcEgaNjabcUfaBjabigdaXiabgUcaRiGbcwgaLjabcIha4jabcchaWjabcUha7jab9Hha4naaDaaaleaacqWGPbqAaeaacqWGubavaaGccqWFYoGycqGHRaWkcqWGObaAcqGGOaakcqqF6bGEdaWgaaWcbaGaemyAaKgabeaakiabcMcaPiabc2ha9jabc2faDjabcMcaPiabgkHiTKqbaoaalaaabaGaeq4UdWgabaGaeGOmaidaaOGaeiiFaWNaeiiFaWNaemiAaGMaeiiFaWNaeiiFaW3aa0baaSqaaiab+TqiinaaBaaameaacqWGlbWsaeqaaaWcbaGaeGOmaidaaaqaaiabdMgaPjabg2da9iabigdaXaqaaiabd6gaUbqdcqGHris5aOGaeiilaWcaaaaa@BE12@

(4)

where λ is a regularization parameter that controls the tradeoff between goodness of fit and complexity of the model. When λ = 0, it fits a saturated model, and when λ = ∞, the model reduces to a simple logistic model logit (μ i ) = x i T β MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaacbmGae8hEaG3aa0baaSqaaiabdMgaPbqaaiabdsfaubaaiiWakiab+j7aIbaa@31C1@ . Note that there are two tuning parameters in the above likelihood function, the regularization parameter λ and kernel parameter ρ Intuitively, λ controls the magnitude of the unknown function while ρ mainly governs the smoothness property of the function.

By the representer theorem [16], the general solution for the nonparametric function h(·) in (4) can be expressed as

h ( z i ) = ∑ i ′ = 1 n α i ′ K ( z i , z i ′ ) = k i T α , MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemiAaGMaeiikaGccbmGae8NEaO3aaSbaaSqaaiabdMgaPbqabaGccqGGPaqkcqGH9aqpdaaeWbqaaiabeg7aHnaaBaaaleaacuWGPbqAgaqbaaqabaGccqWGlbWscqGGOaakcqWF6bGEdaWgaaWcbaGaemyAaKgabeaakiabcYcaSiab=Pha6naaBaaaleaacuWGPbqAgaqbaaqabaGccqGGPaqkaSqaaiqbdMgaPzaafaGaeyypa0JaeGymaedabaGaemOBa4ganiabggHiLdGccqGH9aqpcqWFRbWAdaqhaaWcbaGaemyAaKgabaGaemivaqfaaGGadOGae4xSdeMaeiilaWcaaa@4F0F@

(5)

where k i = {K(z i , z1), ..., K(z i , z n )}_T_and **α**= (_α_1, ⋯, α n )T, an n × 1 vector of unknown parameters.

Substituting (5) into (4) we have

J ( β , α ) = ∑ i = 1 n [ y i ( x i T β + k i T α ) − log { 1 + exp ( x i T β + k i T α ) } ] − 1 2 λ α T K α , MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemOsaOKaeiikaGcccmGae8NSdiMaeiilaWIae8xSdeMaeiykaKIaeyypa0ZaaabCaeaadaWadaqaaiabdMha5naaBaaaleaacqWGPbqAaeqaaOGaeiikaGccbmGae4hEaG3aa0baaSqaaiabdMgaPbqaaiabdsfaubaakiab=j7aIjabgUcaRiab+TgaRnaaDaaaleaacqWGPbqAaeaacqWGubavaaGccqWFXoqycqGGPaqkcqGHsislcyGGSbaBcqGGVbWBcqGGNbWzdaGadaqaaiabigdaXiabgUcaRiGbcwgaLjabcIha4jabcchaWnaabmaabaGae4hEaG3aa0baaSqaaiabdMgaPbqaaiabdsfaubaakiab=j7aIjabgUcaRiab+TgaRnaaDaaaleaacqWGPbqAaeaacqWGubavaaGccqWFXoqyaiaawIcacaGLPaaaaiaawUhacaGL9baaaiaawUfacaGLDbaaaSqaaiabdMgaPjabg2da9iabigdaXaqaaiabd6gaUbqdcqGHris5aOGaeyOeI0scfa4aaSaaaeaacqaIXaqmaeaacqaIYaGmaaGccqaH7oaBcqWFXoqydaahaaWcbeqaaiabdsfaubaakiab+Tealjab=f7aHjabcYcaSaaa@753D@

(6)

where K= K(ρ) is an n × n matrix whose (i, i')th element is K(z i , z i ') and often depends on an unknown parameter ρ.

Since J(β, α) in (6) is a nonlinear function of (β, α), one can use the Fisher scoring or Newton-Raphson iterative algorithm to maximize (6) with respect to β and α. Let (k) denote the _k_th iteration step, then it can be shown (for details see Appendix A.3) that the (k + 1)th update for β and α solves the following normal equation:

[ X T D ( k ) X X T D ( k ) K D ( k ) X τ − 1 I + D ( k ) K ] [ β ( k + 1 ) α ( k + 1 ) ] = [ X T D ( k ) y ˜ ( k ) D ( k ) y ˜ ( k ) ] . MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaWaamWaaeaafaqabeGacaaabaacbmGae8hwaG1aaWbaaSqabeaacqWGubavaaGccqWFebardaahaaWcbeqaaiabcIcaOiabdUgaRjabcMcaPaaakiab=Hfaybqaaiab=HfaynaaCaaaleqabaGaemivaqfaaOGae8hraq0aaWbaaSqabeaacqGGOaakcqWGRbWAcqGGPaqkaaGccqWFlbWsaeaacqWFebardaahaaWcbeqaaiabcIcaOiabdUgaRjabcMcaPaaakiab=Hfaybqaaiabes8a0naaCaaaleqabaGaeyOeI0IaeGymaedaaOGae8xsaKKaey4kaSIae8hraq0aaWbaaSqabeaacqGGOaakcqWGRbWAcqGGPaqkaaGccqWFlbWsaaaacaGLBbGaayzxaaWaamWaaeaafaqabeGabaaabaaccmGae4NSdi2aaWbaaSqabeaacqGGOaakcqWGRbWAcqGHRaWkcqaIXaqmcqGGPaqkaaaakeaacqGFXoqydaahaaWcbeqaaiabcIcaOiabdUgaRjabgUcaRiabigdaXiabcMcaPaaaaaaakiaawUfacaGLDbaacqGH9aqpdaWadaqaauaabeqaceaaaeaacqWFybawdaahaaWcbeqaaiabdsfaubaakiab=reaenaaCaaaleqabaGaeiikaGIaem4AaSMaeiykaKcaaOGaf8xEaKNbaGaadaahaaWcbeqaaiabcIcaOiabdUgaRjabcMcaPaaaaOqaaiab=reaenaaCaaaleqabaGaeiikaGIaem4AaSMaeiykaKcaaOGaf8xEaKNbaGaadaahaaWcbeqaaiabcIcaOiabdUgaRjabcMcaPaaaaaaakiaawUfacaGLDbaacqGGUaGlaaa@7771@

(7)

where y ˜ ( k ) = X β ( k ) + K α ( k ) + D ( k ) − 1 ( y − μ ( k ) ) MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaacbmGaf8xEaKNbaGaadaahaaWcbeqaaiabcIcaOiabdUgaRjabcMcaPaaakiabg2da9iab=HfayHGadiab+j7aInaaCaaaleqabaGaeiikaGIaem4AaSMaeiykaKcaaOGaey4kaSIae83saSKae4xSde2aaWbaaSqabeaacqGGOaakcqWGRbWAcqGGPaqkaaGccqGHRaWkcqWFebardaahaaWcbeqaaiabcIcaOiabdUgaRjabcMcaPmaaCaaameqabaGaeyOeI0IaeGymaedaaaaakiabcIcaOiab=Lha5jabgkHiTiab+X7aTnaaCaaaleqabaGaeiikaGIaem4AaSMaeiykaKcaaOGaeiykaKcaaa@4F08@ , τ = 1/λ, h(k)= Kα(k), and D ( k ) = Diag { μ i ( k ) ( 1 − μ i ( k ) ) } MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaacbmGae8hraq0aaWbaaSqabeaacqGGOaakcqWGRbWAcqGGPaqkaaGccqGH9aqpcqqGebarcqqGPbqAcqqGHbqycqqGNbWzcqGG7bWEcqaH8oqBdaqhaaWcbaGaemyAaKgabaGaeiikaGIaem4AaSMaeiykaKcaaOGaeiikaGIaeGymaeJaeyOeI0IaeqiVd02aa0baaSqaaiabdMgaPbqaaiabcIcaOiabdUgaRjabcMcaPaaakiabcMcaPiabc2ha9baa@4983@ . The estimators β ^ MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaaccmGaf8NSdiMbaKaaaaa@2D8E@ and h ^ MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaacbmGaf8hAaGMbaKaaaaa@2D46@ at convergence are the kernel machine estimators that maximize (6).

The Connection of Logistic Kernel Machine Regression to Logistic Mixed Models

Generalized linear mixed models (GLMMs) have been used to analyze correlated categorical data and have gained much popularity in the statistical literature [10]. Logistic mixed models are a special case of GLMMs. We show in this section that the kernel machine estimator in the semiparametric logistic regression model (3) corresponds to the Penalized Quasi-Likelihood (PQL) [10] estimator from a logistic mixed model, and the regularization parameter τ = 1/λ and kernel parameter ρ can be treated as variance components and estimated simultaneously from the corresponding logistic mixed model. Specifically, consider the following logistic mixed model:

logit ( μ i ) = x i T β + h i , MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeeiBaWMaee4Ba8Maee4zaCMaeeyAaKMaeeiDaqNaeiikaGIaeqiVd02aaSbaaSqaaiabdMgaPbqabaGccqGGPaqkcqGH9aqpieWacqWF4baEdaqhaaWcbaGaemyAaKgabaGaemivaqfaaGGadOGae4NSdiMaey4kaSIaemiAaG2aaSbaaSqaaiabdMgaPbqabaGccqGGSaalaaa@439B@

(8)

where β is a q × 1 vector of fixed effects, and **h**= (_h_1, ..., h n ) is a n × 1 vector of subject-specific random effects following h ~N{0, τ K(ρ)}, and the covariance matrix K(ρ) is the n × n kernel matrix as defined in previous section.

As K is not diagonal or block-diagonal, the random effects h i 's across all subjects are correlated. The _i_th mean response μ i depends on other random effects h i ' (i' ≠ i) through the correlations of h i with other random effects. To estimate the unknown parameters in the logistic mixed model (8), we estimate β and h by maximizing the PQL [10], which can be viewed as a joint log likelihood of (β, h),

∑ i = 1 n [ y i ( x i T β + h i ) − log { 1 + exp ( x i T β + h i ) } ] − 1 2 τ h T K − 1 h . MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaWaaabCaeaadaWadaqaaiabdMha5naaBaaaleaacqWGPbqAaeqaaOGaeiikaGccbmGae8hEaG3aa0baaSqaaiabdMgaPbqaaiabdsfaubaaiiWakiab+j7aIjabgUcaRiabdIgaOnaaBaaaleaacqWGPbqAaeqaaOGaeiykaKIaeyOeI0IagiiBaWMaei4Ba8Maei4zaCMaei4EaSNaeGymaeJaey4kaSIagiyzauMaeiiEaGNaeiiCaaNaeiikaGIae8hEaG3aa0baaSqaaiabdMgaPbqaaiabdsfaubaakiab+j7aIjabgUcaRiabdIgaOnaaBaaaleaacqWGPbqAaeqaaOGaeiykaKIaeiyFa0hacaGLBbGaayzxaaGaeyOeI0scfa4aaSaaaeaacqaIXaqmaeaacqaIYaGmcqaHepaDaaGccqWFObaAdaahaaWcbeqaaiabdsfaubaakiab=TealnaaCaaaleqabaGaeyOeI0IaeGymaedaaOGae8hAaGMaeiOla4caleaacqWGPbqAcqGH9aqpcqaIXaqmaeaacqWGUbGBa0GaeyyeIuoaaaa@6A46@

(9)

Setting τ = 1/λ and **h**= Kα, one can easily see that equations (6) and (9) are identical. It follows that the logistic kernel machine estimators β ^ MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaaccmGaf8NSdiMbaKaaaaa@2D8E@ and h ^ MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaacbmGaf8hAaGMbaKaaaaa@2D46@ can be obtained by fitting the logistic mixed model (8) using PQL. In fact, examination of the kernel machine normal equations (7) shows that they are identical to the normal equations obtained from the PQL (9) [10], where y ˜ MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaacbmGaf8xEaKNbaGaaaaa@2D67@ in (7) is in fact the PQL working vector and D is the PQL working weight matrix.

Note that the estimators of β and h depend on the unknown regularization parameter τ and the kernel parameter ρ. Within the PQL framework, we can estimate these parameters **δ**= (τ, ρ) by maximizing the approximate REML likelihood

ℓ R ( β ^ ( δ ) , δ ) ≈ − 1 2 log | V | − 1 2 log | X T V − 1 X | − 1 2 ( y ˜ − X β ^ ) T V − 1 ( y ˜ − X β ^ ) , MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeS4eHW2aaSbaaSqaaiabdkfasbqabaGccqGGOaakcuaHYoGygaqcaiabcIcaOGGadiab=r7aKjabcMcaPiabcYcaSiab=r7aKjabcMcaPiabgIKi7kabgkHiTKqbaoaalaaabaGaeGymaedabaGaeGOmaidaaOGagiiBaWMaei4Ba8Maei4zaCMaeiiFaWhcbmGae4NvayLaeiiFaWNaeyOeI0scfa4aaSaaaeaacqaIXaqmaeaacqaIYaGmaaGccyGGSbaBcqGGVbWBcqGGNbWzcqGG8baFcqGFybawdaahaaWcbeqaaiabdsfaubaakiab+zfawnaaCaaaleqabaGaeyOeI0IaeGymaedaaOGae4hwaGLaeiiFaWNaeyOeI0scfa4aaSaaaeaacqaIXaqmaeaacqaIYaGmaaGccqGGOaakcuGF5bqEgaacaiabgkHiTiab+Hfayjqb=j7aIzaajaGaeiykaKYaaWbaaSqabeaacqWGubavaaGccqGFwbGvdaahaaWcbeqaaiabgkHiTiabigdaXaaakiabcIcaOiqb+Lha5zaaiaGaeyOeI0Iae4hwaGLaf8NSdiMbaKaacqGGPaqkcqGGSaalaaa@6E2F@

(10)

where **V**= **D**-1 + τ K, and y ˜ MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaacbmGaf8xEaKNbaGaaaaa@2D67@ is the working vector as defined above. The estimator δ ^ MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaaccmGaf8hTdqMbaKaaaaa@2D92@ of δ can be obtained by setting equal to zero the first derivative of (10) with respect to δ. The estimating procedure for β, h, and δ= (τ, ρ) can be summarized as follows: we fit the logistic kernel machine model by iteratively fitting the following working linear mixed model to estimate (β, h) using BLUPs and to estimate δ using REML, until convergence

y ˜ = X β + h + ε , MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaacbmGaf8xEaKNbaGaacqGH9aqpcqWFybawiiWacqGFYoGycqGHRaWkcqWFObaAcqGHRaWkcqGF1oqzcqGGSaalaaa@3732@

where y ˜ MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaacbmGaf8xEaKNbaGaaaaa@2D67@ is the working vector defined below equation (7), h is a random effect vector following N{0, τ K(ρ)}, and ε ~ N(0, D). The covariance of β ^ MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaaccmGaf8NSdiMbaKaaaaa@2D8E@ is estimated by (XT V-1X)-1, and the covariance of h ^ MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaacbmGaf8hAaGMbaKaaaaa@2D46@ is estimated by τ ^ K − τ ^ K P K MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGafqiXdqNbaKaaieWacqWFlbWscqGHsislcuaHepaDgaqcaiab=Tealjab=bfaqjab=Tealbaa@34EE@ , where P= V-1 - V-1X(XT V-1X)-1XT **V**-1 and **V**= V( δ ^ MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaaccmGaf8hTdqMbaKaaaaa@2D92@ ). The covariance of δ ^ MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaaccmGaf8hTdqMbaKaaaaa@2D92@ can be obtained as the inverse of the expected information matrix calculated using the second derivative of (10) with respect to δ. The square roots of the diagonal elements of the estimated covariance matrices give the standard errors of β ^ MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaaccmGaf8NSdiMbaKaaaaa@2D8E@ , h ^ MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaacbmGaf8hAaGMbaKaaaaa@2D46@ , and δ. The above discussion shows that we can easily fit the logistic kernel machine regression using the existing PQL-based mixed model software, such as SAS GLIMMIX and R GLMMPQL.

Test for the Genetic Pathway Effect

It is of significant practical interest to test the overall genetic pathway effect _H_0: h(z) = 0. Assuming h(z) ∈ H K MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaWenfgDOvwBHrxAJfwnHbqeg0uy0HwzTfgDPnwy1aaceaGae83cHG0aaSbaaSqaaiabdUealbqabaaaaa@37D6@ , one can easily see from the logistic mixed model representation (8) that _H_0: h(z) = 0 vs _H_1: h(z) ≠ 0 is equivalent to testing the variance component τ as _H_0: τ = 0 vs _H_1: τ > 0. Note that the null hypothesis places τ on the boundary of the parameter space. Since the kernel matrix K is not block diagonal, unlike the standard case considered by Self and Liang [17], the likelihood ratio for _H_0: τ = 0 does not follow a mixture of χ 0 2 MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeq4Xdm2aa0baaSqaaiabicdaWaqaaiabikdaYaaaaaa@2F99@ and χ 1 2 MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeq4Xdm2aa0baaSqaaiabigdaXaqaaiabikdaYaaaaaa@2F9B@ distribution. We consider instead a score test in this paper.

When conducting statistical tests for pathways, two types of tests could be formulated. The first is called the competitive test and the second the self-contained test [18]. The competitive test compares an interested gene set to all the other genes on a gene chip. An example of the competitive test is the gene set enrichment analysis (GSEA) [1], where an enrichment score of a gene set is defined and a permutation test is used to test for the significance of the gene set based on the enrichment score. The self-contained test compares the gene set to an internal standard which does not involve any genes outside the gene set considered. In other words, the self-contained test examines the null hypothesis that a pathway has no effect on the outcome versus the alternative hypothesis that the pathway has an effect. The variance component test of [13] for the linear pathway effect is a self-contained test. Goeman and Bühlmann [18] pointed out that the self-contained test has a higher power than a competitive test and that its statistical formulation is also consistent for both single gene tests and gene set tests, and the statistical sampling properties of the competitive test can be difficult to interpret.

Our pathway effect hypothesis _H_0: h(z) = 0 vs _H_1: h(z) ≠ 0 is a self-contained hypothesis. We propose in this paper a self-contained test for the pathway effect by developing a kernel machine variance component score test for _H_0: τ = 0 vs _H_0:τ > 0. The proposed test allows for both linear and nonlinear pathway effects and includes the tests by Goeman et al. [13, 14] as a special case. A key advantage of our kernel-based test is that we do not need to explicitly specify the basis functions for h(·), which is often difficult for modeling the joint effects of multiple genes, and we all let the data to estimate the best curvature of h(·).

Zhang and Lin [19] proposed a score test for _H_0: τ = 0 to compare a polynomial model with a smoothing spline. Goeman et al. [14] also proposed a global test against a high dimensional alternative under the empirical Bayesian framework. The variance-covariance matrix used in these tests do not involve any unknown parameters. However, the kernel function K(·, ·) in a kernel machine model usually depends on some unknown parameter ρ. One can easily see from the mixed model representation (8) that under _H_0: τ = 0, the kernel matrix K disappears. This makes the parameter ρ inestimable under the null hypothesis and therefore renders the above tests inapplicable.

Davies [20, 21] studied the problem of a parameter disappearing under _H_0 and proposed a score test by treating the score statistic as a Gaussian process indexed by the nuisance parameter and then obtaining an upper bound to approximate the p-value of the score test. We adopt this line of approaches for our proposed score test.

Using the derivative of (10) with respect to τ, we propose the following score test statistic for _H_0: τ = 0 as

S ( ρ ) = Q τ ( β 0 ∧ , ρ ) − μ Q σ Q , MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaem4uamLaeiikaGIaeqyWdiNaeiykaKIaeyypa0tcfa4aaSaaaeaacqWGrbqudaWgaaqaaiabes8a0bqabaGaeiikaGYaaecaaeaacqaHYoGydaWgaaqaaiabicdaWaqabaaacaGLcmaacqGGSaalcqaHbpGCcqGGPaqkcqGHsislcqaH8oqBdaWgaaqaaiabdgfarbqabaaabaGaeq4Wdm3aaSbaaeaacqWGrbquaeqaaaaakiabcYcaSaaa@4525@

(11)

where

Q τ ( β ^ 0 , ρ ) = ( y ˜ − X β ^ 0 ) T D K ( ρ ) D ( y ˜ − X β ^ 0 ) = ( y − μ ^ 0 ) T K ( y − μ ^ 0 ) , MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemyuae1aaSbaaSqaaiabes8a0bqabaGccqGGOaakdaqiaaqaaGGadiab=j7aIbGaayPadaWaaSbaaSqaaiabicdaWaqabaGccqGGSaalcqaHbpGCcqGGPaqkcqGH9aqpcqGGOaakieWacuGF5bqEgaacaiabgkHiTiab+HfaynaaHaaabaGae8NSdigacaGLcmaadaWgaaWcbaGaeGimaadabeaakiabcMcaPmaaCaaaleqabaGaemivaqfaaOGae4hraqKae43saSKaeiikaGIaeqyWdiNaeiykaKIae4hraqKaeiikaGIaf4xEaKNbaGaacqGHsislcqGFybawdaqiaaqaaiab=j7aIbGaayPadaWaaSbaaSqaaiabicdaWaqabaGccqGGPaqkcqGH9aqpcqGGOaakcqGF5bqEcqGHsislcuWF8oqBgaqcamaaBaaaleaacqaIWaamaeqaaOGaeiykaKYaaWbaaSqabeaacqWGubavaaGccqGFlbWscqGGOaakcqGF5bqEcqGHsislcuWF8oqBgaqcamaaBaaaleaacqaIWaamaeqaaOGaeiykaKIaeiilaWcaaa@6458@

where β ^ 0 MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaWaaecaaeaaiiWacqWFYoGyaiaawkWaamaaBaaaleaacqaIWaamaeqaaaaa@2F5A@ is the MLE of β under _H_0: τ = 0, μ ^ 0 = l o g i t − 1 ( X β ^ 0 ) MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaaccmGaf8hVd0MbaKaadaWgaaWcbaGaeGimaadabeaakiabg2da9iabdYgaSjabd+gaVjabdEgaNjabdMgaPjabdsha0naaCaaaleqabaGaeyOeI0IaeGymaedaaOGaeiikaGccbmGae4hwaGLaf8NSdiMbaKaadaWgaaWcbaGaeGimaadabeaakiabcMcaPaaa@3E8E@ , μ Q = tr{P0K(ρ)}, σ Q 2 MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeq4Wdm3aa0baaSqaaiabdgfarbqaaiabikdaYaaaaaa@2FE2@ = 2tr{P0K(ρ)P0K(ρ)}, and **P**0 = D0 – D0X(XT D0X)-1 XT **D**0, where D 0 = diag { μ ^ i 0 ( 1 − μ i 0 ) } MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaacbmGae8hraq0aaSbaaSqaaiabicdaWaqabaGccqGH9aqpcqqGKbazcqqGPbqAcqqGHbqycqqGNbWzcqGG7bWEcuaH8oqBgaqcamaaBaaaleaacqWGPbqAcqaIWaamaeqaaOGaeiikaGIaeGymaeJaeyOeI0IaeqiVd02aaSbaaSqaaiabdMgaPjabicdaWaqabaGccqGGPaqkcqGG9bqFaaa@4367@ . Note that under _H_0: τ = 0, model (3) reduces to the simple logistic model logit (μ i ) = x i T β MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaacbmGae8hEaG3aa0baaSqaaiabdMgaPbqaaiabdsfaubaaiiWakiab+j7aIbaa@31C1@ . Hence the β ^ 0 MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaWaaecaaeaaiiWacqWFYoGyaiaawkWaamaaBaaaleaacqaIWaamaeqaaaaa@2F5A@ is the MLE of β under this null logistic model.

If the Gaussian kernel is used, then an arbitrary nonlinear pathway effect is implicitly assumed. Our proposed test, which is derived to test for any nonlinear effect, is therefore more powerful than tests based on a parametric assumption. We show in Appendix A.1 that when ρ becomes large in the Gaussian kernel, our test statistic reduces asymptotically to the one based on linearity assumption of genetic effects. Hence our test includes linear model based test as a special case. From (11) it is also clear that our test is invariant to the relative scaling of the kernel function K(·, ·).

Under appropriate regularity conditions similar to those specified in [22], S(ρ) under the null hypothesis can be considered as an approximate Gaussian process indexed by ρ. Using this formulation, we can then apply Davies' results [20, 21] to obtain the upper bound for the p-value of the test. Since a large value of _Q_τ ( β ^ MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaaccmGaf8NSdiMbaKaaaaa@2D8E@ , ρ) would lead to the rejection of _H_0, the p-value of the test corresponds to the up-crossing probability. Following Davies [21], the p-value is upper-bounded by

Φ ( − M ) + W exp ( − 1 2 M 2 ) / 8 π , MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeuOPdyKaeiikaGIaeyOeI0Iaemyta0KaeiykaKIaey4kaSIaem4vaCLagiyzauMaeiiEaGNaeiiCaaNaeiikaGIaeyOeI0scfa4aaSaaaeaacqaIXaqmaeaacqaIYaGmaaGccqWGnbqtdaahaaWcbeqaaiabikdaYaaakiabcMcaPiabc+caVmaakaaabaGaeGioaGJaeqiWdahaleqaaOGaeiilaWcaaa@43C7@

(12)

where Φ (·) is the normal cumulative distribution function, M is the maximum of S(ρ) over the range of ρ, W = | S(_ρ_1) – S(L)| + | S(_ρ_2) – S(_ρ_1) | + ... + | S(U) – S(ρ m ) |, L and U are the lower and upper bound of ρ respectively and ρ l , l = 1, ..., m are the m grid points between L and U. Davies [20] points out that this bound is sharp. For the Gaussian kernel, we suggest to set the bound of ρ as L = 0.1 min i ≠ j ∑ l = 1 p ( z i l − z j l ) 2 MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGagiyBa0MaeiyAaKMaeiOBa42aaSbaaSqaaiabdMgaPjabgcMi5kabdQgaQbqabaGcdaaeWaqaaiabcIcaOiabdQha6naaBaaaleaacqWGPbqAcqWGSbaBaeqaaOGaeyOeI0IaemOEaO3aaSbaaSqaaiabdQgaQjabdYgaSbqabaGccqGGPaqkdaahaaWcbeqaaiabikdaYaaaaeaacqWGSbaBcqGH9aqpcqaIXaqmaeaacqWGWbaCa0GaeyyeIuoaaaa@4801@ and U = 100 max i ≠ j ∑ l = 1 p ( z i l − z j l ) 2 MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGagiyBa0MaeiyyaeMaeiiEaG3aaSbaaSqaaiabdMgaPjabgcMi5kabdQgaQbqabaGcdaaeWaqaaiabcIcaOiabdQha6naaBaaaleaacqWGPbqAcqWGSbaBaeqaaOGaeyOeI0IaemOEaO3aaSbaaSqaaiabdQgaQjabdYgaSbqabaGccqGGPaqkdaahaaWcbeqaaiabikdaYaaaaeaacqWGSbaBcqGH9aqpcqaIXaqmaeaacqWGWbaCa0GaeyyeIuoaaaa@4805@ . For justifications, see Appendix A.2.

Extension to generalized kernel machine model

For simplicity, we focus in this paper on logistic regression for binary outcomes. The proposed semiparametric model (3) can be easily extended to other types of continuous and discrete outcomes, such as normal, count, skewed data, whose distributions are in the exponential family [23]. In this section, we briefly discuss how to generalize our estimation and testing procedures for binary outcomes to other data types within the generalized kernel machine framework and discuss its fitting using generalized linear mixed models.

Suppose the data consist of n independent subjects. For subject i (i = 1, ..., n), y i is a response variable, x i is a q × 1 vector of covariates, z i is a p × 1 vector of gene expressions within a pathway. Suppose y i follows a distribution in the exponential family with density [23]

p ( y i ; θ i , φ ) = exp { y i θ i − a ( θ i ) φ / m i + c ( y i , φ ) } , MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemiCaaNaeiikaGIaemyEaK3aaSbaaSqaaiabdMgaPbqabaGccqGG7aWocqaH4oqCdaWgaaWcbaGaemyAaKgabeaakiabcYcaSiabeA8aMjabcMcaPiabg2da9iGbcwgaLjabcIha4jabcchaWnaacmaabaqcfa4aaSaaaeaacqWG5bqEdaWgaaqaaiabdMgaPbqabaGaeqiUde3aaSbaaeaacqWGPbqAaeqaaiabgkHiTiabdggaHjabcIcaOiabeI7aXnaaBaaabaGaemyAaKgabeaacqGGPaqkaeaacqaHgpGzcqGGVaWlcqWGTbqBdaWgaaqaaiabdMgaPbqabaaaaOGaey4kaSIaem4yamMaeiikaGIaemyEaK3aaSbaaSqaaiabdMgaPbqabaGccqGGSaalcqaHgpGzcqGGPaqkaiaawUhacaGL9baacqGGSaalaaa@5E66@

(13)

where θ i is the canonical parameter, a(·) and c(·) are known functions, φ is a dispersion parameter, and m i is a known weight. The mean of y i satisfies μ i = E(y i ) = a'(θ i ) and Var(y i ) = φ m i a"(θ i ). The generalized kernel machine model is an extension of the generalized linear model [23] by allowing the pathway effect to be modeled nonparametrically using kernel machine as

g ( μ i ) = x i T β + h ( z i ) , MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaem4zaCMaeiikaGIaeqiVd02aaSbaaSqaaiabdMgaPbqabaGccqGGPaqkcqGH9aqpieWacqWF4baEdaqhaaWcbaGaemyAaKgabaGaemivaqfaaGGadOGae4NSdiMaey4kaSIaemiAaGMaeiikaGIae8NEaO3aaSbaaSqaaiabdMgaPbqabaGccqGGPaqkcqGGSaalaaa@413C@

(14)

where g(·) is a known monotone link function, and h(·) is an unknown centered smooth function lying in the function space H K MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaWenfgDOvwBHrxAJfwnHbqeg0uy0HwzTfgDPnwy1aaceaGae83cHG0aaSbaaSqaaiabdUealbqabaaaaa@37D6@ generated by a positive definite kernel function K(·, ·). For binary data, setting g(μ) = logit(μ) = log μ 1 − μ MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGagiiBaWMaei4Ba8Maei4zaCwcfa4aaSaaaeaacqaH8oqBaeaacqaIXaqmcqGHsislcqaH8oqBaaaaaa@35DA@ gives the logistic kernel machine model (3); for count data, g(μ) = log(μ) gives the Poisson kernel machine model; for Gaussian data, g(μ) = μ gives linear kernel machine model [9]. The regression coefficients β and the nonparametric function h(·) in (14) can be obtained by maximizing the penalized log-likelihood function

J ( h , β ) = ∑ i = 1 n ℓ { y i , x i , z i ; β , h ( ⋅ ) } − 1 2 λ | | h | | H K 2 MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemOsaOKaeiikaGIaemiAaGMaeiilaWcccmGae8NSdiMaeiykaKIaeyypa0ZaaabCaeaacqWItecBcqGG7bWEcqWG5bqEdaWgaaWcbaGaemyAaKgabeaakiabcYcaSGqadiab+Hha4naaBaaaleaacqWGPbqAaeqaaOGaeiilaWIae4NEaO3aaSbaaSqaaiabdMgaPbqabaGccqGG7aWocqWFYoGycqGGSaalcqWGObaAcqGGOaakcqGHflY1cqGGPaqkcqGG9bqFaSqaaiabdMgaPjabg2da9iabigdaXaqaaiabd6gaUbqdcqGHris5aOGaeyOeI0scfa4aaSaaaeaacqaIXaqmaeaacqaIYaGmaaGccqaH7oaBcqGG8baFcqGG8baFcqWGObaAcqGG8baFcqGG8baFdaqhaaWcbaWenfgDOvwBHrxAJfwnHbqeg0uy0HwzTfgDPnwy1aaceaGae03cHG0aaSbaaWqaaiabdUealbqabaaaleaacqaIYaGmaaaaaa@6C70@

(15)

where ℓ(·) = ln(p) is the log-likelihood, p is the density function given in (13), and λ is a tuning parameter. Using the kernel expression of h(·) in (5), the generalized kernel machine model (14) can be written as

g ( μ i ) = x i T β + k i T α , MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaem4zaCMaeiikaGIaeqiVd02aaSbaaSqaaiabdMgaPbqabaGccqGGPaqkcqGH9aqpieWacqWF4baEdaqhaaWcbaGaemyAaKgabaGaemivaqfaaGGadOGae4NSdiMaey4kaSIae83AaS2aa0baaSqaaiabdMgaPbqaaiabdsfaubaakiab+f7aHjabcYcaSaaa@40DE@

and the penalized likelihood can be written

J ( β , α ) = ∑ i = 1 n ℓ ( y i , x i , z i ; β , α ) − 1 2 λ α T K α , MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemOsaOKaeiikaGcccmGae8NSdiMaeiilaWIae8xSdeMaeiykaKIaeyypa0ZaaabCaeaacqWItecBcqGGOaakcqWG5bqEdaWgaaWcbaGaemyAaKgabeaakiabcYcaSGqadiab+Hha4naaBaaaleaacqWGPbqAaeqaaOGaeiilaWIae4NEaO3aaSbaaSqaaiabdMgaPbqabaGccqGG7aWocqWFYoGycqGGSaalcqWFXoqycqGGPaqkaSqaaiabdMgaPjabg2da9iabigdaXaqaaiabd6gaUbqdcqGHris5aOGaeyOeI0scfa4aaSaaaeaacqaIXaqmaeaacqaIYaGmaaGccqaH7oaBcqWFXoqydaahaaWcbeqaaiabdsfaubaakiab+Tealjab=f7aHjabcYcaSaaa@59BB@

(16)

where K is an n × n matrix whose (i, j)th element is K(z i , z j ).

One can use the Fisher scoring iteration to solve for β and a. The procedure is virtually the same as that described in Section "The Estimation Procedure". The normal equation takes the same form as (7), except that now μ i is specified under (14) and **D**= diag{var(y i )} under (13). Similar calculations to those in Section "The Connection of Logistic Kernel Machine Regression to Logistic Mixed Models" show that model (14) can be fit using the generalized linear mixed model [10] via PQL

g ( μ i b ) = x i T β + h i , MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaem4zaCMaeiikaGIaeqiVd02aa0baaSqaaiabdMgaPbqaaGqadiab=jgaIbaakiabcMcaPiabg2da9iab=Hha4naaDaaaleaacqWGPbqAaeaacqWGubavaaaccmGccqGFYoGycqGHRaWkcqWGObaAdaWgaaWcbaGaemyAaKgabeaakiabcYcaSaaa@3F5B@

where τ = 1/λ, and **h**= (_h_1 ..., h n ) is an n × n random vector with distribution N {0, τ K(ρ)}. The same PQL statistical software, such as SAS PROC GLIMMIX and R GLMMPQL, can be used to fit this model and obtain the kernel machine estimators of β and h(·).

The score test (11) also has a straightforward extension. The only change is that the elements in matrix D in (11) be replaced by appropriate variance function var(y i ) under the assumed parametric distribution of y i .

Appendix

A.1 Proof of the relationship of the proposed score test and that of Goeman, et al [13] under the linearity assumption

We show in this section when the scale parameter is large, the proposed nonparametric variance component test for the pathway effect using the Gaussian kernel reduces to the linearity-based global test of Goeman et al. [13].

Suppose K(·) is the Gaussian kernel. It can be shown that the score statistic for testing _H_0: τ = 0 satisfies

Q τ ( β ^ , ρ ) = ( y ˜ − X β ^ 0 ) T D K ( ρ ) D ( y ˜ − X β ^ 0 ) = ( y − μ ^ ) T K ( ρ ) ( y − μ ^ 0 ) , MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemyuae1aaSbaaSqaaiabes8a0bqabaGccqGGOaakdaqiaaqaaGGadiab=j7aIbGaayPadaGaeiilaWIaeqyWdiNaeiykaKIaeyypa0JaeiikaGccbmGaf4xEaKNbaGaacqGHsislcqGFybawdaqiaaqaaiab=j7aIbGaayPadaWaaSbaaSqaaiabicdaWaqabaGccqGGPaqkdaahaaWcbeqaaiabdsfaubaakiabhseaejab+TealjabcIcaOiabeg8aYjabcMcaPiabhseaejabcIcaOiqb+Lha5zaaiaGaeyOeI0Iae4hwaG1aaecaaeaacqWFYoGyaiaawkWaamaaBaaaleaacqaIWaamaeqaaOGaeiykaKIaeyypa0JaeiikaGIae4xEaKNaeyOeI0Iaf8hVd0MbaKaacqGGPaqkdaahaaWcbeqaaiabdsfaubaakiab+TealjabcIcaOiabeg8aYjabcMcaPiabcIcaOiab+Lha5jabgkHiTiqb=X7aTzaajaWaaSbaaSqaaiabicdaWaqabaGccqGGPaqkcqGGSaalaaa@6594@

(17)

where μ ^ 0 MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaaccmGaf8hVd0MbaKaadaWgaaWcbaGaeGimaadabeaaaaa@2EBD@ is the MLE of μ under _H_0. The test statistic of Goeman et al. [13] takes the form

( y − μ ^ ) T R ( h − μ ^ ) , MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeiikaGccbmGae8xEaKNaeyOeI0cccmGaf4hVd0MbaKaacqGGPaqkdaahaaWcbeqaaiabdsfaubaakiab=jfasjabcIcaOiab=HgaOjabgkHiTiqb+X7aTzaajaGaeiykaKIaeiilaWcaaa@3B37@

(18)

where **R**= ZZT. We now show when ρ is large relative to max i ≠ j ∑ l = 1 p ( z i l − z j l ) 2 MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGagiyBa0MaeiyyaeMaeiiEaG3aaSbaaSqaaiabdMgaPjabgcMi5kabdQgaQbqabaGcdaaeWaqaaiabcIcaOiabdQha6naaBaaaleaacqWGPbqAcqWGSbaBaeqaaOGaeyOeI0IaemOEaO3aaSbaaSqaaiabdQgaQjabdYgaSbqabaGccqGGPaqkdaahaaWcbeqaaiabikdaYaaaaeaacqWGSbaBcqGH9aqpcqaIXaqmaeaacqWGWbaCa0GaeyyeIuoaaaa@4805@

ρ 2 ( y − μ ^ ) T K ( ρ ) ( y − μ ^ ) ≈ ( y − μ ^ ) T R ( y − μ ^ ) . MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaqcfa4aaSaaaeaacqaHbpGCaeaacqaIYaGmaaGccqGGOaakieWacqWF5bqEcqGHsisliiWacuGF8oqBgaqcaiabcMcaPmaaCaaaleqabaGaemivaqfaaOGae83saSKaeiikaGIaeqyWdiNaeiykaKIaeiikaGIae8xEaKNaeyOeI0Iaf4hVd0MbaKaacqGGPaqkcqGHijYUcqGGOaakcqWF5bqEcqGHsislcuGF8oqBgaqcaiabcMcaPmaaCaaaleqabaGaemivaqfaaOGae8NuaiLaeiikaGIae8xEaKNaeyOeI0Iaf4hVd0MbaKaacqGGPaqkcqGGUaGlaaa@5209@

(19)

Simple Taylor expansions show that

( y − μ ^ ) T K ( ρ ) ( y − μ ^ ) = ∑ i = 1 n ∑ j = 1 n ( y i − μ ^ i ) ( y j − μ ^ j ) exp { − ∑ l = 1 p ( z i l − z j l ) 2 / ρ } = ∑ i = 1 n ( y i − μ ^ i ) 2 + ∑ i ≠ j ( y i − μ ^ i ) ( y j − μ ^ j ) exp { − ∑ l = 1 p ( z i l − z j l ) 2 / ρ } . MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaqbaeaabiWaaaqaaiabcIcaOGqadiab=Lha5jabgkHiTGGadiqb+X7aTzaajaGaeiykaKYaaWbaaSqabeaacqWGubavaaGccqWFlbWscqGGOaakcqaHbpGCcqGGPaqkcqGGOaakcqWF5bqEcqGHsislcuGF8oqBgaqcaiabcMcaPaqaaiabg2da9aqaamaaqahabaWaaabCaeaacqGGOaakcqWG5bqEdaWgaaWcbaGaemyAaKgabeaakiabgkHiTiqbeY7aTzaajaWaaSbaaSqaaiabdMgaPbqabaGccqGGPaqkcqGGOaakcqWG5bqEdaWgaaWcbaGaemOAaOgabeaakiabgkHiTiqbeY7aTzaajaWaaSbaaSqaaiabdQgaQbqabaGccqGGPaqkcyGGLbqzcqGG4baEcqGGWbaCcqGG7bWEcqGHsisldaaeWbqaaiabcIcaOiabdQha6naaBaaaleaacqWGPbqAcqWGSbaBaeqaaOGaeyOeI0IaemOEaO3aaSbaaSqaaiabdQgaQjabdYgaSbqabaGccqGGPaqkdaahaaWcbeqaaiabikdaYaaakiabc+caViabeg8aYjabc2ha9bWcbaGaemiBaWMaeyypa0JaeGymaedabaGaemiCaahaniabggHiLdaaleaacqWGQbGAcqGH9aqpcqaIXaqmaeaacqWGUbGBa0GaeyyeIuoaaSqaaiabdMgaPjabg2da9iabigdaXaqaaiabd6gaUbqdcqGHris5aaGcbaaabaGaeyypa0dabaWaaabCaeaacqGGOaakcqWG5bqEdaWgaaWcbaGaemyAaKgabeaakiabgkHiTiqbeY7aTzaajaWaaSbaaSqaaiabdMgaPbqabaGccqGGPaqkdaahaaWcbeqaaiabikdaYaaakiabgUcaRmaaqafabaGaeiikaGIaemyEaK3aaSbaaSqaaiabdMgaPbqabaGccqGHsislcuaH8oqBgaqcamaaBaaaleaacqWGPbqAaeqaaOGaeiykaKIaeiikaGIaemyEaK3aaSbaaSqaaiabdQgaQbqabaGccqGHsislcuaH8oqBgaqcamaaBaaaleaacqWGQbGAaeqaaOGaeiykaKIagiyzauMaeiiEaGNaeiiCaaNaei4EaSNaeyOeI0YaaabCaeaacqGGOaakcqWG6bGEdaWgaaWcbaGaemyAaKMaemiBaWgabeaakiabgkHiTiabdQha6naaBaaaleaacqWGQbGAcqWGSbaBaeqaaOGaeiykaKYaaWbaaSqabeaacqaIYaGmaaGccqGGVaWlcqaHbpGCcqGG9bqFaSqaaiabdYgaSjabg2da9iabigdaXaqaaiabdchaWbqdcqGHris5aOGaeiOla4caleaacqWGPbqAcqGHGjsUcqWGQbGAaeqaniabggHiLdaaleaacqWGPbqAcqGH9aqpcqaIXaqmaeaacqWGUbGBa0GaeyyeIuoaaaaaaa@C87E@

When max i ≠ j ∑ l = 1 p ( z i l − z j l ) 2 / ρ MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGagiyBa0MaeiyyaeMaeiiEaG3aaSbaaSqaaiabdMgaPjabgcMi5kabdQgaQbqabaGcdaaeWaqaaiabcIcaOiabdQha6naaBaaaleaacqWGPbqAcqWGSbaBaeqaaOGaeyOeI0IaemOEaO3aaSbaaSqaaiabdQgaQjabdYgaSbqabaGccqGGPaqkdaahaaWcbeqaaiabikdaYaaaaeaacqWGSbaBcqGH9aqpcqaIXaqmaeaacqWGWbaCa0GaeyyeIuoakiabc+caViabeg8aYbaa@4AB5@ is small, i.e., when ρ is large relative to max i ≠ j ∑ l = 1 p ( z i l − z j l ) 2 MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGagiyBa0MaeiyyaeMaeiiEaG3aaSbaaSqaaiabdMgaPjabgcMi5kabdQgaQbqabaGcdaaeWaqaaiabcIcaOiabdQha6naaBaaaleaacqWGPbqAcqWGSbaBaeqaaOGaeyOeI0IaemOEaO3aaSbaaSqaaiabdQgaQjabdYgaSbqabaGccqGGPaqkdaahaaWcbeqaaiabikdaYaaaaeaacqWGSbaBcqGH9aqpcqaIXaqmaeaacqWGWbaCa0GaeyyeIuoaaaa@4805@ , we have that exp { − ∑ l = 1 p ( z i l − z j l ) 2 / ρ } ≈ 1 − ∑ l = 1 p ( z i l − z j l ) 2 / ρ MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGagiyzauMaeiiEaGNaeiiCaaNaei4EaSNaeyOeI0YaaabmaeaacqGGOaakcqWG6bGEdaWgaaWcbaGaemyAaKMaemiBaWgabeaakiabgkHiTiabdQha6naaBaaaleaacqWGQbGAcqWGSbaBaeqaaOGaeiykaKYaaWbaaSqabeaacqaIYaGmaaaabaGaemiBaWMaeyypa0JaeGymaedabaGaemiCaahaniabggHiLdGccqGGVaWlcqaHbpGCcqGG9bqFcqGHijYUcqaIXaqmcqGHsisldaaeWaqaaiabcIcaOiabdQha6naaBaaaleaacqWGPbqAcqWGSbaBaeqaaOGaeyOeI0IaemOEaO3aaSbaaSqaaiabdQgaQjabdYgaSbqabaGccqGGPaqkdaahaaWcbeqaaiabikdaYaaaaeaacqWGSbaBcqGH9aqpcqaIXaqmaeaacqWGWbaCa0GaeyyeIuoakiabc+caViabeg8aYbaa@638E@ for any i ≠ j. Hence

( y − μ ^ ) T K ( ρ ) ( y − μ ^ ) = ∑ i = 1 n ( y i − μ ^ i ) 2 + ∑ i ≠ j ( y i − μ ^ i ) ( y j − μ ^ j ) − 1 ρ ∑ i ≠ j ( y i − μ ^ i ) ( y j − μ ^ j ) ∑ l = 1 p ( z i l − z j l ) 2 . MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaqbaeaabiGaaaqaaaqaaiabcIcaOGqadiab=Lha5jabgkHiTGGadiqb+X7aTzaajaGaeiykaKYaaWbaaSqabeaacqWGubavaaGccqWFlbWscqGGOaakcqaHbpGCcqGGPaqkcqGGOaakcqWF5bqEcqGHsislcuGF8oqBgaqcaiabcMcaPaqaaiabg2da9aqaamaaqahabaGaeiikaGIaemyEaK3aaSbaaSqaaiabdMgaPbqabaGccqGHsislcuaH8oqBgaqcamaaBaaaleaacqWGPbqAaeqaaOGaeiykaKYaaWbaaSqabeaacqaIYaGmaaGccqGHRaWkdaaeqbqaaiabcIcaOiabdMha5naaBaaaleaacqWGPbqAaeqaaOGaeyOeI0IafqiVd0MbaKaadaWgaaWcbaGaemyAaKgabeaakiabcMcaPiabcIcaOiabdMha5naaBaaaleaacqWGQbGAaeqaaOGaeyOeI0IafqiVd0MbaKaadaWgaaWcbaGaemOAaOgabeaakiabcMcaPaWcbaGaemyAaKMaeyiyIKRaemOAaOgabeqdcqGHris5aaWcbaGaemyAaKMaeyypa0JaeGymaedabaGaemOBa4ganiabggHiLdGccqGHsisljuaGdaWcaaqaaiabigdaXaqaaiabeg8aYbaakmaaqafabaGaeiikaGIaemyEaK3aaSbaaSqaaiabdMgaPbqabaGccqGHsislcuaH8oqBgaqcamaaBaaaleaacqWGPbqAaeqaaOGaeiykaKIaeiikaGIaemyEaK3aaSbaaSqaaiabdQgaQbqabaGccqGHsislcuaH8oqBgaqcamaaBaaaleaacqWGQbGAaeqaaOGaeiykaKcaleaacqWGPbqAcqGHGjsUcqWGQbGAaeqaniabggHiLdGcdaaeWbqaaiabcIcaOiabdQha6naaBaaaleaacqWGPbqAcqWGSbaBaeqaaOGaeyOeI0IaemOEaO3aaSbaaSqaaiabdQgaQjabdYgaSbqabaGccqGGPaqkdaahaaWcbeqaaiabikdaYaaaaeaacqWGSbaBcqGH9aqpcqaIXaqmaeaacqWGWbaCa0GaeyyeIuoakiabc6caUaaaaaa@9B27@

Since ∑ j = 1 ( y j − μ ^ j ) = 0 MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaWaaabeaeaacqGGOaakcqWG5bqEdaWgaaWcbaGaemOAaOgabeaakiabgkHiTiqbeY7aTzaajaWaaSbaaSqaaiabdQgaQbqabaGccqGGPaqkcqGH9aqpcqaIWaamaSqaaiabdQgaQjabg2da9iabigdaXaqab0GaeyyeIuoaaaa@3C05@ under the PQL, we have ∑ j ≠ 1 ( y j − μ ^ j ) = − ( y i − μ ^ i ) MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaWaaabeaeaacqGGOaakcqWG5bqEdaWgaaWcbaGaemOAaOgabeaakiabgkHiTiqbeY7aTzaajaWaaSbaaSqaaiabdQgaQbqabaGccqGGPaqkcqGH9aqpcqGHsislcqGGOaakcqWG5bqEdaWgaaWcbaGaemyAaKgabeaakiabgkHiTiqbeY7aTzaajaWaaSbaaSqaaiabdMgaPbqabaGccqGGPaqkaSqaaiabdQgaQjabgcMi5kabigdaXaqab0GaeyyeIuoaaaa@45C7@ . Hence

( y − μ ^ ) T K ( ρ ) ( y − μ ^ ) ≈ ∑ i = 1 n ( y i − μ ^ i ) 2 − ∑ i = j n ( y i − μ ^ i ) 2 − 1 ρ ∑ i ≠ j ( y i − μ ^ i ) ( y j − μ ^ j ) ∑ l = 1 p ( z i l 2 − 2 z i l z j l + z j l 2 ) = 2 ρ ∑ i = 1 n ( y i − μ ^ i ) 2 ∑ l = 1 p z i l 2 + 2 ρ ∑ i ≠ j ( y i − μ ^ i ) ( y j − μ ^ j ) ∑ l = 1 p z i l z j l = 2 ρ ( y − μ ^ ) T R ( y − μ ^ ) . MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaqbaeaabqGaaaaabaaabaGaeiikaGccbmGae8xEaKNaeyOeI0cccmGaf4hVd0MbaKaacqGGPaqkdaahaaWcbeqaaiabdsfaubaakiab=TealjabcIcaOiabeg8aYjabcMcaPiabcIcaOiab=Lha5jabgkHiTiqb+X7aTzaajaGaeiykaKcabaGaeyisISlabaWaaabCaeaacqGGOaakcqWG5bqEdaWgaaWcbaGaemyAaKgabeaakiabgkHiTiqbeY7aTzaajaWaaSbaaSqaaiabdMgaPbqabaGccqGGPaqkdaahaaWcbeqaaiabikdaYaaakiabgkHiTmaaqahabaGaeiikaGIaemyEaK3aaSbaaSqaaiabdMgaPbqabaGccqGHsislcuaH8oqBgaqcamaaBaaaleaacqWGPbqAaeqaaOGaeiykaKYaaWbaaSqabeaacqaIYaGmaaaabaGaemyAaKMaeyypa0JaeGymaedabaGaemOBa4ganiabggHiLdaaleaacqWGPbqAcqGH9aqpcqaIXaqmaeaacqWGUbGBa0GaeyyeIuoakiabgkHiTKqbaoaalaaabaGaeGymaedabaGaeqyWdihaaOWaaabuaeaacqGGOaakcqWG5bqEdaWgaaWcbaGaemyAaKgabeaakiabgkHiTiqbeY7aTzaajaWaaSbaaSqaaiabdMgaPbqabaGccqGGPaqkcqGGOaakcqWG5bqEdaWgaaWcbaGaemOAaOgabeaakiabgkHiTiqbeY7aTzaajaWaaSbaaSqaaiabdQgaQbqabaGccqGGPaqkaSqaaiabdMgaPjabgcMi5kabdQgaQbqab0GaeyyeIuoakmaaqahabaGaeiikaGIaemOEaO3aa0baaSqaaiabdMgaPjabdYgaSbqaaiabikdaYaaakiabgkHiTiabikdaYiabdQha6naaBaaaleaacqWGPbqAcqWGSbaBaeqaaOGaemOEaO3aaSbaaSqaaiabdQgaQjabdYgaSbqabaGccqGHRaWkcqWG6bGEdaqhaaWcbaGaemOAaOMaemiBaWgabaGaeGOmaidaaOGaeiykaKcaleaacqWGSbaBcqGH9aqpcqaIXaqmaeaacqWGWbaCa0GaeyyeIuoaaOqaaiabg2da9aqaaKqbaoaalaaabaGaeGOmaidabaGaeqyWdihaaOWaaabCaeaacqGGOaakcqWG5bqEdaWgaaWcbaGaemyAaKgabeaakiabgkHiTiqbeY7aTzaajaWaaSbaaSqaaiabdMgaPbqabaGccqGGPaqkdaahaaWcbeqaaiabikdaYaaakmaaqahabaGaemOEaO3aa0baaSqaaiabdMgaPjabdYgaSbqaaiabikdaYaaaaeaacqWGSbaBcqGH9aqpcqaIXaqmaeaacqWGWbaCa0GaeyyeIuoakiabgUcaRKqbaoaalaaabaGaeGOmaidabaGaeqyWdihaaOWaaabuaeaacqGGOaakcqWG5bqEdaWgaaWcbaGaemyAaKgabeaakiabgkHiTiqbeY7aTzaajaWaaSbaaSqaaiabdMgaPbqabaGccqGGPaqkcqGGOaakcqWG5bqEdaWgaaWcbaGaemOAaOgabeaakiabgkHiTiqbeY7aTzaajaWaaSbaaSqaaiabdQgaQbqabaGccqGGPaqkaSqaaiabdMgaPjabgcMi5kabdQgaQbqab0GaeyyeIuoaaSqaaiabdMgaPjabg2da9iabigdaXaqaaiabd6gaUbqdcqGHris5aOWaaabCaeaacqWG6bGEdaWgaaWcbaGaemyAaKMaemiBaWgabeaakiabdQha6naaBaaaleaacqWGQbGAcqWGSbaBaeqaaaqaaiabdYgaSjabg2da9iabigdaXaqaaiabdchaWbqdcqGHris5aaGcbaGaeyypa0dabaqcfa4aaSaaaeaacqaIYaGmaeaacqaHbpGCaaGccqGGOaakcqWF5bqEcqGHsislcuGF8oqBgaqcaiabcMcaPmaaCaaaleqabaGaemivaqfaaOGae8NuaiLaeiikaGIae8xEaKNaeyOeI0Iaf4hVd0MbaKaacqGGPaqkcqGGUaGlaaaaaa@010D@

This proves the approximate relation (19).

A.2 Calculations of the lower and upper bounds of ρ

Although in theory ρ could take any positive values up to infinity, for computational purpose we would require ρ to be bounded. For the proposed test statistic (11), its value in fact only depends on a finite range of ρ values. We describe why this is the case and how to find this range. For a given data set, the proof in Appendix A.1 shows that when is sufficiently large, the quantity 0.5_ρQ_ τ ( β ^ 0 MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaWaaecaaeaaiiWacqWFYoGyaiaawkWaamaaBaaaleaacqaIWaamaeqaaaaa@2F5A@ , ρ) converges to S 0 = ( y ˜ − μ ^ 0 ) T R ( y ˜ − μ ^ 0 ) MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaem4uam1aaSbaaSqaaiabicdaWaqabaGccqGH9aqpcqGGOaakieWacuWF5bqEgaacaiabgkHiTGGadiqb+X7aTzaajaWaaSbaaSqaaiabicdaWaqabaGccqGGPaqkdaahaaWcbeqaaiabdsfaubaakiab=jfasjabcIcaOiqb=Lha5zaaiaGaeyOeI0Iaf4hVd0MbaKaadaWgaaWcbaGaeGimaadabeaakiabcMcaPaaa@3FEA@ , which is free of ρ.

These arguments suggest that for numerical evaluation, it is not necessary to consider all ρ values up to infinity. Instead, a moderately large enough value would suffice. Now the question comes down to how to decide on appropriate upper and lower bounds for ρ. The proof in Appendix A.1 requires max i ≠ j ∑ l = 1 p ( z i l − z j l ) 2 / ρ MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGagiyBa0MaeiyyaeMaeiiEaG3aaSbaaSqaaiabdMgaPjabgcMi5kabdQgaQbqabaGcdaaeWaqaaiabcIcaOiabdQha6naaBaaaleaacqWGPbqAcqWGSbaBaeqaaOGaeyOeI0IaemOEaO3aaSbaaSqaaiabdQgaQjabdYgaSbqabaGccqGGPaqkdaahaaWcbeqaaiabikdaYaaaaeaacqWGSbaBcqGH9aqpcqaIXaqmaeaacqWGWbaCa0GaeyyeIuoakiabc+caViabeg8aYbaa@4AB5@ be close to 0. Let _C_1 be some large positive number such that 1/_C_1 ≈ 0. If we take the upper bound of ρ to be _C_1 max i ≠ j ∑ l = 1 p ( z i l − z j l ) 2 MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGagiyBa0MaeiyyaeMaeiiEaG3aaSbaaSqaaiabdMgaPjabgcMi5kabdQgaQbqabaGcdaaeWaqaaiabcIcaOiabdQha6naaBaaaleaacqWGPbqAcqWGSbaBaeqaaOGaeyOeI0IaemOEaO3aaSbaaSqaaiabdQgaQjabdYgaSbqabaGccqGGPaqkdaahaaWcbeqaaiabikdaYaaaaeaacqWGSbaBcqGH9aqpcqaIXaqmaeaacqWGWbaCa0GaeyyeIuoaaaa@4805@ , then max i ≠ j ∑ l = 1 p ( z i l − z j l ) 2 / ρ MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGagiyBa0MaeiyyaeMaeiiEaG3aaSbaaSqaaiabdMgaPjabgcMi5kabdQgaQbqabaGcdaaeWaqaaiabcIcaOiabdQha6naaBaaaleaacqWGPbqAcqWGSbaBaeqaaOGaeyOeI0IaemOEaO3aaSbaaSqaaiabdQgaQjabdYgaSbqabaGccqGGPaqkdaahaaWcbeqaaiabikdaYaaaaeaacqWGSbaBcqGH9aqpcqaIXaqmaeaacqWGWbaCa0GaeyyeIuoakiabc+caViabeg8aYbaa@4AB5@ would be close to 0. In practice we suggest taking _C_1 = 100, which would give good approximation. Using a similar idea, we can find a lower bound for ρ. It is clear that when min i ≠ j ∑ l = 1 p ( z i l − z j l ) 2 MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGagiyBa0MaeiyAaKMaeiOBa42aaSbaaSqaaiabdMgaPjabgcMi5kabdQgaQbqabaGcdaaeWaqaaiabcIcaOiabdQha6naaBaaaleaacqWGPbqAcqWGSbaBaeqaaOGaeyOeI0IaemOEaO3aaSbaaSqaaiabdQgaQjabdYgaSbqabaGccqGGPaqkdaahaaWcbeqaaiabikdaYaaaaeaacqWGSbaBcqGH9aqpcqaIXaqmaeaacqWGWbaCa0GaeyyeIuoaaaa@4801@ /ρ → ∲ any non-diagonal element of K(ρ) will be 0 and the kernel matrix reduces to an identity matrix. Hence, if we pick a small enough number _C_2 such that 1/_C_2 → ∲, we can effectively set the lower bound of ρ to be _C_2 min i ≠ j ∑ l = 1 p ( z i l − z j l ) 2 MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGagiyBa0MaeiyAaKMaeiOBa42aaSbaaSqaaiabdMgaPjabgcMi5kabdQgaQbqabaGcdaaeWaqaaiabcIcaOiabdQha6naaBaaaleaacqWGPbqAcqWGSbaBaeqaaOGaeyOeI0IaemOEaO3aaSbaaSqaaiabdQgaQjabdYgaSbqabaGccqGGPaqkdaahaaWcbeqaaiabikdaYaaaaeaacqWGSbaBcqGH9aqpcqaIXaqmaeaacqWGWbaCa0GaeyyeIuoaaaa@4801@ . In practice we suggest take _C_2 = 0.1, which yields a good approximation.

A.3 derivation of normal equation (5)

Taking partial derivative of (6) with respect to β and writing in matrix notation, we have XT(y-μ). Similarly for α, we have K(**y**– μ) – λ Kα. The gradient vector is thus

q = [ X T ( y − μ ) K ( y − μ ) − λ K α ] . MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaacbmGae8xCaeNaeyypa0ZaamWaaeaafaqabeGabaaabaGae8hwaG1aaWbaaSqabeaacqWGubavaaGccqGGOaakcqWF5bqEcqGHsisliiWacqGF8oqBcqGGPaqkaeaacqWFlbWscqGGOaakcqWF5bqEcqGHsislcqGF8oqBcqGGPaqkcqGHsislcqaH7oaBcqWFlbWscqGFXoqyaaaacaGLBbGaayzxaaGaeiOla4caaa@4625@

(20)

Taking derivative of q with respect to β and α, we can get the following hessian matrix

H = − [ X T D X X T D K K D X λ K + K D K ] , MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaacbmGae8hsaGKaeyypa0JaeyOeI0YaamWaaeaafaqabeGacaaabaGae8hwaG1aaWbaaSqabeaacqWGubavaaGccqWFebarcqWFybawaeaacqWFybawdaahaaWcbeqaaiabdsfaubaakiab=reaejab=Tealbqaaiab=Tealjab=reaejab=HfaybqaaiabeU7aSjab=TealjabgUcaRiab=Tealjab=reaejab=TealbaaaiaawUfacaGLDbaacqGGSaalaaa@460E@

(21)

where D= Diag{μ i (1 – μ i )}. The Newton-Raphson iteration states that the parameter value at the (k + 1)th iteration can be updated by the following relationshipδ(k+1)= δ(k)– (H(k))-1 q(k),

where δ= (βT, αT)T. Substitute (20) and (21) into (22), we arrive at normal equation (7).

Availability

References

1. Subramanian A, Tamayo P, Mootha V, Mukherjee S, Ebert B, Gillette M, Paulovich A, Pomeroy S, Golub T, Lander E, Mesirov J: Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles. Proceedings of the National Academy of Sciences 2005, 102: 15545–15550. 10.1073/pnas.0506580102
   Article CAS Google Scholar
2. Eisenberg D, Graeber TG: Bioinformatic identification of potential autorine signaling loops in cancers from gene expression profiles. Nature Genetics 2001, 29: 295–300. 10.1038/ng718
   Article PubMed Google Scholar
3. Raponi M, Belly R, Karp J, Lancet J, Atkins D, Wang Y: Microarray analysis reveals genetic pathways modulated by tipifarnib in acute myeloid leukemia. BMC Cancer 2004, 4: 56. 10.1186/1471-2407-4-56
   Article PubMed Central PubMed Google Scholar
4. Dahlquist KD, Salomonis N, Vranizan K, Lawlor SC, Conklin BR: GenMAPP, a new tool for viewing and analyzing microarray data on biological pathways. Nature Genetics 2002, 31: 19–20. 10.1038/ng0502-19
   Article CAS PubMed Google Scholar
5. Grosu P, Twonsend JP, Hartl DL, Cavalieri D: Pathway Processor: A tool for integrating whole-genome expression results into metabolic networks. Genome Research 2002, 12: 1121–1126. 10.1101/gr.226602
   Article PubMed Central CAS PubMed Google Scholar
6. Doniger SW, Salomonis N, Dahlquist KD, Vranizan K, Lawlor SC, Conklin BR: MAPPFinder: using Gene Ontology and GenMAPP to create a global gene-expression profile from microarray data. Genome Biology 2003, 4: R7. 10.1186/gb-2003-4-1-r7
   Article PubMed Central PubMed Google Scholar
7. Vapnik V: Statistical Learning Theory. New York: Wiley; 1998.
   Google Scholar
8. Schölkopf B, Smola A: Learning with Kernels. Cambridge, Massachusetts: MIT press; 2002.
   Google Scholar
9. Liu D, Lin X, Ghosh D: Semiparametric regression of multi-dimensional genetic pathway data: least squares kernel machines and linear mixed models. Biometrics 2007, 63(4):1079–1088.
   Article PubMed Central PubMed Google Scholar
10. Breslow N, Clayton D: Approximate inference in generalized linear mixed models. Journal of the American Statistical Association 1993, 88: 9–25. 10.2307/2290687
    Google Scholar
11. Wei Z, Li H: Nonparametric pathway-based regression models for analysis of genomic data. Biostatistics 2007, 8(2):265–284. 10.1093/biostatistics/kxl007
    Article PubMed Google Scholar
12. Sprague R, Ed: Proceedings of the 39th Annual Hawaii International Conference on System Sciences. Los Alamitos: IEEE; 2006. [CD ROM version]
    Google Scholar
13. Goeman JJ, Geer SA, de Kort F, van Houwelingen HC: A global test for groups of genes: testing association with a clinical outcome. Bioinformatics 2004, 20: 93–99. 10.1093/bioinformatics/btg382
    Article CAS PubMed Google Scholar
14. Goeman JJ, Geer SA, van Houwelingen HC: Testing against a high dimensional alternative. Journal of the Royal Statistical Society: Series B 2006, 68: 477–493. 10.1111/j.1467-9868.2006.00551.x
    Article Google Scholar
15. Dhanasekaran S, Barrette T, Ghosh D, Shah R, Varambally S, Kurachi K, Pienta K, Rubin M, Chinnaiyan A: Delineation of prognostic biomarkers in prostate cancer. Nature 2001, 412(6849):822–6. 10.1038/35090585
    Article CAS PubMed Google Scholar
16. Kimeldorf G, Wahba G: Some results on Tchebycheffian spline functions. Journal of Mathematical Analysis and Applications 1970, 33: 82–95. 10.1016/0022-247X(71)90184-3
    Article Google Scholar
17. Self SG, Liang KY: Asymptotic properties of maximum likelihood estimators and likelihood ratio tests under non-standard conditions. Journal of the American Statistical Association 1987, 82: 605–610. 10.2307/2289471
    Article Google Scholar
18. Goeman JJ, Bühlmann P: Analyzing gene expression data in terms of gene sets: methodological issues. Bioinformatics 2007, 23: 980–987. 10.1093/bioinformatics/btm051
    Article CAS PubMed Google Scholar
19. Zhang D, Lin X: Hypothesis testing in semiparametric additive mixed models. Biostatistics 2002, 4: 57–74. 10.1093/biostatistics/4.1.57
    Article CAS Google Scholar
20. Davies R: Hypothesis testing when a nuisance parameter is present only under the alternative. Biometrika 1977, 64: 247–254. 10.2307/2335690
    Article Google Scholar
21. Davies R: Hypothesis testing when a nuisance parameter is present only under the alternative. Biometrika 1987, 74: 33–43.
    Google Scholar
22. le Cessie S, van Houwelingen J: Goodness of fit tests for generalized linear models based on random effect models. Biometrics 1995, 51(2):600–614. 10.2307/2532948
    Article CAS PubMed Google Scholar
23. McCullagh P, Nelder J: Generalized Linear Models. New York: Chapman & Hall; 1989.
    Book Google Scholar

Download references

Acknowledgements

Liu and Lin's research was supported by a grant from the National Cancer Institute (R37 CA-76404). Ghosh's research was supported by a grant from the National Institute of Health (R01 GM-72007). We thank A. Chinnaiyan for providing the prostate cancer data. We also thank the two anonymous referees for their comments which helped improve the manuscript.

Author information

Authors and Affiliations

1. Center for Statistical Sciences, Brown University, Providence, RI, 02912, USA
   Dawei Liu
2. Departments of Statistics and Public Health Sciences, Pennsylvania State University, University Park, PA, 16802, USA
   Debashis Ghosh
3. Department of Biostatistics, Harvard School of Public Health, 655 Huntington Avenue, Boston, MA, 02115, USA
   Xihong Lin

Authors

1. Dawei Liu
   You can also search for this author inPubMed Google Scholar
2. Debashis Ghosh
   You can also search for this author inPubMed Google Scholar
3. Xihong Lin
   You can also search for this author inPubMed Google Scholar

Corresponding author

Correspondence toDawei Liu.

Additional information

Authors' contributions

DL performed statistical analysis. All authors participated in development of the methods and preparation of the manuscript. All authors read and approved the final manuscript.

Rights and permissions

This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Reprints and permissions

About this article

Cite this article

Liu, D., Ghosh, D. & Lin, X. Estimation and testing for the effect of a genetic pathway on a disease outcome using logistic kernel machine regression via logistic mixed models.BMC Bioinformatics 9, 292 (2008). https://doi.org/10.1186/1471-2105-9-292

Download citation

• Received: 19 February 2008
• Accepted: 24 June 2008
• Published: 24 June 2008
• DOI: https://doi.org/10.1186/1471-2105-9-292

Keywords