# A Multinomial Ordinal Probit Model with Singular Value Decomposition Method for a Multinomial Trait

Authors: Soonil Kwon; Mark O. Goodarzi; Kent D. Taylor; Jinrui Cui; Y.-D. Ida Chen; Jerome I. Rotter; Willa Hsueh; Xiuqing Guo
Journal: Journal of Probability and Statistics
 (2012)

---

## Abstract

We developed a multinomial ordinal probit model with singular value decomposition for testing a large number of single nucleotide polymorphisms (SNPs) simultaneously for association with multidisease status when sample size is much smaller than the number of SNPs. The validity and performance of the method was evaluated via simulation. We applied the method to our real study sample recruited through the Mexican-American Coronary Artery Disease study. We found 3 genes (SORCS1, AMPD1, and PPARα) to be associated with the development of both IGT and IFG, while 5 genes (AMPD2, PRKAA2, C5, TCF7L2, and ITR) with the IGT mechanism only and 6 genes (CAPN10, IL4, NOS3, CD14, GCG, and SORT1) with the IFG mechanism only. These data suggest that IGT and IFG may indicate different physiological mechanism to prediabetes, via different genetic determinants.


---

## Body

## 1. Introduction

Genome-wide association studies (GWASs) examine genetic variants across the entire genome to improve the understanding of genetic components underlying complex human disease. With whole-genome genotyping techniques that allow GWAS to involve hundreds of thousands of single nucleotide polymorphisms (SNPs), many studies have successfully identified novel genetic components for many diseases or related quantitative traits. However, the sample size is often limited due to the difficulty of recruiting patients and/or the cost of research. This leads to the situation that the number of SNPs (m) is much larger than the samples (n) available, that is, m≫n, and makes the traditional statistical methods unsuitable for analyzing multiple SNPs simultaneously. Most of current GWASs deal with this shortcoming by performing single SNP association test, which analyzes one SNP at a time and results in a huge multiple testing problem. These motivated us to develop methods that can avoid the multiple testing problem, in other words, methods that can evaluate multiple SNPs simultaneously when m≫n. We first introduced the iterative Bayesian variable selection (IBVS) method [1], which analyzes all SNPs simultaneously when m≫n and uses the Bayesian variable selection [2] iteratively to find SNPs that are associated with disease. The method was successfully applied to the simulated rheumatoid arthritis data provided by the Genetic Analysis Workshop 15 (GAW15). We later introduced the Bayesian classification with singular value decomposition (BCSVD) method [3]. The method applies the singular value decomposition (SVD) to the covariate matrix, which is usually the genotype data in GWAS and reduces the dimension of parameters to be estimated to the number of samples. This makes the method feasible to handle multiple SNPs simultaneously when m≫n. The validation of the method was demonstrated by applying to the simulated data provided by GAW16. We now extend our method from binary disease to the analysis of polytomous ordinal response variables. We propose here a multinomial ordinal probit model with singular value decomposition method. We show the validity of the newly developed method by applying it to simulated data sets as well as to a real study sample to identify genes contributing to two different mechanisms for prediabetes, namely, impaired glucose tolerance (IGT) and impaired fasting glucose (IFG). With the simulated data sets, we demonstrate that this new method is superior to single SNP analysis method and, with the real data, identify different genes for each mechanism.


## 2. Method and Materials

### 2.1. Multinomial Probit Model with Singular Value Decomposition

