A Generative Approach to Joint Modeling of Quantitative and Qualitative Responses

In many scientific areas, data with quantitative and qualitative (QQ) responses are commonly encountered with a large number of predictors. By exploring the association between QQ responses, existing approaches often consider a joint model of QQ responses given the predictor variables. However, the dependency among predictive variables also provides useful information for modeling QQ responses. In this work, we propose a generative approach to model the joint distribution of the QQ responses and predictors. The proposed generative model provides efficient parameter estimation under a penalized likelihood framework. It achieves accurate classification for qualitative response and accurate prediction for quantitative response with efficient computation. Because of the generative approach framework, the asymptotic optimality of classification and prediction of the proposed method can be established under some regularity conditions. The performance of the proposed method is examined through simulations and real case studies in material science and genetics.


Introduction
Analyzing data with heterogeneous types of responses has been an important topic with broad applications. Such heterogeneous data often involve both quantitative and qualitative (QQ) responses. For example, Klein et al. (2019) described a human health study examining the risk factors of adverse birth outcomes, which contains a qualitative response "presence/absence of low birth weight" and a quantitative response "gestational age". In material science, the properties of a material are often characterized by QQ measures. As shown in the case study of Section 5 on Heusler compounds, two metrics, the mixing enthalpy (quantitative) and the global stability based on hull energy (qualitative) are used to determine the thermodynamic stability of a full Heusler compound. The QQ responses is a special case of "mixed outcomes" in the literature. In this paper, we focus on two types of mixed responses: the quantitative continuous response and the multi-class qualitative response.
In the literature of mixed outcomes, particularly on QQ responses, it has been known that overlooking the relationship between the QQ responses is inappropriate. Researches on the joint model of the QQ responses include early ones such as Fitzmaurice and Laird (1995); Moustaki and Knott (2000); Dunson (2000); Gueorguieva and Agresti (2001); Dunson (2003) and recent ones such as Deng and Jin (2015); Kürüm et al. (2016); Kang et al. (2018); Amini et al. (2018); Klein et al. (2019). These works mostly consider a joint regression model conditioned on the predictor variables. Based on our best knowledge, we can group them into two categories.
The first group of methods considers a conditional regression model, in which one type of the QQ responses is treated as the response of the model, and the other type is treated as the regressor. For instance, Fitzmaurice and Laird (1997) introduced a marginal regression model of quantitative response conditioned on the qualitative response. Song et al. (2009) used Gaussian copulas to integrate separate one-dimensional generalized linear models into a joint regression model for mixed outcomes. Lin et al. (2010) developed a conditional mixedeffects model to analyze clustered data containing QQ responses. These methods are suitable for data with a small number of predictor variables. To handle the high-dimensional input variables, Deng and Jin (2015) proposed a conditional model that encourages model sparsity through a constrained likelihood estimation. However, inferences and asymptotic properties of their method have not been explored due to the complicated constrained likelihood estimation. Kang et al. (2018) considered a Bayesian estimation for the conditional model of Deng and Jin (2015) to obtain proper inferences of model parameters. Nevertheless, their work is not designated for studying the asymptotic properties of the proposed estimator.
The second group of methods considers a continuous latent variable for the qualitative response, and then jointly models the latent variable and the quantitative response (Sammel et al., 1997;Dunson and Herring, 2005;Bello et al., 2012). For example, Gueorguieva and Agresti (2001) studied a probit model with a latent variable and developed a Monte Carlo expectationconditional maximization algorithm for parameter estimation. Klein et al. (2019) introduced the idea of latent variable into the framework of copula regressions, constructing a latent continuous representation of binary regression models. However, the use of latent variables often involves considerable computation in the parameter estimation. It also makes the investigation of theoretical properties difficult. Moreover, most of these works focus on the binary qualitative response and their model assumptions may not be easily extended to the multi-class qualitative response cases.
In this work, we propose a novel approach to jointly model the QQ responses based on the generative approach. The proposed generative model considers the joint distribution of the high-dimensional input variables, the quantitative responses, and the multi-class qualitative response. It is a very unique and different perspective from the existing literature and also brings advantages in both theoretical and computational aspects. The proposed method can accommodate multi-class qualitative response and multivariate quantitative responses with attractive theoretical properties. We call the proposed method GAQQ, a Generative Approach for QQ responses.
The key contributions of this work are summarized as follows. First, based on the generative model framework, we are able to establish the asymptotic properties of the proposed estimators with respect to both the classification accuracy of the qualitative response and the prediction accuracy of the quantitative response under some regularity conditions. Such conditions are commonly used in the regularized estimation framework (Shao et al., 2011;Zhao and Yu, 2006). The classification of the qualitative response enjoys the asymptotic optimality of the resulting linear discriminate classification rule. The mean squared error (MSE) of prediction for the quantitative response is as good as the optimal prediction under the Bayes risk. Second, an efficient procedure for parameter estimation is developed via the regularized log-likelihood function of the joint distribution of input variables and QQ responses. Specifically, we impose regularization on both the mean differences and the covariance matrix from the joint distribution to achieve sparsity for high-dimensional predictor variables. Third, the use of the generative approach leads to an effective prediction procedure by inferring the conditional distribution of QQ responses conditioned on the predictor variables. That is, the quantitative response is predicted through the property of conditional multivariate normal distribution, and the linear discriminant analysis (LDA) is employed for classification of the qualitative response. Fourth, the proposed generative model allows the parameters related to QQ responses to be mutually learned from each other, which is different from existing methods in which only modeling one type of QQ responses attempts to benefit from the information of the other type of QQ responses.
The remainder of this paper is organized as follows. Section 2 details the proposed method. The main theoretical results are presented in Section 3. Simulation and real data analysis are conducted in Sections 4 and 5, respectively. Section 6 concludes this work with some discussion. Technical proofs are in the Appendix.

