Smooth-threshold multivariate genetic prediction withunbiased model selection
Masao Ueki∗ and Gen Tamiya†
Abstract — We develop a new genetic prediction method, smooth-threshold multivariate genetic prediction, using single
nucleotide polymorphisms (SNPs) data in genome-wide association studies (GWASs). Our method consists of two stages.
At the first stage, unlike the usual discontinuous SNP screening as used in the gene score method, our method continuously
screens SNPs based on the output from standard univariate analysis for marginal association of each SNP. At the second
stage, the predictive model is built by a generalized ridge regression simultaneously using the screened SNPs with SNP
weight determined by the strength of marginal association. Continuous SNP screening by the smooth-thresholding not only
makes prediction stable but also leads to a closed form expression of generalized degrees of freedom (GDF). The GDF leads
to the Stein’s unbiased risk estimation (SURE) which enables data-dependent choice of optimal SNP screening cutoﬀ without
using cross-validation. Our method is very rapid because computationally expensive genome-wide scan is required only
once in contrast to the penalized regression methods including lasso and elastic net. Simulation studies which mimic real
GWAS data with quantitative and binary traits demonstrate that the proposed method outperforms the gene score method
and genomic best linear unbiased prediction (GBLUP), and also shows comparable or sometimes improved performancewith the lasso and elastic net being known to have good predictive ability but with heavy computational cost. Application
to whole-genome sequencing (WGS) data from the Alzheimer’s Disease Neuroimaging Initiative (ADNI) exhibits that the
proposed method shows higher predictive power than the gene score and GBLUP methods.

Keyword: Genetic prediction; marginal association screening; model selection; smooth-thresholding.
1

Introduction