Logit and probit models are statistical models that are widely used for the analysis of categorical (ordinal/nominal) data. The difference between these two models is the choice of the link function relating the linear predictor to the expected value; the probit model uses the inverse normal cumulative distribution, and the logit model uses the logit transformation. As discussed by Greene [4], in most cases, the choice of the link function is largely a matter of taste. We utilized the probit model here to analyze data with polytomous ordinal response variables. In general, the multinomial ordinal probit model can be expressed by latent (unobserved) continuous variables associated with categorical responses. Let us assume that responses y1,y2,…,yn are observed, where yi takes one of the J-ordered categories and θ1,…,θJ are real numbers of bin boundaries, which satisfy that -∞=θ1≤⋯≤θJ=∞. As discussed by Albert and Chib [5], we denote that z1,z2,…,zn are latent continuous random variables. We assume that the latent variable (zi) associated with a categorical outcome (yi) can be explained in terms of an underlying linear model, and that the observed response yi has category j if and only if zi falls between θj-1 and θj. The multinomial ordinal probit model is equivalent to the following model:
(2.1)zi=xiβ+ϵi,ϵi~N(0,1),i=1,…,n,yi=j⟺θj-1<zi≤θj,j=1,…,J,
where xi is a 1×m vector of the explanatory variables for the ith sample, and β is a m×1 vector of parameters to be estimated. In vector-matrix notation, we can have the multinomial ordinal probit model
(2.2)z=Xβ+ϵ,
where z is the n×1 vector of latent variables, X is the n×m matrix of the explanatory variables, β is the m×1 vector of unknown regression coefficients, and ϵ follows an independent standard multivariate normal distribution,ϵ~N(0,In). By applying SVD to the matrixX in (2.2), when rank(X)=n, the matrix can be expressed as X′=ADF′, where A is the m×n singular value factor loading matrix with orthonormal columns so that A′A=In,F is the n×n SVD orthogonal factor matrix with F′F=FF′=In, andD=diag(d1,…,dn), the diagonal matrix of positive singular values, ordered as d1≥⋯≥dn>0. When rank(X)=r<n, the smallest n-r singular values in D are replaced with 0, that is, d1≥⋯≥dr>dr+1=⋯=dn=0. Therefore, in the product X′=ADF′, the last n-r columns of both A and F for which dr+1=⋯=dn=0 are ignored since they interact with the block of zeros in D. Hence, this leads to another form of SVD, X′=ArDrFr′, that is, the product of the first r columns of A, the upper r×r block of D, and the first r columns of F. Since the difference between the both scenario is only in dimension of matrices in SVD, we assume that rank(X)=n in the rest part of the paper for convenience. Thus, the model in (2.2) with the SVD of X can be written as follows:
(2.3)z=Xβ+ϵ=(ADF′)′β+ϵ=FDA′β+ϵ=Lγ+ϵ,
where L=FD and γn×1=A′n×mβm×1. Therefore, z, the n×1 vector of latent variables in (2.3), has a multivariate normal distribution, that is, z~N(Lγ,In). As shown in (2.3), γ is expressed by a linear combination of the original parameters (β). Hence, we call γ as the vector of superfactors. The model in (2.3) represents a massive dimension reduction from m to n parameters. The regression model with m parameters reduced to that with n parameters derived from the SVD of the covariate matrix X. Therefore, the statistical inference on the original parameter turns into the superfactors. Let pi=(pi1,…,piJ) denote the vector of probabilities associated with the assignment of the ith sample into categories 1,…,J, where pij denote the probability that a sample falls into category j. From (2.1) and (2.3), it follows that
(2.4)pij=∫θj-1θjϕ(z-liγ)dz=Pr(θj-1<Zi<θj)=Φ(θj-liγ)-Φ(θj-1-liγ),
where ϕ(·) and Φ(·) denote the probability density function and the cumulative density function of the standard normal distribution, respectively, and li is the ith row of matrix L=FD. Let y=(y1,…,yn) denote the vector of responses observed for all samples. Then, the probability of observing data y is given as follow:
(2.5)Pr(y∣pi)=∏i=1n∏j=1JpijI{yi=j}.
From (2.4) and (2.5), the log likelihood function for (γ,θ) can be written as
(2.6)ln⁡L(γ,θ)=∑i=1n∑j=1JI{yi=j}ln⁡[Φ(θj-liγ)-Φ(θj-1-liγ)].


### 2.2. Model Fitting with Maximum Likelihood Estimation

The maximum likelihood estimates (MLEs) of the superfactors,γ, in (2.3) can be obtained by the iteratively reweighted least squares (IRLSs) procedure [6] using the log likelihood function for (γ,θ) in (2.6). The procedure can be briefly described as follows. Let η denote the vector of all model parameters, that is, η=(θ2,…,θJ-1,γ1,…,γn-J+2). Note that θ1 and θJ are not included in this vector because their values are assumed to be 0 and ∞, respectively, for the purpose of model identifiability. Also note that the (J-2) smallest singular values together with their corresponding factors are dropped from the parameters since the number of parameters must not exceed the number of samples. Assuming thatJ=4, define that
(2.7)Ci=[100-1100-1100-1],Li=[100-li010-li001-li],
and ℋi=diag(fi1,…,fiJ-1), where fij denotes the derivative of the standard normal cumulative distribution function at θj-liγ. Take Wi=diag(pi), where pi is the J×1 vector of probabilities that the ith individual falls in each category, that is, pi=(pi1,…,piJ)′, and let 𝒩i be a J×1 vector of observation, that is, 𝒩i=(I{yi=1},…,I{yi=J})′. After initialization of all elements, the iteration s+1(s=1,2,…) can be written as
(2.8)ηs+1=ηs+(∑i=1nLi′Hi(s)Ci′Wi-1(s)CiHi(s)Li)-1∑i=1nLi′Hi(s)Ci′Wi-1(s)(Ni-pi(s)).
The MLE of η can be found by performing the process recursively until the change between ηs+1 and ηs is negligible.


### 2.3. General Solution for the Original Parameters

We have discussed how to estimate the superfactor (γ) in (2.3) thus far. Since the primary interest is to find SNPs that are significantly associated with a disease, it is necessary to transform the superfactor (γ) to the original parameters (β) in (2.1). The equation γ=A′β in (2.3) can be utilized for the transformation even though A is m×n nonsquare matrix. As discussed by Graybill [7], the unique solution for β can be achieved by taking the generalized inverse matrix of A′ as A since A′A=In. Therefore, the unique solution for SNP effect (β) can be calculated by β=Aγ.