The Proposed Model
Suppose that the variables of interest are denoted by (X, y, Z) where X = (X 1 , . . . , X p−1 ) ′ is a (p − 1) dimensional vector of predictor variables, y is a quantitative response variable and Z ∈ {1, 2} is a qualitative response variable. From a generative modeling perspective, we consider the data generation mechanism as p(X, y, Z) = p(X, y|Z)p(Z), indicating that data are from two classes G 1 and G 2 under (X, y)|Z. Assume that W = (X ′ , y) ′ follows multivariate norm distributions with different means for two classes, but sharing the same covariance matrix as follows (1) Denote the observed data w 1 , · · · , w n 1 , w n 1 +1 , · · · , w n 1 +n 2 with the first n 1 observations from G 1 and the rest n 2 observations from G 2 , where w i = (x ′ i , y i ) ′ , i = 1, 2, . . . , n 1 + n 2 . Let n = n 1 + n 2 . The log-likelihood function can be written as up to some constant, where C = Σ −1 is the inverse covariance matrix. Let π 1 and π 2 be the prior probability of w belonging to classes G 1 and G 2 , respectively. Hence, the LDA assigns where δ = µ 1 − µ 2 . Otherwise, w is classified to G 2 . The estimates of π 1 and π 2 are the empirical proportions of data from each class. The parameters µ 1 , µ 2 and C can be estimated by maximizing the log-likelihood function of (2).
For high-dimensional data when p ≥ n, the regularization is often needed to ensure the proper estimation of inverse covariance matrix C and mean difference δ. We thus propose to penalize C = (c ij ) 1≤i,j≤p and δ simultaneously, resulting in the following optimization where ||C|| 1 = i =j |c ij |, and |α| 1 = i |α i | with α i being the ith entry of vector α. Here λ 1 ≥ 0 and λ 2 ≥ 0 are two tuning parameters. By applying such regularization, the proposed model can encourage the sparse structures in C and δ at the same time. Note that similar spirits of regularizing both C and δ are used in several works on the LDA (Shao et al., 2011;Cai and Liu, 2012).
To estimate the parameters, we develop an iterative procedure to solve the sub-optimization problem with respect to C and δ respectively. Define δ 2 = (µ 1 − µ 2 )/2 as well as γ = (µ 1 + µ 2 )/2, then accordingly we have µ 1 = δ 2 + γ and µ 2 = γ − δ 2 . As a result, the optimization problem (4) is re-written as It is thus easy to obtain the maximum likelihood estimate of γ from (5) aŝ wherew = 1 n n i=1 w i is the overall mean. Then pluggingγ back into (5) yields (δ 2 ,Ĉ) = arg min In this manner, solving the optimization problem (4) is equivalent to solving the optimization problem (7). Next, we show that (7) can be decomposed as a graphical lasso model (Glasso) (Yuan and Lin, 2007;Deng and Yuan, 2009) in terms of C and a Lasso regression (Tibshirani, 1996) in terms of δ 2 with the other parameter fixed, such that these two parameters can be estimated iteratively. To be more precise, for a given value of δ 2 , the minimization problem (7) with respect to C is min C −n ln |C| + tr(CS) + λ 1 ||C|| 1 , It has the same form as the graphical lasso, which has been extensively studied in literature by Yuan and Lin (2007);Friedman et al. (2008); Lam and Fan (2009);Raskutti et al. (2008); Liu et al. (2020), and many others. On the other hand, when the inverse covariance matrix C is fixed, the minimization problem (7) regarding δ 2 becomes min δ 2 i∈G 1 which is equivalent to min δ 2 whereỹ = 1 2n 1 n 2 C 1/2 (n 2 i∈G 1 w i − n 1 i∈G 2 w i ). A detailed derivation of (10) from (9) is provided in the Appendix. We solve the minimization problem (10) by the Lasso technique.
Consequently, solving the complicated optimization problem (7) is decomposed to the simple tasks of iteratively solving a Glasso estimate for C and a Lasso estimate for δ 2 until both of them are converged. We summarize the above estimation procedure for the proposed model in Algorithm 1.
Step 0: Set an initial value of δ 2 .
Note that there are two tuning parameters λ 1 and λ 2 in the optimization problem (7).
To choose their optimal values, we minimize a BIC-type criterion proposed by Wang et al. (2007) as where v(δ 2 ) and v(Ĉ) stand for the number of nonzero entries in the estimatesδ 2 andĈ, respectively. This criterion enjoys consistency properties and has been commonly used in literature (Zou and Zhang, 2009;Lv and Fan, 2009;Armagan et al., 2013).

