A Powerful Method To Test Associations Between Ordinal Traits and Genotypes

The methods commonly used to test the associations between ordinal phenotypes and genotypes often treat either the ordinal phenotype or the genotype as continuous variables. To address limitations of these approaches, we propose a model where both the ordinal phenotype and the genotype are viewed as manifestations of an underlying multivariate normal random variable. The proposed method allows modeling the ordinal phenotype, the genotype and covariates jointly. We employ the generalized estimating equation technique and M-estimation theory to estimate the model parameters and deduce the corresponding asymptotic distribution. Numerical simulations and real data applications are also conducted to compare the performance of the proposed method with those of methods based on the logit and probit models. Even though there may be potential limitations in Type I error rate control for our method, the gains in power can prove its practical value in case of exactly ordinal phenotypes.

With the development of high throughput biologic technology, increasingly more genotypes and data with complex traits have been generated and deposited in public databases. It is urgently required to develop new statistical testing methods to investigate the associations between these and extract useful information to understand the underlying occurrence and development mechanisms of diseases and traits.
Genome-wide association studies aim to identify associations between phenotypes and genotypes. In these studies, genotypes are often treated as predictors and phenotypes as outcomes. If the phenotype of interest is continuous, then the classic linear regression model is commonly employed. When the phenotype is ordinal, the multinomial logit model (McCullagh 1980; or ordered probit model (Daykin and Moffatt 2002;Wang 2014) should be recommended. All these models regress phenotype values or their distribution-based transformations on genotypes, with the assumptions that genotype values are continuous (Korse et al. 2009;Bedogni et al. 2010) and the probability of having a disease increases linearly with the genotype value. However, the continuity assumption on genotype values and the linearity assumption between a phenotype and genotype are difficult to verify in practice. If these two assumptions are violated, the corresponding Wald testing statistics may severely decrease in power. To overcome this, some researchers treated genotypes as ordinal variables and reversed the regression process by regressing genotypes on phenotypes (O'Reilly et al. 2012). When a phenotype is a continuous variable, this new method is indeed useful for removing or relaxing the continuity and linearity assumptions. However, this does not work when a phenotype is exactly ordinal, such as in the above-mentioned example of liver steatosis. Therefore, we propose a new method to deal with this problem.
In this work, we treat genotypes as ordinal variables and propose a new procedure to assess the association between an ordinal phenotype and ordinal genotype after adjusting for covariates. Rather than regressing the phenotype on the genotype or regressing the genotype on the phenotype using existing methods, we jointly model the phenotype and genotype by introducing a latent variable following a multivariate normal distribution. The phenotype and genotype are regarded as manifestation values of the latent variable. The relationships between phenotypes, genotypes, and covariates of interest are elaborately described by the covariance matrix. Taking advantage of the framework of generalized estimation equations (Hanley et al. 2003;Zhang et al. 2014) and M-estimation theory (Huber 1981;Stefanski and Boos 2002), we construct a Wald test statistic for an equivalent transformation of the original null hypothesis, and prove that it asymptotically follows the standard normal distribution under the null hypothesis. Numerical simulations are conducted to compare the proposed method with other methods. Our simulation results show that the proposed method can suitably maintain Type I error control and may achieve considerable statistical power compared to existing methods in various scenarios. Finally, we apply the proposed method to anticyclic citrullinated protein antibody data for rheumatoid arthritis studies, to further demonstrate its performance.

Notations
Let random variables G and Y denote a genotype and ordinal phenotype, respectively, and Z ¼ ðZ 1 ; Z 2 ; ⋯ ; Z k22 Þ t be a ðk 2 2Þ-dimensional continuous covariate of interest, with k $ 3, where t denotes the transpose of a matrix or vector. Without loss of generality, we assume that there are two alleles at a genetic locus, with one being the risk allele and the other being the reference allele. The value of the random variable G represents the number of the risk allele at a locus, which means that G can take three values: 0, 1, or 2. Suppose that Y takes m ordinal values: 1; 2; ⋯ ; m. The null hypothesis states that after adjusting for covariates, the phenotype is not related to the genotype, i.e., the phenotype and genotype are conditionally independent given a set of covariates, which can be denoted by H 01 : PrðG; YjZÞ ¼ PrðGjZÞPrðYjZÞ: (1) Assume that n subjects are enrolled in the genetic association study. Further, let fg 1 ; g 2 ; ⋯ ; g n g, fy 1 ; y 2 ; ⋯ ; y n g, and fz 1 ; z 2 ; ⋯ ; z n g be the n observations of G, Y, and Z corresponding to the n subjects, respectively.
Equivalent statement of H 01 by introducing a latent variable U We assume that the combined vector ðG; Y; Z t Þ t is generated from a k-dimensional random variable U ¼ ðU 1 ; U 2 ; ⋯ ; U k Þ t following a multivariate normal distribution, with mean vector 0 k ¼ ð0; 0; ⋯ ; 0Þ t k · 1 (the k-dimensional column vector with all units being zero) and covariance matrix D, where D ¼ ðd k1k2 Þ k·k and all its diagonal elements are equal to one. That is, ðG; YÞ represents the manifestation of U 01 ¼ ðU 1 ; U 2 Þ and Z ¼ ðU 3 ; ⋯ ; U k Þ t . We rewrite D in the partitioned matrix form where D 11 is a 2 · 2 matrix. Then, Z follows a multivariate normal distribution with mean vector 0 k22 ¼ ð0; 0; ⋯ ; 0Þ t ðk22Þ · 1 and covariance matrix D 22 . Then, G and Y can be obtained as follows: Based on the theory of the conditional normal distribution, we have that We define z 1 ¼ ðd 13 ; d 14 ; ⋯ ; d 1k Þ t ; z 2 ¼ ðd 23 ; d 24 ⋯ ; d 2k Þ t : Then, the conditional covariance matrix above can be further expressed as Now, denote d 0 ¼ ðd 12 ; z t 1 ; z t 2 Þ t , .ðd 0 Þ ¼ d 12 2 z t 1 D 21 22 z 2 . By introducing the latent variable U and taking advantage of its distribution property, we can state that the original hypothesis H 01 is exactly equivalent to H 02 : .ðd 0 Þ ¼ 0: Proposed statistical test In this subsection, we construct a test statistic to test H 02 based on the generalized estimating equation technique and M-estimation theory.
Recall that the joint distribution of U 1 and Z t is Hence, the conditional distribution of U 1 given Z is Similarly, the conditional distribution of U 2 given Z is It should be noted that the marginal density function of Z and the joint density functions of each pair of variables among G; Y, and Z are as follows: where IðEÞ is an indicator function, i.e., IðEÞ is one if the event E holds and zero otherwise. The unknown parameters D 22 ; fj j ; j ¼ 1; 2g; fh r ; r ¼ 1; ⋯ ; m 2 1g and d 0 can be estimated via the following procedures.
First, based on the marginal density function (11) of Z, we have the likelihood function By maximizing LðD 22 Þ on D 22 , we can obtain the MLE of D 22 (denoted byD 22 ), which is the sample covariance matrix of the observed data fz 1 ; z 2 ; ⋯ ; z n g.
Second, in our model both G and Y are ordinal variables, whose realizations are determined by the intervals in which the values of two standard normal random variables U 1 and U 2 may fall in, respectively. Specifically, we employ the distribution properties of G and Y to intuitively estimate fj j ; j ¼ 1; 2g and fh r ; r ¼ 1; ⋯ ; m 2 1g, respectively. We define , and v r ¼ #fyi¼r;i¼1;2 ⋯ ;ng n for r ¼ 1; 2; ⋯ ; m based on the observed data fg 1 ; g 2 ; ⋯ ; g n g and fy 1 ; y 2 ; ⋯ ; y n g of G and Y, respectively. Recall that j 0 ¼ h 0 ¼ 2N. Then, we can estimate j j asĵ j , j ¼ 1; 2 and h r asĥ r , r ¼ 1; ⋯ ; m 2 1 by solving the following equations: where FðÁÞ is the cumulative distribution function of the standard normal variable. The parameter d 0 can be estimated using the generalized estimating equation technique and M-estimation theory. First, let The function vector consisting of the first-order partial derivatives of hðd 0 Þ with respect to each parameter in d 0 is Then, the estimatord 0 ¼ ðd 12 ;ẑ t 1 ;ẑ t 2 Þ t of d 0 is the root of the following generalized estimating equation After estimating all the unknown parameters, the estimate of .ðd 0 Þ can be expressed as .ðd 0 Þ ¼d 12 2ẑ t 1D 21 22ẑ2 . To construct a Wald-type test statistic, we need to derive the asymptotical distribution of .ðd 0 Þ. Based on the classical M-estimation theory,d 0 asymptotically follows a multivariate normal distribution. That is, where @uðgi;yi;zi;d0Þ @d t 0 : According to the delta method, .ðd 0 Þ is also asymptotically multivariate normal. Namely, where Now, we can propose a new Wald test statistic for the null hypothesis H 02 as follows: Nð0; 1Þ as n/N:

Data availability
The authors affirm that all data necessary for confirming the conclusions of the article are present within the article, figures, and tables. Supplemental material available at FigShare: https://doi.org/10.25387/ g3.8226650.

Simulation results
In this subsection, we present a series of simulation studies to investigate the performance of our proposed latent variable model (abbreviated as lvm), and compare it with the probit and logit models, which both regressing an ordinal phenotype on a genotype. We compared these under multiple simulation scenarios, so that different modeling assumptions would be favored. Two types of data generation mechanism were considered throughout our simulation studies: (i) generating data from a multivariate normal random variable (named the ND mechanism) and (ii) generating data under the proportional odds model (named the PO mechanism). In addition, three genetic models (co-dominant, dominant, and recessive models) were considered. The specifics of our simulation data generation scenarios are as follows. For simplicity, the dimension of the covariate Z is one, and the number of levels for a phenotype Y is set to five. Under the ND mechanism, the three-dimensional latent variable U was generated from a multivariate normal distribution, with the genotype G being a manifestation of U 1 , the ordinal phenotype Y being a manifestation of U 2 , and the covariate Z being equal to U 3 . Each marginal distribution of U 1 ; U 2 ; and U 3 followed a standard normal distribution. Note that the distribution of G varied according to the type of the true genetic model. Let A (major allele) and a (minor allele) denote two alleles at the single biallelic locus corresponding to the genotype G. Under the Hardy-Weinberg equilibrium (HWE) conditions, the expected genotype frequencies of G being AA, Aa, and aa would be ð12pÞ 2 ; 2pð1 2 pÞ, and p 2 , respectively, with the minor allele frequency (MAF) p taking on values from f0:1; 0:2; 0:3; 0:4; 0:5g. If the HWE assumption did not hold, we directly set three kinds of multinomial distribution (P) for (AA, Aa, aa) with different parameter structures P 1 ¼ ð0:2; 0:3; 0:5Þ; P 2 ¼ ð0:2; 0:5; 0:3Þ, and P 3 ¼ ð0:2; 0:7; 0:1Þ. When a co-dominant model was assumed, the three genotypes AA, Aa, and aa were coded as 0, 1, and 2, respectively. In a model where a dominant effect was assumed, the genotype AA was coded as 0 while both Aa and aa were coded as 1. Accordingly, scores of 0 for both AA and Aa and 1 for aa were employed in a recessive model. In addition, under each genetic model the probabilities of Y being 1; 2; 3; 4, and 5 were always 0:7; 0:1; 0:08; 0:07, and 0.05, respectively. We set the covariance matrix Σ of U as to investigate the Type I error rate, and let Σ be to compare the statistical power under different alternatives, depending on the parameter u, whose range was the set f20:2; 2 0:1; 0:1; 0:2g. Under the PO mechanism, the ordinal phenotype Y was related to the genotype G and covariate Z through the proportional odds model. Specifically, the distribution of G under different genetic models would remain the same as that under the ND mechanism, regardless of whether the HWE held. In this case, Z still followed the standard normal distribution, and Y was generated with five levels f1; 2; 3; 4; 5g using the proportional odds model where j ¼ 1; ⋯ ; 4, ða 1 ; a 2 ; a 3 ; a 4 Þ ¼ ð0:85; 1:39; 1:99; 2:94Þ, ðb 1v ; b 2v ; b 3v Þ ¼ v · ð0; 0:15; 0:3Þ, v ¼ 0; 5. It should be noted that PrðY # 5Þ ¼ 1. When v ¼ 0, the null hypothesis was true, and the derived simulation data were used to compare the Type I error rates for the three models of interest, i.e., the lvm, probit, and logit models.
When v ¼ 5, we explored the power of these three models under different alternatives. As previously described, we considered 5 · 5 · 3 þ 2 · 5 · 3þ 5 · 3 · 3 þ 2 · 3 · 3 ¼ 168 scenarios. For each simulation scenario, we generated 1000 datasets, each consisting of 300 subjects. P-values were calculated for each dataset using the three respective models. The nominal level of the tests was set to 0.05, and all simulations were performed using the R language (https://www.r-project.org/). The empirical Type I error rates and power estimates were calculated using the percentage of rejection in each scenario. The results are presented side-by-side in Tables 1-4. Note that the blanks (marked with -) in these four tables are a result of unavailability of the simulation data under the corresponding scenarios. The reason for this is that in these parameter setting conditions, the mean number for G amounting to 1 is three, which can easily lead to samples with all the G values being 0, such that none of the three considered models can be applied. Table 1 presents the results for the ND mechanism when the HWE holds. The first five rows suggest that all three methods can control the Type I error rate at the nominal level of 0.05 under the three different genetic models. Furthermore, the remaining rows show that the lvm model enhances the statistical power over the probit and logit models. In some cases, the power gain can be as high as 0.07 to 0.1. For example, when u ¼ 2 0:2 and p ¼ 0:1, the empirical powers for the probit, logit, and lvm models are 0.514, 0.481, and 0.557 under the co-dominant genetic model, and are 0.503, 0.459, and 0.544 under the dominant genetic model, respectively. When the true genetic model is recessive and u ¼ 2 0:1 and p ¼ 0:3, the power estimates for the probit, logit, and lvm models are 0.126, 0.094, and 0.166, respectively.
The results of the three methods under the ND mechanism when the HWE does not hold are presented in Table 2. We can observe that the proposed method achieves a greater power than the other two methods, even though the distribution of the genotype G does not satisfy the HWE conditions. Specifically, the power gain can be as high as 0.07 when u ¼ 0:2 and the distribution of G is P 2 .
The corresponding results when the data are generated under the PO mechanism and the HWE is assumed are displayed in Table  3. It follows that all three models can control the Type I error rates under the null hypothesis. Even though the data are generated using the proportional odds model, the lvm method still performs better than the other two methods in detecting an alternative hypothesis, and can achieve a power gain of up to 0.058 in some cases, such as the setting with p ¼ 0:3 for the recessive genetic model. It is worth noting that the advantage of the proposed lvm method is more obvious under the recessive model.
The simulation results under the PO mechanism when the HWE does not hold are presented in Table 4. We observe that the advantage of the proposed lvm method is not as significant as that when the HWE holds, but the model is superior to the probit model in all scenarios. Moreover, the logit model has a slightly greater power under the co-dominant and dominant models, while the lvm method outperforms it in the recessive model.
n From all of the four tables, it can be seen that the proposed method might have potential limitations in controlling Type I error rates in a few situations, while the power gains in almost all of simulation scenarios indeed indicate its efficiency for practical applications.
Application to anticyclic citrullinated protein antibody data for rheumatoid arthritis study It is well known that rheumatoid arthritis (RA) is significantly associated with some genetic variants (Carlton et al. 2005;Ruiz-Larrañaga et al. 2016). The anticyclic citrullinated protein antibody (anti-CCP) can be an auxiliary diagnosis indicator for RA, and the specificity of anti-CCP lies between 87.8% and 96.4% (Coenen et al. 2007). Besides, the genomic region of 6p21.33 has been reported to be associated with RA (Zhang et al. 2009;. The aim here is to check whether the single nucleotide polymorphisms (SNPs) in the 6p21.33 region are associated with anti-CCP, taking advantage of the proposed lvm test.
Note that there are a total of 45 SNPs in the region of 6p21.33 according to the Genetic Analysis Workshop 16 Data, all of which meet the quality control rule of the MAF being more than 5%, the missing rate being smaller than 15%, and the least genotype frequency being no less than five. The anti-CCP measure takes four values, 1, 2, 3, and 4, and the number of subjects who have these four values were 1195, 103, 66, and 698, respectively. The total number of subjects was 2062. Five principal components coordinated by applying the multi-dimensional scaling method  to the 12747 population structure information SNPs ) were used to adjust for population stratification effects. Before conducting association analysis, we run chi-square tests to check whether HWE holds for each of these 45 SNPs in controls. The P-value results are summarized in Table 5. At a 0.05 level of significance, we can state that the HWE law holds for all SNPs in controls on the basis of their Bonferroni-corrected P-values. Then we apply the probit model, the logit model, and the lvm model to these 45 SNPs to test their association with anti-CCP in sequence. The results are presented in Table 5. It shows that after Bonferroni correction, the SNPs rs2246986, rs3093998, rs2071596, and rs2844509 were found to be significant under the probit and logit models, while the SNPs rs2516398, rs2844494, rs3130637, rs3093993, rs3095227, rs2259435, and rs2844509 were identified as significant using the lvm method. Though these two groups of SNPs overlap at only one SNP rs2844509, each of other SNPs found by the probit (logit) model is physically close to one or two SNPs found by the lvm model. For example, the SNP rs2246986 of the first group is 677 kb away from rs2516398 on one side and is 1212 kb away from rs2844494 on the other side, while both of these two SNPs rs2516398 and rs2844494 are included in the second group. In addition, the SNP rs3093998 (in the first group) is 2971 kb away from rs3130637 (in the second group). The distances are so short that it is reasonable to infer that the SNPs rs2246986, rs2516398, and rs2844494 contain similar information. So do another two SNPs rs3093998 and rs3130637. In short, for detecting the association between anti-CPP and the genomic region 6p21.33, the proposed lvm method n is more powerful than the methods based on logit and probit models.

DISCUSSION
In this work, we have shown that the idea of treating a genotype variable as ordinal without assuming linearity can result in a more powerful and robust test, via introducing a joint multivariate normal distribution for the group of genotypes, traits, and covariates. Meanwhile, we have also demonstrated that the proposed lvm test can provide appropriate Type I error rates. The important strength of our method is that it does not make an assumption on the type of relationship between a phenotype and genotype; nor does it treat the genotype as a continuous variable. Rather, our approach only introduces a latent multivariate normal variable to characterize the relationship between the two, which is very reasonable, and generally considerably more useful.
Besides the simulations with respect to significance level of 0.05, we also conducted simulation studies with a lower significance level 0.005. The results are given in Supplementary Table S1. We found that the proposed method can reasonably control Type I error rates and achieve power gains at this lower significance level, similar to those results in Tables 1-4. It is worth mentioning that our proposed test model can also be applied to other situations where the outcome is continuous. In such a situation, we can still employ a joint multivariate normal distribution to model outcomes, genotypes, and covariates simultaneously. Even though the proposed method might have potential limitations in Type I error rate control in some situations, the power gains prove its efficiency in practical applications.
In population-based genetic association studies, hundreds of thousands of subjects are often enrolled to achieve optimal power. It is inevitable that there exists a population stratification effect in such large-scale studies. Not considering the effect of population stratification could lead to many false positive findings, and therefore adjusting for its effect represents the basis for conducting a genetic association analysis (Price et al. 2006;Li and Yu 2008). In this study, to characterize the influence of population stratification when investigating the relationships between ordinal traits and genotypes, we treat these effects as covariates in the proposed lvm method. The numerical results of our simulation studies and real data applications have demonstrated that the strategy is feasible and effective.