### 2.4. Selection of Significant SNPs

Finding significant SNPs is the same as testing if each SNP effect(βi,i=1,…,m) is statistically significant, that is, testing the hypothesis: H0: βi=0 versus H1: βi≠0, i=1,…,m. The simple method is to use Wald’s test statistic, which forms (β̂-β)/se(β̂) and assumes a normal distribution. However, when m≫n, it is hard to calculate se(β̂) directly from the data. We therefore utilized permutation test to select significant SNPs. The rationale behind the test is that, under the null hypothesis, the estimate of β obtained from the raw (unpermuted) data is similar to the estimate of β obtained from the permuted data. That is, the difference between two estimates is closed to zero under H0. With this idea, we can construct Wald’s test statistic as follows. Let β̂i (i=1,…,m) be the estimate of the ith SNP effect from the raw data and β̂ik (k=1,…,K) be the estimate of the ith SNP effect from the kth-permuted data. Let us define βidk as the difference between β̂i and β̂ik, that is, βidk=β̂i-β̂ik. Then, Wald’s test statistic can be as follows:
(2.9)Λi=β¯idse(β¯id),i=1,…,m,
where β¯id is the sample mean of βidk’s, which is β¯id=(1/K)∑k=1Kβidk, and se(β¯id) is the standard error of β¯id. Under the null hypothesis, the statistic Λi defined in (2.9) follows approximately standard normal distribution when k is large. P value for rejecting the null hypothesis at a significance level α=0.05 can be utilized to identify significant SNPs.


### 2.5. Application of the Multinomial Probit Model with SVD

#### 2.5.1. Simulated Multinomial Ordinal Data

The validity of the proposed method was evaluated using simulated data sets. The procedure of data generation was composed of three steps: generating genotype data with certain genetic model, generating the latent variable, and defining the disease status variable by applying the predefined bin boundaries. The brief scheme of each step is as follows: we first generated 10 sets of the simulated genotype data under an additive genetic model, each set consists of 100 samples and 1000 SNPs. From (2.1), we can notice that the latent variable (zi) consists of two parts: the expected value (xiβ) and the random error (ϵi). In order to generate the expected value, we assumed that, for each sample, 9 out of the 1000 SNPs (every 101th SNP, except the last one) contribute to disease status xiβ=β1·SNP101+β2·SNP201+⋯+β1·SNP801+β2·SNP901, where β1 and β2 are set as −1 and 1, respectively. Hence, the latent variable can be obtained from the sum of the expected values (xiβ) and the random error generated from standard normal distribution. We then generated disease status variable (yi) assuming 3 disease development stages. Therefore, when applying the proposed method to the simulated data sets, we would expect 9 strong signals corresponding to each of the 9 disease-associated SNPs. We also compared results obtained from the proposed method with that from single SNP analysis with multinomial ordinal probit model.


### 2.6. Mexican-American Coronary Artery Disease (MACAD) Study

We also applied the proposed method to study sample recruited through the Mexican-American Coronary Artery Disease (MACAD) study [8, 9]. The study population consists of probands who are Mexican American aged between 45 and 75 with coronary artery disease: spouses of probands, adult offspring (≥18), and their spouses. For the offspring generation, we performed oral glucose tolerance test and genotyped 132 SNPs in 32 genes selected based on a prior relationship to insulin physiology. The goal of the study herein was to identify genes involved in the development of IGT and/or IFG, where IGT was defined as a 2 hr glucose level between 140 and 199 mg/dL and IFG defined as a fasting glucose level between 100 and 125 mg/dL. In order to identify and compare genes affecting the development of IGT and/or IFG, we generated two study samples, for which each sample has 3 disease stages (D1) both 2 hr and fasting glucoses normal (N/N)(n1=60), IGT only (IGT/N) (n2=31) and IGT and IFG (IGT/IFG) (n3=15) (D2) both 2 hr and fasting glucoses normal (N/N)(n1=60), IFG only (N/IFG) (n2′=34) and IGT and IFG (IGT/IFG) (n3=15).


## 2.1. Multinomial Probit Model with Singular Value Decomposition