Model Prediction
In this section, we demonstrate how to conduct model prediction by the proposed method.
For convenience, let us write where µ 1X and µ 2X are p−1 dimensional vectors representing the means of variable X for two classes, and Σ X is the (p−1)×(p−1) covariance matrix of X. The estimatesμ 1 ,μ 2 ,Ĉ,Σ can be partitioned accordingly. From model assumption (1) as well as the property of multivariate Therefore, the prediction for the quantitative variable y from a new observation x iŝ Note thatΣ ′ Xy whereĉ 2 y is a scalar, implying that the sparsity ofĈ Xy will lead to the sparse model for the prediction of y.
On the other hand, the prediction for the qualitative variable Z by the proposed model is naturally based on the estimated LDA classification rule of (3) as From Equations (11) and (12), however, we note that the prediction of one response variable depends on the information of the other. To address this issue, we propose to calculate two candidate values of y for a new observation x by Equation (11) for two different classes, denoted byŷ 1 andŷ 2 . Then the conditional probability densities p(W = (x ′ ,ŷ 1 ) ′ |G 1 ) and p(W = (x ′ ,ŷ 2 ) ′ |G 2 ) can be estimated via the density functions of N(μ 1 ,Σ) and N(μ 2 ,Σ).
Denote such two values asp 1 andp 2 . The prediction of y at this new observation is then obtained asŷ k corresponding to the larger value ofπ kpk , k = 1, 2. To express it clearly, we describe the above steps of the model prediction in Algorithm 2 for a new observation x.
It is seen that in Algorithm 2, we obtain the prediction of y first, and then predict Z with the value ofŷ. One would argue that it is not a unique way of making predictions on QQ responses, as we may also predict Z first and then variable y. The following proposition provides an interesting insight into this issue. Proposition 1 implies that we can predict the response variable Z by simply comparing values ofπ kpk instead of employing LDA. Therefore, the order of which response variable to be predicted first is not a concern. Actually, the Step 2a and Step 2b are equivalent to the following Step 2 as Step 2: Ifπ 1p1 >π 2p2 , letŷ =ŷ 1 andẐ = 1; otherwise letŷ =ŷ 2 andẐ = 2.