Genome-wide association study (GWAS) is a popular tool for discovering disease-susceptibility genes
using large number of single nucleotide polymorphisms (SNPs) without prior knowledge. Apart fromdiscovery of susceptibility genes, prediction of individual’s phenotype from high-dimensional genetic
information, termed as a genetic prediction, is an important task for personalized medicine. Currently,
researchers are exploring the most eﬀective way of building genetic prediction models (Purcell et al.
2009). In this paper, we develop a new statistical approach, smooth-threshold multivariate genetic
prediction, for building genetic predictive models with input of large-scale genome-wide SNPs data.
We consider standard multiple regression model but with high-dimensional predictor variables. To
be specific, y = (y1 , . . . , yn )T represent response variables of individual’s phenotype data modeled by aconditional distribution given predictor variables X = (X1 , . . . , Xp ) observed for n individuals, in which
X j = (x1,j , . . . , xn,j )T for j ∈ M = {1, . . . , p}. The conditional expectation of yi given xi = (xi,1 , . . . , xi,p )is assumed to be a linear combination η{E(yi |xi )} = xi β, where η is some known monotone function
and β is a vector of regression coeﬃcients. In this paper, we consider linear regression with identity
map η for quantitative trait such as clinical characteristics, and logistic regression with logit function
η for binary trait such as aﬀected/unaﬀected status. Each X j is either genotype at a SNP site or othercovariate such as sex, age, body mass index (BMI), smoking status, alcohol consumption and principal
components for population stratification. Each SNP can take one of three possible genotypes, gg, gG
and GG, where g and G denote minor and major alleles at the SNP site, respectively. If X j represents theobserved count of minor allele g at a SNP site, X j takes a value from {0, 1, 2}. Under the Hardy–Weinberg
equilibrium (HWE), the observed count of minor allele g at each SNP follows a binomial distribution
with parameter f ∈ [0, 0.5] called a minor allele frequency (MAF), i.e. frequency of the minor allele
g in general population. Quality controls (QCs) are often conducted to remove low-quality SNPs bychecking HWE and missing rates as well as low MAF SNPs. Even after those QCs, large number of
SNPs still remain. Since sample sizes are usually far less than the number of SNPs, the predictive∗ Biostatistics Center, Kurume University, 67 Asahi-Machi, Kurume, Fukuoka 830-0011, Japan, E-mail: uekimrsd@nifty.com† Tohoku Medical Megabank Organization, Tohoku University, 2-1 Seiryo-Machi, Aoba-Ku, Sendai 980-8573, Japanmodeling in GWAS faces the p ≫ n problem. The p ≫ n condition hampers multiple regression thatfits simultaneously using p predictors X.
Standard GWAS analysis conducts marginal association scan between y and each X j independently,
i.e. a univariate analysis which tests the slope parameter in univariate regression model, followed bymultiple test using a Bonferroni correction with a stringent significance level (e.g. p-value less than
5 × 10−8 ) in order to control the rate of false positive findings. Meanwhile, suppose that X does notinclude covariates and consists of SNPs only. Let T j (y, X) represent a non-negative test statistic for
testing association between jth SNP X j and y as a function of y and X, and the corresponding inclusion
threshold be t > 0. For example, t is a chi-squared quantile at a given p-value cutoﬀ for chi-squared
test statistics T j (y, X). The resulting SNP set from a marginal association screening at a threshold t isdefined by A = {j ∈ M : T j (y, X) > t}. Purcell et al. (2009) proposed a gene score method which simply
averages each genotype data weighted by estimated eﬀect size for each SNP in A. Warren et al. (2013)
consider multiple regression for SNPs in A, called a multivariate gene score method.
In the purpose of prediction, the cutoﬀ t can be chosen in terms of prediction ability. However,
evaluating prediction ability is not straightforward unlike in traditional setting without screening. It
is known that, the screening invalidates traditional statistical procedures, called an winner’s curse
eﬀect. Analogous problem arises in the context of prediction modeling. Actually, simulation studies
as well as examination on real GWAS datasets reported that screening leads to overfitting. In Ueki
and Tamiya (2016), we show that the screening can deflate the residual sum of squares (RSS) comparedwith the RSS without screening, so that the RSS becomes too optimistic. Since screening complicatesthe behavior of RSS, naive use of RSS is unwarranted in measuring prediction ability. Instead, we canuse cross-validation (or sample splitting) which divides the training data into two parts, one of which
is used for ranking SNPs and remaining is used to construct a predictive model. Purcell et al. (2009)choose an optimal inclusion cutoﬀ by cross-validation.
Although cross-validation takes into account of the screening, reduced sample sizes in training
stage may lose predictive power, which is a severe concern when sample sizes are small. Five or
ten-folds cross-validation is commonly used in model selection. For example, the SparSNP program(Abraham et al., 2012) implementing penalized regression methods, the lasso and elastic net, searches
for entire genome-wide SNPs data without SNP screening. SparSNP selects the tuning parameter by
k-fold cross-validation with default setting of k = 10. Repeated genome-wide scans needed at eachcandidate tuning parameter and multiple runs of model fitting in each fold increase computationalcost. For large-scale data such as the whole-genome sequencing (WGS) data, heavy computational
cost critically limits the applicability although penalized methods are known to give better predictive
power than the simpler gene score method.
In this paper, we develop a new predictive modeling approach, a smooth-threshold multivariate
genetic prediction, which is really applicable to large-scale genome-wide data such as WGS data while
preserving high prediction ability. Our method consists of two stages. At the first stage, our method
continuously screens SNPs based on the output from standard univariate analysis for marginal association of each SNP. At the second stage, the predictive model is built by a generalized ridge regression
simultaneously using the screened SNPs with SNP weight determined by the strength of marginalassociation reflecting the uncertainty of inclusion. Since the final predictive model is essentially built
in multiple regression model as in the sure independence screening, the correlations between predictor
variables are accounted for (See also Warren et al. (2013)). Marginal association signal is used only forpenalizing each regression coeﬃcient. Our method is very rapid because computationally expensive
genome-wide scan is required only once in contrast to the penalized methods which need genome-wide
scan several times. Our proposal can be seen as a smoothed version of multiple regression after single
SNP-GWAS screening of predictor variables at some p-value cutoﬀ, in which the discontinuous process in screening is replaced by a continuous function. The resulting continuity makes the prediction
stable in the sense of Breiman (1996). The continuity in SNP screening also leads to a closed form
expression of generalized degrees of freedom (GDF; Ye, 1998), and allows an application of Stein’s
unbiased risk estimation (SURE). While the Mallows’ Cp (Mallows, 1973) with the usual degrees offreedom is no longer unbiased model selection criterion due to the eﬀect of screening, we can readily
construct an unbiased Cp -type model selection criterion using the GDF. It allows data-dependent choiceof optimal SNP inclusion cutoﬀ without relying on cross-validation. The eﬀect of screening is properly
accounted for by the SURE’s unbiasedness. Since no cross-validation is needed, computationally expensive genome-wide scan is required only once in ranking SNPs. We also extend to generalized linear
models and propose a loglikelihood-based Cp -type model selection criterion. Simulation studies whichmimic real SNP-GWAS data for both quantitative and binary traits show that the proposed method
gives better performance than gene score and genomic best linear unbiased prediction (GBLUP) andattains a comparable or sometime improved prediction performance with the lasso and elastic net
in SparSNP program. Application to large-scale WGS data from Alzheimer’s Disease Neuroimaging
Initiative (ADNI) exhibits that the proposed method gives higher predictive performance than both thegene score and GBLUP methods.
2