Logit and probit models are statistical models that are widely used for the analysis of categorical (ordinal/nominal) data. The difference between these two models is the choice of the link function relating the linear predictor to the expected value; the probit model uses the inverse normal cumulative distribution, and the logit model uses the logit transformation. As discussed by Greene [4], in most cases, the choice of the link function is largely a matter of taste. We utilized the probit model here to analyze data with polytomous ordinal response variables. In general, the multinomial ordinal probit model can be expressed by latent (unobserved) continuous variables associated with categorical responses. Let us assume that responses y1,y2,…,yn are observed, where yi takes one of the J-ordered categories and θ1,…,θJ are real numbers of bin boundaries, which satisfy that -∞=θ1≤⋯≤θJ=∞. As discussed by Albert and Chib [5], we denote that z1,z2,…,zn are latent continuous random variables. We assume that the latent variable (zi) associated with a categorical outcome (yi) can be explained in terms of an underlying linear model, and that the observed response yi has category j if and only if zi falls between θj-1 and θj. The multinomial ordinal probit model is equivalent to the following model:
(2.1)zi=xiβ+ϵi,ϵi~N(0,1),i=1,…,n,yi=j⟺θj-1<zi≤θj,j=1,…,J,
where xi is a 1×m vector of the explanatory variables for the ith sample, and β is a m×1 vector of parameters to be estimated. In vector-matrix notation, we can have the multinomial ordinal probit model
(2.2)z=Xβ+ϵ,
where z is the n×1 vector of latent variables, X is the n×m matrix of the explanatory variables, β is the m×1 vector of unknown regression coefficients, and ϵ follows an independent standard multivariate normal distribution,ϵ~N(0,In). By applying SVD to the matrixX in (2.2), when rank(X)=n, the matrix can be expressed as X′=ADF′, where A is the m×n singular value factor loading matrix with orthonormal columns so that A′A=In,F is the n×n SVD orthogonal factor matrix with F′F=FF′=In, andD=diag(d1,…,dn), the diagonal matrix of positive singular values, ordered as d1≥⋯≥dn>0. When rank(X)=r<n, the smallest n-r singular values in D are replaced with 0, that is, d1≥⋯≥dr>dr+1=⋯=dn=0. Therefore, in the product X′=ADF′, the last n-r columns of both A and F for which dr+1=⋯=dn=0 are ignored since they interact with the block of zeros in D. Hence, this leads to another form of SVD, X′=ArDrFr′, that is, the product of the first r columns of A, the upper r×r block of D, and the first r columns of F. Since the difference between the both scenario is only in dimension of matrices in SVD, we assume that rank(X)=n in the rest part of the paper for convenience. Thus, the model in (2.2) with the SVD of X can be written as follows:
(2.3)z=Xβ+ϵ=(ADF′)′β+ϵ=FDA′β+ϵ=Lγ+ϵ,
where L=FD and γn×1=A′n×mβm×1. Therefore, z, the n×1 vector of latent variables in (2.3), has a multivariate normal distribution, that is, z~N(Lγ,In). As shown in (2.3), γ is expressed by a linear combination of the original parameters (β). Hence, we call γ as the vector of superfactors. The model in (2.3) represents a massive dimension reduction from m to n parameters. The regression model with m parameters reduced to that with n parameters derived from the SVD of the covariate matrix X. Therefore, the statistical inference on the original parameter turns into the superfactors. Let pi=(pi1,…,piJ) denote the vector of probabilities associated with the assignment of the ith sample into categories 1,…,J, where pij denote the probability that a sample falls into category j. From (2.1) and (2.3), it follows that
(2.4)pij=∫θj-1θjϕ(z-liγ)dz=Pr(θj-1<Zi<θj)=Φ(θj-liγ)-Φ(θj-1-liγ),
where ϕ(·) and Φ(·) denote the probability density function and the cumulative density function of the standard normal distribution, respectively, and li is the ith row of matrix L=FD. Let y=(y1,…,yn) denote the vector of responses observed for all samples. Then, the probability of observing data y is given as follow:
(2.5)Pr(y∣pi)=∏i=1n∏j=1JpijI{yi=j}.
From (2.4) and (2.5), the log likelihood function for (γ,θ) can be written as
(2.6)ln⁡L(γ,θ)=∑i=1n∑j=1JI{yi=j}ln⁡[Φ(θj-liγ)-Φ(θj-1-liγ)].


## 2.2. Model Fitting with Maximum Likelihood Estimation

The maximum likelihood estimates (MLEs) of the superfactors,γ, in (2.3) can be obtained by the iteratively reweighted least squares (IRLSs) procedure [6] using the log likelihood function for (γ,θ) in (2.6). The procedure can be briefly described as follows. Let η denote the vector of all model parameters, that is, η=(θ2,…,θJ-1,γ1,…,γn-J+2). Note that θ1 and θJ are not included in this vector because their values are assumed to be 0 and ∞, respectively, for the purpose of model identifiability. Also note that the (J-2) smallest singular values together with their corresponding factors are dropped from the parameters since the number of parameters must not exceed the number of samples. Assuming thatJ=4, define that
(2.7)Ci=[100-1100-1100-1],Li=[100-li010-li001-li],
and ℋi=diag(fi1,…,fiJ-1), where fij denotes the derivative of the standard normal cumulative distribution function at θj-liγ. Take Wi=diag(pi), where pi is the J×1 vector of probabilities that the ith individual falls in each category, that is, pi=(pi1,…,piJ)′, and let 𝒩i be a J×1 vector of observation, that is, 𝒩i=(I{yi=1},…,I{yi=J})′. After initialization of all elements, the iteration s+1(s=1,2,…) can be written as
(2.8)ηs+1=ηs+(∑i=1nLi′Hi(s)Ci′Wi-1(s)CiHi(s)Li)-1∑i=1nLi′Hi(s)Ci′Wi-1(s)(Ni-pi(s)).
The MLE of η can be found by performing the process recursively until the change between ηs+1 and ηs is negligible.