Extension to Multi-Class Qualitative Response
The proposed generative modeling approach also has the advantage to enable the GAQQ to deal with the qualitative response with multiple classes, i.e., the qualitative variable Z ∈ {1, 2, . . . , K}. In such cases, the GAQQ method is extended and expressed as G k : Based on a baseline class G 1 , we regularize on the difference between means through µ k − µ 1 for k = 2, 3, . . . , K. The objective function is thus formulated as The subsequent derivation follows similar steps as that described in Section 2.1. With a little abuse of notation, let Kδ k = µ k − µ 1 for k = 1, 2, . . . , K and Kγ = K k=1 µ k , then we have µ k = γ − K g=2 δ g + Kδ k , k = 1, 2, . . . , K. As a result, the optimization problem (13) can be re-written as min (δ 2 ,...,δ K ,γ,C) Let n k represent the number of observations belonging to class G k . The maximum likelihood estimator of γ from (14) isγ =w + K g=2 δ g − K n K g=2 n g δ g . Consequently, the optimization problem (14) becomes (15) can be decomposed as one Glasso problem min C −n ln |C| + tr(CS) + λ 1 ||C|| 1 , and min (δ 2 ,...,δ K ) The optimization problem (16) is equivalent to the following K − 1 Lasso regressions sepa- The detailed derivation from (16) to (17) is provided in the Appendix. Therefore, the parameters δ k and C can be solved iteratively until convergence following the spirit of Algorithm 1. The optimal values of tuning parameters are chosen by the BIC-type criterion extended for the multi-class problem as where v(δ) represents the number of nonzero entries in all the estimatesδ k .
For a new observation x, the quantitative response y is predicted, similarly as in Algo- The class label is estimated asẐ = arg max kπ kpk , or equivalently by the LDA rule aŝ

Theoretical Properties
In this section, we will investigate the asymptotic optimality of the classification rule by the proposed GAQQ method in Theorem 1 to Theorem 3. The asymptotic consistency properties of the prediction of y by the GAQQ method are established in Theorem 4. For the proposed classification rule, we first establish the theoretical results for the multi-class problem and then provide a thorough discussion of the two-class case. We use the same definition of asymptotic optimality for a classification rule as defined in Shao et al. (2011).
Denote by R Bayes and R P ROP (T ) the Bayes error and the conditional misclassification rate of the proposed rule, where T denotes the training samples. The asymptotic optimality for a classification rule is defined as follows.