Materials and Methods

Here we consider linear multiple regression model, y = µ + ϵ, where µ = E(y|X) = Xβ, ϵ ∼ N(0, σ2 In ), X
is a p-dimensional design matrix and β is the corresponding p regression coeﬃcients. Since p is muchlarger than n in typical GWAS data, some dimensionality reduction is required. Sparsity assuming thatmany components of β are zero would be a realistic assumption. If susceptible SNPs show relatively
large marginal signal, marginal association screening eﬀectively reduces the dimensionality. The genescore method (Purcell et al., 2009) and its multivariate generalization (Warren et al., 2013) use upperranked SNPs in marginal association, A = {j ∈ M : T j (y, X) > t}, for a given cutoﬀ value t > 0. Althoughdimensionality is eﬀectively reduced, discontinuity in y present in the screening process in A may incur
instability of prediction, i.e. small change in data can make large changes in the prediction (Breiman,
1996). To address the discontinuity issue, we use a smooth-thresholding proposed by Ueki (2009). To
be specific, we propose to estimate the regression coeﬃcients by
(
) (
)
β̌A
ǦA (I|A| − ĎA )XAT y
β̌ =
=
,
(2.1)
0
β̌Ac
where Ac indicates the complement set of A, ǦA = {(I|A| − ĎA )(ΣAA + λI|A| ) + τĎA }−1 , Σ = XT X, ΣAA =(Σ jk ) j∈A,k∈A , γ and τ are non-negative tuning parameters and λ > 0 is a small constant to avoid singularityof ǦA . The corresponding prediction of yi is then µ̌i (y) = XiT β̌. Here Ď j is an adaptive lasso smooththresholding function defined by
1+γĎ j = min[1, {t/T j (y, X)} 2 ].

(2.2)

Since Ď j = 1 if and only if T j (y, X) ≤ t, the screened set A with Ď j is the same as that with D̂ j = 1{T j (y,X)≤t} ,where 1{·} denotes the indicator function. It can be seen that Ď j replaces the discontinuous screening
process D̂ j by a continuous function. As a result, µ̌i (y) turns out to be continuous in y.
The regression coeﬃcient for the screened set in (2.1), β̌A , can be seen as a solution to
XAT (XA β̌A − y) + WA β̌A = 0,

(2.3)