## 2.3. General Solution for the Original Parameters

We have discussed how to estimate the superfactor (γ) in (2.3) thus far. Since the primary interest is to find SNPs that are significantly associated with a disease, it is necessary to transform the superfactor (γ) to the original parameters (β) in (2.1). The equation γ=A′β in (2.3) can be utilized for the transformation even though A is m×n nonsquare matrix. As discussed by Graybill [7], the unique solution for β can be achieved by taking the generalized inverse matrix of A′ as A since A′A=In. Therefore, the unique solution for SNP effect (β) can be calculated by β=Aγ.


## 2.4. Selection of Significant SNPs

Finding significant SNPs is the same as testing if each SNP effect(βi,i=1,…,m) is statistically significant, that is, testing the hypothesis: H0: βi=0 versus H1: βi≠0, i=1,…,m. The simple method is to use Wald’s test statistic, which forms (β̂-β)/se(β̂) and assumes a normal distribution. However, when m≫n, it is hard to calculate se(β̂) directly from the data. We therefore utilized permutation test to select significant SNPs. The rationale behind the test is that, under the null hypothesis, the estimate of β obtained from the raw (unpermuted) data is similar to the estimate of β obtained from the permuted data. That is, the difference between two estimates is closed to zero under H0. With this idea, we can construct Wald’s test statistic as follows. Let β̂i (i=1,…,m) be the estimate of the ith SNP effect from the raw data and β̂ik (k=1,…,K) be the estimate of the ith SNP effect from the kth-permuted data. Let us define βidk as the difference between β̂i and β̂ik, that is, βidk=β̂i-β̂ik. Then, Wald’s test statistic can be as follows:
(2.9)Λi=β¯idse(β¯id),i=1,…,m,
where β¯id is the sample mean of βidk’s, which is β¯id=(1/K)∑k=1Kβidk, and se(β¯id) is the standard error of β¯id. Under the null hypothesis, the statistic Λi defined in (2.9) follows approximately standard normal distribution when k is large. P value for rejecting the null hypothesis at a significance level α=0.05 can be utilized to identify significant SNPs.


## 2.5. Application of the Multinomial Probit Model with SVD

### 2.5.1. Simulated Multinomial Ordinal Data

The validity of the proposed method was evaluated using simulated data sets. The procedure of data generation was composed of three steps: generating genotype data with certain genetic model, generating the latent variable, and defining the disease status variable by applying the predefined bin boundaries. The brief scheme of each step is as follows: we first generated 10 sets of the simulated genotype data under an additive genetic model, each set consists of 100 samples and 1000 SNPs. From (2.1), we can notice that the latent variable (zi) consists of two parts: the expected value (xiβ) and the random error (ϵi). In order to generate the expected value, we assumed that, for each sample, 9 out of the 1000 SNPs (every 101th SNP, except the last one) contribute to disease status xiβ=β1·SNP101+β2·SNP201+⋯+β1·SNP801+β2·SNP901, where β1 and β2 are set as −1 and 1, respectively. Hence, the latent variable can be obtained from the sum of the expected values (xiβ) and the random error generated from standard normal distribution. We then generated disease status variable (yi) assuming 3 disease development stages. Therefore, when applying the proposed method to the simulated data sets, we would expect 9 strong signals corresponding to each of the 9 disease-associated SNPs. We also compared results obtained from the proposed method with that from single SNP analysis with multinomial ordinal probit model.


## 2.5.1. Simulated Multinomial Ordinal Data

The validity of the proposed method was evaluated using simulated data sets. The procedure of data generation was composed of three steps: generating genotype data with certain genetic model, generating the latent variable, and defining the disease status variable by applying the predefined bin boundaries. The brief scheme of each step is as follows: we first generated 10 sets of the simulated genotype data under an additive genetic model, each set consists of 100 samples and 1000 SNPs. From (2.1), we can notice that the latent variable (zi) consists of two parts: the expected value (xiβ) and the random error (ϵi). In order to generate the expected value, we assumed that, for each sample, 9 out of the 1000 SNPs (every 101th SNP, except the last one) contribute to disease status xiβ=β1·SNP101+β2·SNP201+⋯+β1·SNP801+β2·SNP901, where β1 and β2 are set as −1 and 1, respectively. Hence, the latent variable can be obtained from the sum of the expected values (xiβ) and the random error generated from standard normal distribution. We then generated disease status variable (yi) assuming 3 disease development stages. Therefore, when applying the proposed method to the simulated data sets, we would expect 9 strong signals corresponding to each of the 9 disease-associated SNPs. We also compared results obtained from the proposed method with that from single SNP analysis with multinomial ordinal probit model.