Definition 1. Let T be a classification rule with conditional misclassification rate R T (T ),
given the training samples T .
(1) T is asymptotically optimal if R T (T )/R Bayes P → 1. ( Note that if lim n→∞ R Bayes > 0, then the asymptotically sub-optimality is the same as the asymptotically optimality. To facilitate the construction of theoretical results, we need to introduce some notation and make assumptions on the true model. Define the true values of µ k , Σ, C and δ k as µ 0 k , Σ 0 , C 0 and δ 0 kX is a p − 1 dimensional vector representing the true mean difference of variable X between classes G 1 and G k . Denote the true inverse covariance matrix of variable X by C 0 , which is the set containing location indices of the nonzero entries in δ 0 k . Lets k be the cardinalities of set S δ k . Define s k =s k if δ 0 ky = 0; otherwise s k =s k − 1. That is, s k is the number of nonzero entries of δ 0 kX . Additionally, we use the same sparsity measure on Σ 0 = (σ 0 ij ) 1≤i,j≤p as in Bickel and Levina (2008), which is S h;p = max i≤p p j=1 |σ 0 ij | h where 0 ≤ h < 1 and 0 0 is defined to be 0. Hence firstly, S 0;p equals the maximum of the numbers of nonzero entries in each row of the matrix Σ 0 . In this case, a smaller value of S 0;p compared with p implies a sparse structure in matrix Σ 0 . Secondly, if S h;p is smaller than p for 0 < h < 1, it indicates that many entries of matrix Σ 0 are very small. Moreover, we assume the following regularity conditions.
where λ min (C 0 ) and λ max (C 0 ) are the minimum and maximum eigenvalues of matrix , and J c represents the complement set of J.
• (C4) Irrepresentable condition: without loss of generality, write δ 0 where 1 is a p −s k dimensional unit vector, and the inequality holds element-wise.
By conditions (C1) and (C2), Rothman et al. (2008) and Lam and Fan (2009) derived the convergence rate of Glasso estimate. We thus have where d n = S h;p ( log p n ) (1−h)/2 , and A is the matrix spectral norm defined as the squared root of the maximum eigenvalue of matrix A ′ A. The conditions (C2) and (C3) are used in Bühlmann and Van De Geer (2011) to study the theoretical property of Lasso estimate, and we have where b Under conditions (C4) and (C5), Zhao and Yu (2006) showed that the Lasso estimate is model selection consistency, which will be used for investigating δ k in (17). Condition (C6) requires that all the classes should be separated from each other. Also note that condition (C6) is equivalent to that ∆ k is bounded away from 0. The condition (C7) guarantees a balanced sample size for each class, which is commonly used in literature to bound the term log π k π 1 in the LDA rule for establishing the properties of classification rules. Based on the above results, we present the following theories on the consistency of the classification rule by the proposed method.

Theorem 1. Assume that conditions (C1) -(C7) hold, and
Then the proposed rule for the multi-class problem is asymptotically sub-optimal if either one of the following two conditions is satisfied Theorem 1 establishes the sub-optimality property of the proposed classification rule for the multi-class problem. In the case of two-class problem, the Bayes error can be expressed in a closed form of R Bayes = Φ(−∆ 2 /2) when the data are from normal distribution, where Φ represents the cumulative distribution function of N(0, 1), and ). Accordingly, in Theorems 2 and 3, we can compute the convergence rate of the proposed rule for the two-class problem, and subsequently investigate its properties.
Theorem 3. Assuming all the conditions in Theorem 2 are satisfied, we have (1) if ∆ 2 is bounded, then the proposed rule is asymptotically optimal and R P ROP (T ) (2) if ∆ 2 → ∞, then the proposed rule is asymptotically sub-optimal; (3) if ∆ 2 → ∞ and ξ n ∆ 2 2 → 0, then the proposed rule is asymptotically optimal.
Theorem 2 provides the convergence rate of the proposed classification rule for the twoclass problem with respect to ξ n . Base on such a result, Theorem 3 demonstrates that the property of the proposed classification rule (optimality or sub-optimality) depends on the scenarios of the true model's ∆ 2 . Specifically, (1) when ∆ 2 is bounded, i.e. lim n→∞ R Bayes > 0, then R P ROP (T ) converges in probability to the same limit as R Bayes .
R Bayes → 0, then R P ROP (T ) P → 0; in this case, if we further have ξ n ∆ 2 2 → 0, then R P ROP (T ) and R Bayes have the same convergence rate.
Next, we derive the consistency property for the proposed estimate of y. Denote byŷ P the predicted value of y obtained from the proposed model. Defineŷ B to be the predicted value of y for x when all parameters are known. Specifically, first obtain the class label k via the Bayes Now we establish the theoretical results ofŷ P in Theorem 4. This result compares the MSE of the proposed estimate of y with that from the optimal Bayes rule (under all parameters known). Since the classification errors from a classification rule might be larger than 0, the MSE ofŷ may not converge to 0 even though the sample size n is sufficiently large. Here we adopt the MSE Bayes as a reasonable performance benchmark to evaluate the property of the proposed model with respect to y. Theorem 4 states that the difference of MSE between the proposed model and the Bayes method converges to 0 in probability.