with WA = diag(W j : j ∈ A) where W j = λ + τĎ j /(1 − Ď j ), which is the minimizer of a generalized∑
ridge regression loss, ||y − XA βA ||2 + j∈A β2j W j , with respect to βA . Ridge weight for each predictor
variable, W j , represents uncertainty of marginal association screening. If the marginal association isvery weak, we have Ď j ≈ 1 and large W j , then the corresponding regression coeﬃcient is strongly
shrunken towards zero. If the marginal association is strong, we have Ď j ≈ 0 and W j ≈ λ, then thecorresponding regression coeﬃcient is less penalized. From the fact that the winner’s curse eﬀect
produces larger selection bias for small regression coeﬃcient, it is expected that the above penalization
decreases the selection bias.

Predictive power largely depends on the choice of t. It may be done using cross-validation by
dividing a dataset into test and training samples (Warren et al., 2013). Cross-validation takes into
account sampling variability due to the screening. However, repeated genome-wide scans to obtainthe screened set A needed in cross-validation incurs computational burden. It is also concerned thatthe reduction in training sample sizes decreases the predictive power of the model. Instead of crossvalidation, we propose a Cp -type criterion based on SURE using GDF. The continuity of µ̌i (y) in yleads to a closed-form expression of GDF. In what follows, we consider p-value cutoﬀ α instead of t
by a one-to-one transformation t = F−1 (1 − α), where F−1 is a quantile function of the distribution of
T j (y, X) under the null hypothesis of no marginal association such as F or χ2 distribution. An optimal
α is determined by minimizing the Cp -type criterion within a range of α for search, [αmin , αmax ]. The
proposed procedure is outlined at the end of this section. It is noteworthy that the computational
intensive genome-wide scan is required only once in the single-SNP association screening at Step 1.
Subsequent Steps 2–4 are performed on the reduced set of SNPs whose single-SNP association p-value
is less than αmax . More details of this section including formulas, derivations, extension to generalized
linear models and additional descriptions are given in Ueki and Tamiya (2016).
Outline of algorithm
Step 1. Perform single-SNP association analysis for p SNPs.
Step 2. Screen SNPs whose single-SNP association p-value is less than αmax .
Step 3. Fix γ and τ as suggested in main text, and select an optimal α from candidate values in [αmin , αmax ]∑by minimizing the Cp -type criterion: C(α) = ni=1 {yi − µ̌i (α)}2 + 2σ2 GDF(α). Explicit formulas are
given in Ueki and Tamiya (2016).
Step 4. Compute β̌ by (2.1) using the selected α.
3

Acknowledgment

This work was carried out under the ISM General Cooperative Research 1 (2015-ISM-CRP-1013) andwas partially supported by a Grant-in-Aid for Young Scientist (B) (25870074) and Grants-in-Aid for
Scientific Research (C) (25330049 and 25460403).
References[1] Abraham G., Kowalczyk A., Zobel J. and Inouye, M. (2012). SparSNP: fast and memory-eﬃcient
analysis of all SNPs for phenotype prediction. BMC Bioinfomatics, 13:88.
[2] Breiman L. (1996). Heuristics of instability and stabilization in model selection. Annals of Statistics,
24:2350–2383.
[3] Mallows C. L. (1973). Some comments on Cp . Technometrics, 15:661–675.
[4] Purcell SM., Wray N. R., Stone J. L., Visscher P. M., O’Donovan M. C., Sullivan P. F., Sklar P. andConsortium International Schizophrenia. (2009). Common polygenic variation contributes to riskof schizophrenia and bipolar disorder. Nature, 460:748–52.
[5] Ueki M. (2009). A note on automatic variable selection using smooth-threshold estimating equation.
Biometrika, 96:1005–1011.
[6] Ueki M. and Tamiya G. (2016). Smooth-threshold multivariate genetic prediction with unbiased
model selection. Genetic Epidemiology, 40:233–243.
[7] Warren H., Casas J. P., Hingorani A., Dudbridge F. and Whittaker J. (2013). Genetic prediction ofquantitative lipid traits: comparing shrinkage models to gene scores. Genetic Epidemiology, 38:72–83.
[8] Ye J. (1998). On measuring and correcting the eﬀects of data mining and model selection. Journal of
the American Statistical Association, 93, 120–131.