## 2.6. Mexican-American Coronary Artery Disease (MACAD) Study

We also applied the proposed method to study sample recruited through the Mexican-American Coronary Artery Disease (MACAD) study [8, 9]. The study population consists of probands who are Mexican American aged between 45 and 75 with coronary artery disease: spouses of probands, adult offspring (≥18), and their spouses. For the offspring generation, we performed oral glucose tolerance test and genotyped 132 SNPs in 32 genes selected based on a prior relationship to insulin physiology. The goal of the study herein was to identify genes involved in the development of IGT and/or IFG, where IGT was defined as a 2 hr glucose level between 140 and 199 mg/dL and IFG defined as a fasting glucose level between 100 and 125 mg/dL. In order to identify and compare genes affecting the development of IGT and/or IFG, we generated two study samples, for which each sample has 3 disease stages (D1) both 2 hr and fasting glucoses normal (N/N)(n1=60), IGT only (IGT/N) (n2=31) and IGT and IFG (IGT/IFG) (n3=15) (D2) both 2 hr and fasting glucoses normal (N/N)(n1=60), IFG only (N/IFG) (n2′=34) and IGT and IFG (IGT/IFG) (n3=15).


## 3. Results and Discussion

### 3.1. Simulated Multinomial Data

Figure1 summarizes the results of association analyses when applying the multinomial ordinal probit model with SVD to the simulated data sets. All numbers shown in the figures are the average of the estimates obtained from the 10 simulated data sets. As mentioned previously, we expected 9 strong signals corresponding to the 9 SNPs designed to be associated with disease development when generating the simulated data sets, and 9 were observed in our analysis. Similar results from the single SNP analysis were shown in Figure 2.Analysis of the simulated data sets with multinomial ordinal probit model with SVD.
(a)
Estimates of SNP effects by the proposed method(b)
P values of testing SNP effects (-log10(Pvalue))Analysis of the simulated data sets with single SNP analysis.
(a)
Estimates of SNP effects by single SNP analysis(b)
P values of testing SNP effects (-log10(Pvalue))Figure1(a) summarizes MLEs of SNP effects calculated with the multinomial ordinal probit model with SVD. The figure shows that almost all MLEs except 9 were between −0.1 and 0.1, while there were 9 large MLEs (4 around 0.3, 5 around −0.3) corresponding to the 9 SNPs contributed to disease status. Figure 1(b) gives P values in -log10 scale for testing SNP effects. The line in Figure 1(b) corresponds to significance level α=0.05. 9 SNPs were clearly separated from the rest and had -log10(Pvalue)>1.3.Figure2(a) summarizes MLEs of SNP effects obtained by the single SNP analysis and shows that no signal was strong enough to be distinguished from all other signals. The P values are given in Figure 2(b) in -log10 scale. SNP501 in the middle of the figure had a relatively strong signal compared to all others. However, the -log10(P value) was much less than 1.3, which corresponds to significance level α=0.05. Thus, no SNPs were identified as statistically significant from the single SNP analysis method. In contrast to the fact that no SNP was identified as statistically significant by the single SNP analysis, the multinomial ordinal probit model with SVD method was able to identify all 9 SNPs contributing to disease status as statistically significant at significance level α=0.05. These results indicated that the proposed method should be reliable for the analysis of large-scale genome-wide association data that have polytomous ordinal responses when m≫n.


### 3.2. Mexican-American Coronary Artery Disease (MACAD) Study