Two-class Settings of the Qualitative Response
In this section, we evaluate the performance of the proposed GAQQ method for a binary response Z under different inverse covariance matrices C and mean differences δ 2 . The proposed GAQQ model is compared with several benchmark methods, denoted as GLDA, CL, and ENET, which use the predictor variables X to predict Z and y. The GLDA employs the LDA classification rule for Z using the generalized inverse of the sample covariance matrix of X when p > n. The CL method applies the LPD technique introduced by Cai and Liu (2012) to predict the response variable Z based on X. With their estimated class label of Z, the GLDA and CL predict y by Equation (11). The ENET method uses the elastic-net logistic model (Zou and Hastie, 2005) on predictor variables X to fit the qualitative response Z and hence predicts Z for the testing data. For the quantitative response y, the ENET separately fits two elastic-net linear regressions for two classes using training data and then predicts y in the testing data based on its estimated Z. The tuning parameters of the CL and ENET methods are chosen by cross-validation.
Regarding the inverse covariance matrix C, we consider the following five structures in the simulation, which are commonly used in the literature (Yuan and Lin, 2007;Kang and Deng, 2020).
• Model 3. C 3 is generated by randomly permuting rows and corresponding columns of the matrix C 2 . This matrix has more sparsity as the dimensionality increases. Models 2 and 3 are relatively dense matrices, and they also become more sparse when the dimensionality increases. All of these four matrices have sparse structures to some extent, while Model 5 is a general sparse matrix with no structure, which is similarly used in Bien and Tibshirani (2011).
For the mean difference δ 2 , we consider two different levels of sparsity. The µ 1 is the vector with all elements zeros. Then generate µ 2 such that (S1): 25% of the elements in The same procedure is employed to generate the testing data, which is used to evaluate the prediction performance of y and Z for different compared methods. We consider the root mean squared prediction error RMSPE = 1 n n i=1 (y i −ŷ i ) 2 to measure the prediction accuracy for the quantitative response y, whereŷ i represents the predicted value. The prediction performance of the qualitative response Z is measured by the misclassification error ME = 1 n n i=1 I(z i =ẑ i ), whereẑ i is the predicted value of z i and I(·) is an indicator function.
Tables 1 and 2 report the averaged MEs in percentage and averaged RMPSE, as well as their corresponding standard errors in parenthesis for each approach over 100 replications. It can be seen from Table 1 that the proposed method generally outperforms other approaches with respect to MEs. Such an advantage becomes more significant as the underlying models are more sparse. Specifically, in the scenario of S1 = 25% and p = 40, the proposed GAQQ method does not perform as well as others, since the underlying models in this scenario are the least sparse, especially for the dense models 2 and 3. In contrast, the proposed method produces relatively better comparison results in the scenario of S2 = 75% and p = 40, where the true mean difference is more sparse. Furthermore, this advantage of the proposed method is well evidenced in the scenario of p = 80, and even more notable when p = 200 with its substantially lower MEs than other methods.
From Table 2, we observe that the proposed method generally gives superior performance over other compared approaches for each scenario in predicting the quantitative response y.
The possible explanations are in two folds. First, the proposed GAQQ method provides an accurate classification of the qualitative response Z. Second, the proposed GAQQ has a proper estimation of C by the regularization that is used in the prediction of quantitative response y according to (11), resulting in an improvement of the prediction accuracy. It is also seen that the CL and GLDA methods are comparable in some cases, possibly because both of them use the generalized inverse of the sample covariance of X forΣ −1 X in the prediction of quantitative response y in (11). But the CL method is generally better since it has more accurate classification results than the GLDA in Table 1.

Multi-class Settings of the Qualitative Response
Now, we examine the performance of the proposed GAQQ method for multi-class settings of the qualitative response. We consider p = 200 and K = 4 classes of qualitative response Z with training sizes n 1 = n 2 = n 3 = n 4 = 30 for Models 1 -5 of inverse covariance matrix C. Let µ kj represent the jth entry of the mean value µ k . Generate µ kj = 0.5 * k + u kj for j = 2k − 1, 2k, 2k + 1, . . . , 2k + 6, otherwise µ kj = 0, where u kj is from Unif (−1, 1). The training data are generated from N(µ k , C −1 ), and the testing data follow the same generation procedure. We compare the proposed method with the GLDA, as well as the estimators where the latter two methods are designed for multi-class problems. We use the WT and CHWE models to first predict the class label Z for the testing data, and then the response y is estimated, by the multivariate normal property, asμ ky +Σ ′ The results of performance measures, ME and RMSPE are summarized in Table 3 based on 100 replications. One can see that the GAQQ method performs better than the GLDA as well as the WT method, and is comparable with the CHWE in terms of the MEs. Besides, the GAQQ method gives the best performance among the compared approaches with significantly lower values of RMSPE.

Case Studies
In this section, we apply the proposed GAQQ method to two real-data case studies. The first one is from the study of Heusler compounds in material science and the second one is from the study of molecular diagnostics of Ulcerative colitis and Crohn's disease. Although from different fields, both problems contain QQ responses with high-dimensional predictors, and the proposed GAQQ method appears to have much better performance in terms of prediction accuracy compared with other methods.
The case study on material sciences is regarding the Heusler compounds, which are a large family of intermetallics with more than 1000 known members. Many Heusler compounds have shown exotic properties, such as superconductivity and topological band structures, which have promising applications for quantum computing. Understanding the thermody- namic stability of Heusler compounds lays the foundation for exploiting the large chemical space to discover and design new functional Heusler materials (Liu et al., 2016). To determine the thermodynamic stability of Heusler compounds, there are two key metrics: the mixing enthalpy (quantitative response) and the global stability based on hull energy (binary qualitative response). The comprehensive database of 180628 full Heusler structures was built by collecting the relevant structural and energetic data from the Materials Project (Jain et al., 2013), OQMD (Saal et al., 2013), and AFLOW (Curtarolo et al., 2012). These data were calculated using first-principles methods based on density functional theory, and it was extremely computationally expensive (taking hours) to generate one entry of the data.
Therefore, a statistical model that can accurately predict the thermodynamic stability for any elemental and compound features is a useful surrogate of the first-principle computation models.
Since there is an intrinsic relationship between two QQ responses, the proposed GAQQ method is suitable to improve the prediction accuracy by jointly fitting them together. To demonstrate the GAQQ method in the scenario when the number of predictors is large relative to the size of the data, we randomly choose 150 samples from each class of the binary response. We delete the predictor variables whose standard deviations are less than 1.0e −6 , resulting in 157 predictors of elemental and compound features. To examine the prediction performance of the GAQQ method and other comparison methods, we randomly divide data into a training set with a size of 200 and a testing set with a size of 100. Table   4 reports the prediction performance results based on 50 random splits of the Heusler data.
From the results, it is seen that the proposed GAQQ performs much better than other methods in comparison, with the smallest values for the misclassification error (ME) and the root mean squared prediction error (RMSPE). The second data for the case study considers the multi-class settings of the qualitative response. The IBD gene data (Burczynski et al., 2006)