We analyzed the data sets D1 and D2 (see methods) generated from a subsample of subjects recruited through a coronary artery disease proband in the Mexican-American Coronary Artery Disease Project as described in the method section, using both the multinomial ordinal probit model with SVD method and the single SNP analysis method.Figure3 summarizes the analysis results with the data D1 (N/N-IGT/N-IGT/IFG) using the single SNP analysis. Figure 3(a) gives MLEs of SNP effects. Figure 3(b) plots P values of association analysis in -log10 scale. With Sidak correction, which is often used to correct multiple testing problem, the adjusted significance level should be 1-(1-α)1/m, where α is significance level, and m represents the number of tests. Thus, the corrected -log10(P value) threshold for significance level α=0.05 is 3.4, which corresponds to the line in Figure 3(b). We applied the adjusted significance level to the P values in Figure 3(b) since the P values are before correcting multiple testing problem. No SNP was identified as statistically significant (Figure 3(b)).Analysis of genes for IGT/IFG through IGT pathway (Data Set D1) with single-SNP analysis.
(a)
Estimates of parameters by single-SNP analysis(b)
P values of testing SNP effects (-log10(Pvalue))The data set D2 (N/N-N/IFG-IGT/IFG) was analyzed with the same method, and analysis results are given in Figure4. Figure 4(a) plots MLEs of the SNP effects. Since P values in Figure 4(b) are before the the multiple testing correction, we used 3.4 as the -log10(P value) threshold corresponding to 0.05 significance level after the multiple testing correction. Two SNPs corresponding to SORC1 and SORT1 were found significant.Analysis of genes for IGT/IFG through IFG pathway (Data Set D2) with single SNP-analysis.
(a)
Estimates of parameters by single SNP-analysis(b)
P values of testing SNP effects (-log10(Pvalue))We then also analyzed D1 and D2 with the multinomial ordinal probit model with SVD method. Figure5 summarizes the analysis results for data D1. Figure 5(a) plots MLEs of SNP effects. Figure 5(b) plots P values in -log10 scale for testing SNP effects. The multiple testing correction does not need to be applied now since the method tests all SNPs simultaneously. With the 1.3 P value threshold, which corresponds to 0.05 significance level, we identified that 8 out of the 32 candidate genes (SORCS1, AMPD1, PPARα, AMPD2, PRKAA2, C5, TCF7L2, and ITR) were associated with the disease path defined in D1.Analysis of Genes for IGT/IFG through IGT pathway (Data Set D1) with multinomial ordinal probit model with SVD.
(a)
Estimates of parameters by the proposed method(b)
P values of testing SNP effects (-log10(Pvalue))The multinomial ordinal probit model with SVD method was applied to data set D2 as well. The results are shown in Figure6. In Figure 6(a), MLEs of the SNP effects were summarized. Figure 6(b) plots the P values in -log10 scale for testing the SNP effects. It showed that 11 SNPs corresponding to 9 out of 32 candidate genes (SORCS1, AMPD1, PPARα, CAPN10, IL4, NOS3, CD14, GCG, and SORT1) have -log10(P value) greater than the 1.3 P value threshold. From the analyses of D1 and D2, we found that SNPs in 3 genes (SORCS1, AMPD, and PPARα) were associated with both IGT and IFG; SNPs in 5 genes (AMPD2, PRKAA2, C5, TCF7L2, and ITR) were associated with IGT only; SNPs in 6 genes (CAPN, IL4, NOS3, CD14, GCG, and SORT1) were associated with IFG only. These results suggest that IGT and IFG may indicate different pathways to diabetes, with different genetic determinants.Analysis of genes for IGT/IFG through IFG pathway (Data Set D2) with multinomial ordinal probit model with SVD.
(a)
Estimates of parameters by the proposed method(b)
P values of testing SNP effects (-log10(Pvalue))Thus, using both simulated data and a real study sample, we demonstrated that multinomial ordinal probit model with SVD method can be utilized to identify associated markers involved in disease development when multidisease stages are considered. For relatively small size of data set used in the paper, which is 100 samples and 1000 SNPs for the simulation study, the computation took about less than 10 minutes to complete. However, the computation time might be a concern when applying this method to large data set, such as GAWS with millions of SNPs and thousands of samples.


## 3.1. Simulated Multinomial Data

Figure1 summarizes the results of association analyses when applying the multinomial ordinal probit model with SVD to the simulated data sets. All numbers shown in the figures are the average of the estimates obtained from the 10 simulated data sets. As mentioned previously, we expected 9 strong signals corresponding to the 9 SNPs designed to be associated with disease development when generating the simulated data sets, and 9 were observed in our analysis. Similar results from the single SNP analysis were shown in Figure 2.Analysis of the simulated data sets with multinomial ordinal probit model with SVD.
(a)
Estimates of SNP effects by the proposed method(b)
P values of testing SNP effects (-log10(Pvalue))Analysis of the simulated data sets with single SNP analysis.
(a)
Estimates of SNP effects by single SNP analysis(b)
P values of testing SNP effects (-log10(Pvalue))Figure1(a) summarizes MLEs of SNP effects calculated with the multinomial ordinal probit model with SVD. The figure shows that almost all MLEs except 9 were between −0.1 and 0.1, while there were 9 large MLEs (4 around 0.3, 5 around −0.3) corresponding to the 9 SNPs contributed to disease status. Figure 1(b) gives P values in -log10 scale for testing SNP effects. The line in Figure 1(b) corresponds to significance level α=0.05. 9 SNPs were clearly separated from the rest and had -log10(Pvalue)>1.3.Figure2(a) summarizes MLEs of SNP effects obtained by the single SNP analysis and shows that no signal was strong enough to be distinguished from all other signals. The P values are given in Figure 2(b) in -log10 scale. SNP501 in the middle of the figure had a relatively strong signal compared to all others. However, the -log10(P value) was much less than 1.3, which corresponds to significance level α=0.05. Thus, no SNPs were identified as statistically significant from the single SNP analysis method. In contrast to the fact that no SNP was identified as statistically significant by the single SNP analysis, the multinomial ordinal probit model with SVD method was able to identify all 9 SNPs contributing to disease status as statistically significant at significance level α=0.05. These results indicated that the proposed method should be reliable for the analysis of large-scale genome-wide association data that have polytomous ordinal responses when m≫n.