Discussion
In this work, we propose a generative modeling approach to jointly model the data with QQ responses, which is a new perspective different from existing methods in the literature.
By fully exploring the joint distribution of the QQ responses and predictor variables, the proposed method enables efficient parameter estimation, accurate prediction, and lays a good foundation for investigating the asymptotic properties. The proposed model can be naturally extended to the situation for multiple quantitative responses.
One further research direction is to accommodate a more flexible structure on the joint distribution of QQ responses and predictor variables. For example, one can extend the LDA for the classification of the qualitative response to the quadratic discriminant analysis (QDA). The QDA is more flexible with different covariance structures in each class, but its estimation for high-dimensional data would encounter more difficulty due to a large number of parameters. Besides, the derivation of its asymptotic properties is much more technically complicated (Li and Shao, 2015). Another research direction is to apply the generative modeling approach for the data with semi-continuous responses (Wang et al., 2020), or the ordinal and quantitative responses. One may employ the ordinal regression for the ordinal response, and then derive its joint likelihood function with appropriate regularization.

Appendix
Derivation from (9) to (10). Let C denote a generic constant thereafter.
Derivation from (16) to (17). For δ j , j = 2, 3, . . . , K, where Hence, formula (20) is equal to , where x and y are multivariate variables. For a given value of Proof. We need to search for y to minimize where C is a constant not depending on y. Taking derivative of L(y) and setting to zero yields ∂L(y) ∂y This, together with a property of block matrix that Ω ′ XY = −Ω Y Σ ′ XY Σ −1 X , completes the proof.
The inequality p(x ∈ G 1 |X = x) > p(x ∈ G 2 |X = x) in Proposition 2 indicates that the LDA rule assigns x to G 1 . Therefore, Proposition 2 implies that Step 2b of Algorithm 2 is equivalent to applying the LDA classification rule directly on x instead of w = (x ′ ,ŷ) ′ .
This fact enables us to give theoretical proof for the consistency properties of the proposed classification rule based on variable X rather than W = (X ′ , y) ′ . Before the proof of Theorem 1, we present Lemmas 2 -4.
Write µ 0 k = ((µ 0 kX ) ′ , µ 0 ky ) ′ , where µ 0 kX is the true mean value of variable X for class G k . Correspondingly, writeμ k = (μ ′ kX ,μ ky ) ′ . Let a n ≍ b n represent two sequences a n and b n to be the same order. Now we state Lemma 3.
Proof. It is not difficult to derive Consequently, applying the Cauchy-Schwarz inequality together with Equations (25) and (26), we can obtain Similarly, As a result, according to Lemma 3, we have From Lemma 4, note that when t = k, we have (µ 0 With Lemmas 2 -4, we are ready to complete the proof of Theorem 1.
LetẐ P ROP andẐ Bayes denote the predicted class labels obtained by the proposed model and the Bayes rule, respectively. For simplicity, we assume π 1 = π 2 = . . . = π K instead of condition (C7) in the proofs of Theorems 1 -3 with no influence on the theoretical results, since condition (C7) is only used to bound the term log π k π 1 of the LDA rule. Define XδkX for a new sample x. Then for any ǫ > 0, ≤ Pr(|θ k − ϑ k | ≥ ǫ 2 for some k) + Pr(|ϑ k − ϑ l | ≤ ǫ for some k, l).
where the last inequality is obtained by condition (C6). Secondly, we bound the term Pr(|θ k −ϑ k | ≥ ǫ 2 for some k).
The conditional misclassification rate is By the result (18), we havê From the result (19), together with E[(δ 0 Hence we havê Therefore, by Lemma 3, we obtain Similarly, we haveδ which proves theory. To establish the theoretical results in Theorem 3, we need a lemma from Shao et al. (2011). We state it here for completeness, and then prove Theorem 3.

Lemma 5. Let a
(1) n and a (2) n be two sequences of positive numbers such that a (1) n → ∞ and a Proof. See the proof of Lemma 1 in Shao et al. (2011).