## 3.2. Mexican-American Coronary Artery Disease (MACAD) Study

We analyzed the data sets D1 and D2 (see methods) generated from a subsample of subjects recruited through a coronary artery disease proband in the Mexican-American Coronary Artery Disease Project as described in the method section, using both the multinomial ordinal probit model with SVD method and the single SNP analysis method.Figure3 summarizes the analysis results with the data D1 (N/N-IGT/N-IGT/IFG) using the single SNP analysis. Figure 3(a) gives MLEs of SNP effects. Figure 3(b) plots P values of association analysis in -log10 scale. With Sidak correction, which is often used to correct multiple testing problem, the adjusted significance level should be 1-(1-α)1/m, where α is significance level, and m represents the number of tests. Thus, the corrected -log10(P value) threshold for significance level α=0.05 is 3.4, which corresponds to the line in Figure 3(b). We applied the adjusted significance level to the P values in Figure 3(b) since the P values are before correcting multiple testing problem. No SNP was identified as statistically significant (Figure 3(b)).Analysis of genes for IGT/IFG through IGT pathway (Data Set D1) with single-SNP analysis.
(a)
Estimates of parameters by single-SNP analysis(b)
P values of testing SNP effects (-log10(Pvalue))The data set D2 (N/N-N/IFG-IGT/IFG) was analyzed with the same method, and analysis results are given in Figure4. Figure 4(a) plots MLEs of the SNP effects. Since P values in Figure 4(b) are before the the multiple testing correction, we used 3.4 as the -log10(P value) threshold corresponding to 0.05 significance level after the multiple testing correction. Two SNPs corresponding to SORC1 and SORT1 were found significant.Analysis of genes for IGT/IFG through IFG pathway (Data Set D2) with single SNP-analysis.
(a)
Estimates of parameters by single SNP-analysis(b)
P values of testing SNP effects (-log10(Pvalue))We then also analyzed D1 and D2 with the multinomial ordinal probit model with SVD method. Figure5 summarizes the analysis results for data D1. Figure 5(a) plots MLEs of SNP effects. Figure 5(b) plots P values in -log10 scale for testing SNP effects. The multiple testing correction does not need to be applied now since the method tests all SNPs simultaneously. With the 1.3 P value threshold, which corresponds to 0.05 significance level, we identified that 8 out of the 32 candidate genes (SORCS1, AMPD1, PPARα, AMPD2, PRKAA2, C5, TCF7L2, and ITR) were associated with the disease path defined in D1.Analysis of Genes for IGT/IFG through IGT pathway (Data Set D1) with multinomial ordinal probit model with SVD.
(a)
Estimates of parameters by the proposed method(b)
P values of testing SNP effects (-log10(Pvalue))The multinomial ordinal probit model with SVD method was applied to data set D2 as well. The results are shown in Figure6. In Figure 6(a), MLEs of the SNP effects were summarized. Figure 6(b) plots the P values in -log10 scale for testing the SNP effects. It showed that 11 SNPs corresponding to 9 out of 32 candidate genes (SORCS1, AMPD1, PPARα, CAPN10, IL4, NOS3, CD14, GCG, and SORT1) have -log10(P value) greater than the 1.3 P value threshold. From the analyses of D1 and D2, we found that SNPs in 3 genes (SORCS1, AMPD, and PPARα) were associated with both IGT and IFG; SNPs in 5 genes (AMPD2, PRKAA2, C5, TCF7L2, and ITR) were associated with IGT only; SNPs in 6 genes (CAPN, IL4, NOS3, CD14, GCG, and SORT1) were associated with IFG only. These results suggest that IGT and IFG may indicate different pathways to diabetes, with different genetic determinants.Analysis of genes for IGT/IFG through IFG pathway (Data Set D2) with multinomial ordinal probit model with SVD.
(a)
Estimates of parameters by the proposed method(b)
P values of testing SNP effects (-log10(Pvalue))Thus, using both simulated data and a real study sample, we demonstrated that multinomial ordinal probit model with SVD method can be utilized to identify associated markers involved in disease development when multidisease stages are considered. For relatively small size of data set used in the paper, which is 100 samples and 1000 SNPs for the simulation study, the computation took about less than 10 minutes to complete. However, the computation time might be a concern when applying this method to large data set, such as GAWS with millions of SNPs and thousands of samples.



---

*Source: 419832-2012-05-29.xml*