# NetDA: An R Package for Network-Based Discriminant Analysis Subject to Multilabel Classes

Authors: Li-Pang Chen
Journal: Journal of Probability and Statistics
 (2022)

---

## Abstract

In this paper, we introduce the R package NetDA, which aims to deal with multiclassification with network structures in predictors accommodated. To address the natural feature of network structures, we apply Gaussian graphical models to characterize dependence structures of the predictors and directly estimate the precision matrix. After that, the estimated precision matrix is employed to linear discriminant functions and quadratic discriminant functions. The R package NetDA is now available on CRAN, and the demonstration of functions is summarized as a vignette in the online documentation.


---

## Body

## 1. Introduction

Multiclassification, known as a classification problem that the number of classes is greater than two, is a great challenge in data science. In supervised learning, discriminant analysis has been a useful method to do classification. In the conventional method (e.g., Hastie et al. [1]; Section 4.3), linear discriminant functions, which are formulated in terms of mean vectors and the inverse of covariance matrices of the predictors, are used to classify subjects. In addition, some advanced methods have been proposed to address complex settings in the past literature. For example, Guo et al. [2] discussed the LDA method and its application in microarray data analysis. Safo and Ahn [3] studied generalized sparse linear discriminant analysis for multilabel responses. In the presence of high-dimensional predictors, several advanced approaches have also been explored (e.g., Clemmensen et al. [4]; Witten and Tibshirani [5).However, network structures of predictors, which reflect (pairwise) dependence among predictors, are ubiquitous in data analysis [6]. In the recent developments, Chen et al. [7] proposed a graphical-based logistic regression model. He et al. [8] proposed surrogate variables that were transformed from network structures and implemented them to the support vector machine. Regarding the framework of discriminant analysis, Cai et al. [9] and Liu et al. [10] developed graph-based linear discriminant analysis, but their approaches are restricted to binary responses. Moreover, for general data analysts, it is important for them to directly implement existing software and do data analysis. However, rare software related to classification with network structure accommodated has been available. While some R packages related to discriminant analysis exist, such as MASS, sparseLDA, and penalizedLDA, they are not able to handle network structures in predictors.Motivated by these concerns and to address these challenges, we follow the strategy proposed by Chen [11] and develop an R package, which is called NetDA. Under the normality assumption for predictors, we apply the graphical lasso method to estimate precision matrices and the corresponding network structures for the predictors. Since precision matrices are the inverse of covariance matrices, it motivates us to directly implement them to linear/quadratic discriminant functions. This strategy is different from the conventional linear discriminant analysis that simply employs empirical estimates of covariance matrices. Moreover, the other issue is prediction. Based on fitted models and predicted values, we also develop a function that contains several commonly used criteria to assess the performance of classification and prediction.The article is organized as follows. Section 2 introduces the data structure and outlines the methodology in our package. Section 3 describes the usage of the package NetDA. Section 4 illustrates the package by a real dataset. We finally conclude the article in Section 5.


## 2. Overview of Methodology

In this section, we primarily overview the data structure and network-based discriminant analysis proposed by Chen [11].


### 2.1. Data Structure

Suppose that the data containn subjects that come from I classes, where I is a fixed integer greater than 2 and the classes are nominal. Let ni be the size in class iwith i=1,…,I, and hence, n=∑i=1Ini. Let Y denote the n-dimensional vector of responses with the jth component being Yj=i, which reflects the class membership that the jth subject is in the ith class for i=1,…,I and j=1,…,n.Letp>1 denote the dimension of predictors for each subject. Define X=Xj,k as the n×p matrix of predictors for j=1,…,n and k=1,…,p, where the component Xj,k represents the kth predictor for the jth subject. Furthermore, let X⋅k=X1,k,…,Xn,k⊤ represent the n-dimensional vector of the kth predictor in the kth column of X, and let Xj⋅=Xj,1,…,Xj,p⊤ denote the p-dimensional predictor vector for the jth subject in the jth row of X. Let Xj⋅,Yj:j=1,…,n denote an independent and identically distributed (i.i.d.) sample. We let lower case letters represent realized values for the corresponding random variables. For example, xj⋅ stands for a realized value of Xj⋅.


### 2.2. Gaussian Graphical Models

Fori=1,…,I and j=1,…,n, let fj|ixj⋅ denote the conditional probability density function of the predictor Xj⋅ taking a value xj⋅ given that subject j comes from the i th class. We particularly consider the case where the conditional distribution Xj⋅ given Yj=i is assumed to be a multivariate normal distribution with a mean vector μi and a positive-definite covariance matrix Σi. Then, the conditional probability density function is given by(1)fj|ixj⋅=12πp/2Σi1/2exp−12xj⋅−μi⊤Σi−1xj⋅−μi.Moreover, by the suitable reparametrization (e.g., Hastie et al. [12]; p. 246), we can transfer (1) to the Gaussian graphical model (GGM) based on class i. The exact formulation is given by(2)fxj⋅;βi,Θi=exp∑r∈Vβirxj,r−12∑s,t∈Eθi,stxj,sxj,t−12logdetΘi2π,where V=1,…,p includes all the indices and E⊂V×V contains all pairs with unequal coordinates, which yields a graph G≜V,E (e.g., Hastie et al. [1]; Chapter 17); Θi=θi,st is the p×p symmetric precision matrix with Θi≜Σi−1; and βi=βi1,…,βip⊤ is the p-dimensional vector of parameters with μi≜−Θi−1βi. Our main interest is to estimate Θi, since, as we will see later, the main concern in discriminant analysis is to estimate the inverse of the covariance matrix. On the other hand, from the perspective of graphical models, nonzero θi,st implies that Xj,s and Xj,t are conditionally dependent given other variables in class i, while zero value of θi,st gives conditional independence of Xj,s and Xj,t given other variables. Thus, the precision matrix Θi reflects the network structure of the predictors.In the past literature, graphical LASSO (GLASSO) [13] is a common method to estimate Θi. The key idea of GLASSO is based on the likelihood function. To see this, we follow the similar discussion in page 247 of Hastie et al. [12] and write the log-likelihood function of Θi based on (2) with βi=0:(3)ℒΘi=1ni∑j=1nilogfxj⋅;Θi=logdetΘi−traceSiΘi,where Si=1/ni∑j=1nixj⋅xj⋅⊤, and trace⋅ is the sum of diagonal entries for a square matrix, and(4)logdetΘi=∑k=1plogλΘi,ifΘiis positive definite,−∞,otherwise,with λΘi being the jth eigenvalue of Θi. Assume that the precision matrix Θi is sparse. To estimate θi,st and identify network structures by retaining dependent pairs of vertices and removing independent ones, we apply the L1-norm as a constraint to achieve the desired result. In other words, the estimator of Θi can be obtained by the following optimization:(5)Θ^i=argmaxΘilogdetΘi−traceSiΘi−ζρΘi,where ρΘi≜∑s≠tθi,st is the penalty function and ζ is a tuning parameter. The optimization problem (5) is called GLASSO [13]. The detailed algorithm can be found in page 248 of Hastie et al. [12], and the estimator of ζ can be determined by Bayesian information criterion (BIC). By the similar discussion in Yuan and Lin [14], the estimated network structure determined by (5) is equal to the true graph with probability approaching one under suitable conditions. For the computation, the R package glasso can be implemented to derive the estimate Θ^i.


### 2.3. Discriminant Analysis

The idea of discriminant analysis is to model the distribution of the predictorsXj⋅ separately for each of the response classes Yj, and then to use the Bayes theorem to describe the conditional probabilities PYj=i|Xj⋅=xj⋅ (e.g., James et al. [15]).Specifically, letπi≜PYj=i denote the probability that the j th subject is randomly selected from class i so that ∑i=1Iπi=1. Moreover, applying the Bayes theorem to the conditional density function fj|ixj⋅ and π gives the posterior probability as(6)PYj=i|Xj⋅=xj⋅=fj|ixj⋅πi∑l=1Ifj|lxj⋅πl,for i=1,…,I and j=1,…,n.To compare two classesi and l with i≠l, we calculate the log-ratio of (6), given by(7)logPYj=i|Xj⋅=xj⋅PYj=l|Xj⋅=xj⋅=logfj|ixj⋅fj|lxj⋅+logπiπl.With a distribution assumption (1) imposed, we now have two scenarios for the first term in the right-hand side of (7), say logfj|ixj⋅/fj|lxj⋅.Scenario 1.
If the covariance matricesΣi in (1) are assumed to be common, that is, Σi=Σ for every i with Σ being a positive definite matrix, (7) becomes(8)logπiπl−12μi+μl⊤Σ−1μi+μl+xj⋅⊤Σ−1μi+μl.
If equation (8) >0, then PYj=i|Xj⋅=xj⋅>PYj=l|Xj⋅=xj⋅, showing that subject j with predictors Xj⋅=xj⋅ is more likely selected from class i than from class l. Consequently, (8) defines a boundary between classes i and l in the sense that there is a linear function in xj⋅ separating classes i and l.
Motivated by the form of (8), we consider a linear function in x as(9)δix≜logπi−12μi⊤Σ−1μi+x⊤Σ−1μi.
Moreover,μi and πi can be empirically estimated, respectively, as(10)μ^i=1ni∑Yj=ixj⋅,π^i=nin.
For the estimation ofΣ−1, or equivalently Θ, we adopt (5) by pooling all subjects in the dataset and denote Θ˜ as the estimator of Θ. Therefore, (9) can be estimated as(11)δ˜ix=logπ^i−12μ^i⊤Θ˜μ^i+x⊤Θ˜μ^i,and we call (11) the network-based linear discriminant function (NetLDA) and it is used to determine the class label for a new observation. For the prediction of a new subject with the predictor x˜, we first calculate δ˜ix˜ using (11) for i=1,…,I. Next, we find i∗ that is defined as(12)i∗=argmaxi=1,…,Iδ˜ix˜;and the class label for this subject is then predicted as i∗.Scenario 2.
We allowΣi≠Σl, or equivalently, Θi≠Θl, for any i≠l and i,l=1,…,I. Then, under a distribution assumption (1), we have(13)logfj|ixj⋅fj|lxj⋅=log1/2πp/2Σi1/2exp−1/2xj⋅−μi⊤Θixj⋅−μi1/2πp/2Σl1/2exp−1/2xj⋅−μl⊤Θlxj⋅−μl=logΘl−1/2Θi−1/2+12xj⋅−μl⊤Θlxj⋅−μl−xj⋅−μi⊤Θixj⋅−μi.
Replacing the first term in the right-hand side of (7) by (13) yields(14)logπiπl+logΘl−1/2Θi−1/2+12xj⋅−μl⊤Θlxj⋅−μl−xj⋅−μi⊤Θixj⋅−μi.
Therefore, based on (14), we further define a quadratic function of x based on the class i:(15)φix≜logπi+12logΘi−12x−μi⊤Θix−μi.
Fori=1,…,I, the estimator of the precision matrix Θi, denoted as Θ^i, is obtained by (5) based on the predictor information in class i; μi and πi can be estimated by (10). Therefore, (15) can be estimated by(16)φ^ix=logπ^i+12logΘ^i−12x−μ^i⊤Θ^ix−μ^i.
The function (16) is called the network-based quadratic discriminant function (NetQDA) and is used to determine the class label for a new observation. For the prediction of a new subject with the predictor x˜, we first calculate φ^ix˜ using (16) for i=1,…,I. Next, we find i∗ that is defined as(17)i∗=argmaxi=1,…,Iφ^ix˜;and the class label for this subject is then predicted as i∗.


### 2.4. Assessment of Classification and Prediction

We first introduce micro-averaged metrics (e.g., Chen et al. [7]). Let V and T represent the classes of the subject indexes for validation and training datasets, respectively. Let n=T and m=V denote the sizes of the training and validation data, respectively. We use the training data T to fit models and then apply fitted models to compute predicted values Y^j for j∈V. After that, for i=1,…,I, we calculate the number of the true positives (TP), the number of the false positives (FP), and the number of the false negatives (FN) under the validation data V, respectively:(18)TPi=∑j∈VIYj=i,Y^j=i,FPi=∑j∈VIYj≠i,Y^j=i,FNi=∑j∈VIYj=i,Y^j≠i,where I⋅ is the indicator function.We define precision (PRE) and recall (REC) under the validation dataV, respectively, as(19)PRE=∑i=1ITPi∑i=1ITPi+∑i=1IFPi,REC=∑i=1ITPi∑i=1ITPi+∑i=1IFNi.Then, micro-F-score is defined as(20)F=2×PRE×RECPRE+REC.According to definitions in (19), when all subjects are correctly classified, FP and FN are equal to zero, yielding that PRE and REC are equal to one; if all subjects are falsely classified, then TP is equal to zero, and thus, PRE and REC are equal to zero. Therefore, values of PRE and REC are between zero and one. Moreover, under the range 0,1, the F-score falls in 0,1 as well by treating 0/0 as zero. In principle, the higher values of PRE, REC, and F-score reflect the better performance and the more accurate classification (e.g., Chen et al. [7]).In addition to criteria above, the other commonly used criterion is the adjusted Rand index (ARI). Fori,i′=1,…,I and under the validation data V, we define(21)nii′=∑j∈VIYj=i,Y^j=i′.Moreover, we defineai=∑i′=1Inii′ for i=1,…,I and bi′=∑i=1Inii′ for i′=1,…,I. Then, ARI under the validation data V is defined as Hubert and Arabie [16].(22)ARI=∑i,i′nii′2−∑iai2∑i′bi′2/n2∑iai2+∑i′bi′2/2−∑iai2∑i′bi′2/n2.As mentioned in Hubert and Arabie [16], ARI is bounded above by one, and the higher value of ARI indicates the more accurate classification.


### 2.5. Benchmark of NetDA

The conventional linear discriminant methods (e.g., MASS, sparseLDA, and penalizedLDA) aim to adopt (9) with Θ estimated by the inverse of the empirical estimator of Σ. However, this approach may encounter cumbersome computation or possible singularity when calculating inverse matrices. In addition, if Θ is sparse, all entries in empirical estimators of Θ are nonzero, and it indicates that some unconnected pairs of predictors may be falsely included. As a result, imprecise estimator of Θ may implicitly affect the performance of classification.Unlike existing methods, the first contribution of NetDA is to estimateΘ directly. The graphical lasso method is a tool to identify zero entries in Θ and estimate nonzero ones. This approach enables us to retain connected pairs of predictors and exclude unconnected ones. Moreover, our implementation can avoid computing inverse of Σ.The second contribution of NetDA is to handle heterogeneous network structure that stratifies the predictor information by class when characterizing the predictor network structures, and is able to deal with multiclassification by adopting network structure in each class.


## 2.1. Data Structure

Suppose that the data containn subjects that come from I classes, where I is a fixed integer greater than 2 and the classes are nominal. Let ni be the size in class iwith i=1,…,I, and hence, n=∑i=1Ini. Let Y denote the n-dimensional vector of responses with the jth component being Yj=i, which reflects the class membership that the jth subject is in the ith class for i=1,…,I and j=1,…,n.Letp>1 denote the dimension of predictors for each subject. Define X=Xj,k as the n×p matrix of predictors for j=1,…,n and k=1,…,p, where the component Xj,k represents the kth predictor for the jth subject. Furthermore, let X⋅k=X1,k,…,Xn,k⊤ represent the n-dimensional vector of the kth predictor in the kth column of X, and let Xj⋅=Xj,1,…,Xj,p⊤ denote the p-dimensional predictor vector for the jth subject in the jth row of X. Let Xj⋅,Yj:j=1,…,n denote an independent and identically distributed (i.i.d.) sample. We let lower case letters represent realized values for the corresponding random variables. For example, xj⋅ stands for a realized value of Xj⋅.


## 2.2. Gaussian Graphical Models

Fori=1,…,I and j=1,…,n, let fj|ixj⋅ denote the conditional probability density function of the predictor Xj⋅ taking a value xj⋅ given that subject j comes from the i th class. We particularly consider the case where the conditional distribution Xj⋅ given Yj=i is assumed to be a multivariate normal distribution with a mean vector μi and a positive-definite covariance matrix Σi. Then, the conditional probability density function is given by(1)fj|ixj⋅=12πp/2Σi1/2exp−12xj⋅−μi⊤Σi−1xj⋅−μi.Moreover, by the suitable reparametrization (e.g., Hastie et al. [12]; p. 246), we can transfer (1) to the Gaussian graphical model (GGM) based on class i. The exact formulation is given by(2)fxj⋅;βi,Θi=exp∑r∈Vβirxj,r−12∑s,t∈Eθi,stxj,sxj,t−12logdetΘi2π,where V=1,…,p includes all the indices and E⊂V×V contains all pairs with unequal coordinates, which yields a graph G≜V,E (e.g., Hastie et al. [1]; Chapter 17); Θi=θi,st is the p×p symmetric precision matrix with Θi≜Σi−1; and βi=βi1,…,βip⊤ is the p-dimensional vector of parameters with μi≜−Θi−1βi. Our main interest is to estimate Θi, since, as we will see later, the main concern in discriminant analysis is to estimate the inverse of the covariance matrix. On the other hand, from the perspective of graphical models, nonzero θi,st implies that Xj,s and Xj,t are conditionally dependent given other variables in class i, while zero value of θi,st gives conditional independence of Xj,s and Xj,t given other variables. Thus, the precision matrix Θi reflects the network structure of the predictors.In the past literature, graphical LASSO (GLASSO) [13] is a common method to estimate Θi. The key idea of GLASSO is based on the likelihood function. To see this, we follow the similar discussion in page 247 of Hastie et al. [12] and write the log-likelihood function of Θi based on (2) with βi=0:(3)ℒΘi=1ni∑j=1nilogfxj⋅;Θi=logdetΘi−traceSiΘi,where Si=1/ni∑j=1nixj⋅xj⋅⊤, and trace⋅ is the sum of diagonal entries for a square matrix, and(4)logdetΘi=∑k=1plogλΘi,ifΘiis positive definite,−∞,otherwise,with λΘi being the jth eigenvalue of Θi. Assume that the precision matrix Θi is sparse. To estimate θi,st and identify network structures by retaining dependent pairs of vertices and removing independent ones, we apply the L1-norm as a constraint to achieve the desired result. In other words, the estimator of Θi can be obtained by the following optimization:(5)Θ^i=argmaxΘilogdetΘi−traceSiΘi−ζρΘi,where ρΘi≜∑s≠tθi,st is the penalty function and ζ is a tuning parameter. The optimization problem (5) is called GLASSO [13]. The detailed algorithm can be found in page 248 of Hastie et al. [12], and the estimator of ζ can be determined by Bayesian information criterion (BIC). By the similar discussion in Yuan and Lin [14], the estimated network structure determined by (5) is equal to the true graph with probability approaching one under suitable conditions. For the computation, the R package glasso can be implemented to derive the estimate Θ^i.


## 2.3. Discriminant Analysis

The idea of discriminant analysis is to model the distribution of the predictorsXj⋅ separately for each of the response classes Yj, and then to use the Bayes theorem to describe the conditional probabilities PYj=i|Xj⋅=xj⋅ (e.g., James et al. [15]).Specifically, letπi≜PYj=i denote the probability that the j th subject is randomly selected from class i so that ∑i=1Iπi=1. Moreover, applying the Bayes theorem to the conditional density function fj|ixj⋅ and π gives the posterior probability as(6)PYj=i|Xj⋅=xj⋅=fj|ixj⋅πi∑l=1Ifj|lxj⋅πl,for i=1,…,I and j=1,…,n.To compare two classesi and l with i≠l, we calculate the log-ratio of (6), given by(7)logPYj=i|Xj⋅=xj⋅PYj=l|Xj⋅=xj⋅=logfj|ixj⋅fj|lxj⋅+logπiπl.With a distribution assumption (1) imposed, we now have two scenarios for the first term in the right-hand side of (7), say logfj|ixj⋅/fj|lxj⋅.Scenario 1.
If the covariance matricesΣi in (1) are assumed to be common, that is, Σi=Σ for every i with Σ being a positive definite matrix, (7) becomes(8)logπiπl−12μi+μl⊤Σ−1μi+μl+xj⋅⊤Σ−1μi+μl.
If equation (8) >0, then PYj=i|Xj⋅=xj⋅>PYj=l|Xj⋅=xj⋅, showing that subject j with predictors Xj⋅=xj⋅ is more likely selected from class i than from class l. Consequently, (8) defines a boundary between classes i and l in the sense that there is a linear function in xj⋅ separating classes i and l.
Motivated by the form of (8), we consider a linear function in x as(9)δix≜logπi−12μi⊤Σ−1μi+x⊤Σ−1μi.
Moreover,μi and πi can be empirically estimated, respectively, as(10)μ^i=1ni∑Yj=ixj⋅,π^i=nin.
For the estimation ofΣ−1, or equivalently Θ, we adopt (5) by pooling all subjects in the dataset and denote Θ˜ as the estimator of Θ. Therefore, (9) can be estimated as(11)δ˜ix=logπ^i−12μ^i⊤Θ˜μ^i+x⊤Θ˜μ^i,and we call (11) the network-based linear discriminant function (NetLDA) and it is used to determine the class label for a new observation. For the prediction of a new subject with the predictor x˜, we first calculate δ˜ix˜ using (11) for i=1,…,I. Next, we find i∗ that is defined as(12)i∗=argmaxi=1,…,Iδ˜ix˜;and the class label for this subject is then predicted as i∗.Scenario 2.
We allowΣi≠Σl, or equivalently, Θi≠Θl, for any i≠l and i,l=1,…,I. Then, under a distribution assumption (1), we have(13)logfj|ixj⋅fj|lxj⋅=log1/2πp/2Σi1/2exp−1/2xj⋅−μi⊤Θixj⋅−μi1/2πp/2Σl1/2exp−1/2xj⋅−μl⊤Θlxj⋅−μl=logΘl−1/2Θi−1/2+12xj⋅−μl⊤Θlxj⋅−μl−xj⋅−μi⊤Θixj⋅−μi.
Replacing the first term in the right-hand side of (7) by (13) yields(14)logπiπl+logΘl−1/2Θi−1/2+12xj⋅−μl⊤Θlxj⋅−μl−xj⋅−μi⊤Θixj⋅−μi.
Therefore, based on (14), we further define a quadratic function of x based on the class i:(15)φix≜logπi+12logΘi−12x−μi⊤Θix−μi.
Fori=1,…,I, the estimator of the precision matrix Θi, denoted as Θ^i, is obtained by (5) based on the predictor information in class i; μi and πi can be estimated by (10). Therefore, (15) can be estimated by(16)φ^ix=logπ^i+12logΘ^i−12x−μ^i⊤Θ^ix−μ^i.
The function (16) is called the network-based quadratic discriminant function (NetQDA) and is used to determine the class label for a new observation. For the prediction of a new subject with the predictor x˜, we first calculate φ^ix˜ using (16) for i=1,…,I. Next, we find i∗ that is defined as(17)i∗=argmaxi=1,…,Iφ^ix˜;and the class label for this subject is then predicted as i∗.


## 2.4. Assessment of Classification and Prediction

We first introduce micro-averaged metrics (e.g., Chen et al. [7]). Let V and T represent the classes of the subject indexes for validation and training datasets, respectively. Let n=T and m=V denote the sizes of the training and validation data, respectively. We use the training data T to fit models and then apply fitted models to compute predicted values Y^j for j∈V. After that, for i=1,…,I, we calculate the number of the true positives (TP), the number of the false positives (FP), and the number of the false negatives (FN) under the validation data V, respectively:(18)TPi=∑j∈VIYj=i,Y^j=i,FPi=∑j∈VIYj≠i,Y^j=i,FNi=∑j∈VIYj=i,Y^j≠i,where I⋅ is the indicator function.We define precision (PRE) and recall (REC) under the validation dataV, respectively, as(19)PRE=∑i=1ITPi∑i=1ITPi+∑i=1IFPi,REC=∑i=1ITPi∑i=1ITPi+∑i=1IFNi.Then, micro-F-score is defined as(20)F=2×PRE×RECPRE+REC.According to definitions in (19), when all subjects are correctly classified, FP and FN are equal to zero, yielding that PRE and REC are equal to one; if all subjects are falsely classified, then TP is equal to zero, and thus, PRE and REC are equal to zero. Therefore, values of PRE and REC are between zero and one. Moreover, under the range 0,1, the F-score falls in 0,1 as well by treating 0/0 as zero. In principle, the higher values of PRE, REC, and F-score reflect the better performance and the more accurate classification (e.g., Chen et al. [7]).In addition to criteria above, the other commonly used criterion is the adjusted Rand index (ARI). Fori,i′=1,…,I and under the validation data V, we define(21)nii′=∑j∈VIYj=i,Y^j=i′.Moreover, we defineai=∑i′=1Inii′ for i=1,…,I and bi′=∑i=1Inii′ for i′=1,…,I. Then, ARI under the validation data V is defined as Hubert and Arabie [16].(22)ARI=∑i,i′nii′2−∑iai2∑i′bi′2/n2∑iai2+∑i′bi′2/2−∑iai2∑i′bi′2/n2.As mentioned in Hubert and Arabie [16], ARI is bounded above by one, and the higher value of ARI indicates the more accurate classification.


## 2.5. Benchmark of NetDA

The conventional linear discriminant methods (e.g., MASS, sparseLDA, and penalizedLDA) aim to adopt (9) with Θ estimated by the inverse of the empirical estimator of Σ. However, this approach may encounter cumbersome computation or possible singularity when calculating inverse matrices. In addition, if Θ is sparse, all entries in empirical estimators of Θ are nonzero, and it indicates that some unconnected pairs of predictors may be falsely included. As a result, imprecise estimator of Θ may implicitly affect the performance of classification.Unlike existing methods, the first contribution of NetDA is to estimateΘ directly. The graphical lasso method is a tool to identify zero entries in Θ and estimate nonzero ones. This approach enables us to retain connected pairs of predictors and exclude unconnected ones. Moreover, our implementation can avoid computing inverse of Σ.The second contribution of NetDA is to handle heterogeneous network structure that stratifies the predictor information by class when characterizing the predictor network structures, and is able to deal with multiclassification by adopting network structure in each class.


## 3. Details of NetDA

In this section, we first overview the technique that we need in the existing package. After that, we describe functions in our package.


### 3.1. Library Overview

In our package, we use the package glasso in R software. Specifically, the package glasso follows from the graphical lasso method proposed by Friedman et al. [13]. The purpose of glasso is to detect network structures of random vectors that follow multivariate normal distributions. In particular, under multivariate normal distributions and operations in glasso, the detection of network structures is equivalent to the estimation of precision matrix.


### 3.2. WineData

In our package, we take the wine dataset as an example, which is available inhttps://archive.ics.uci.edu/ml/datasets/wine. These data were collected based on a chemical analysis of wines grown in the same region in Italy but derived from three different cultivars. In this dataset, there are three types of wines and 13 constituents, including alcohol (Alcohol), malic acid (Malic acid), ash (Ash), alkalinity of ash (Alkalinity), magnesium (Magnesium), total phenols (phenols), flavanoids (Flavanoids), nonflavanoid phenols (Nonflavanoid), proanthocyanins (Proanthocyanins), color intensity (Color), hue (Hue), OD280/OD315 of diluted wines (OD280), and proline (Proline).In the following analysis, the response is types of wines that are labeled as 1, 2, and 3; constituents are treated as predictors that are continuous. The goal is to adopt the information of constituents to construct predictive models and then use them to classify type of wines for a given subject.


### 3.3. NetDA

NetDA contains two methods. The first method is called NetLDA, which aims to estimate the precision matrix by pooling all individuals in the data, and the corresponding discriminant function is given by (11). The second approach is called NetQDA, whose strategy is to estimate precision matrices based on individuals in different classes, and then use class-dependent estimated precision matrices to define quadratic discriminant functions in (16). Unlike NetLDA, NetQDA takes possibly class-dependent network structures of the predicted variables into account and uses network structures in different classes to determine which classes individuals belong to. When either linear discriminant functions or quadratic discriminant functions are obtained, they can be used to determine the class for a new subject.To implement the NetLDA and NetQDA methods, we use the following command:NetDA (X, Y, method, X_test) where the meaning of each argument is described as follows:(i)
X: this is an n×p matrix of the predictors from the training data(ii)
Y: this is an n-dimensional vector of the response from the training data, whose elements are positive integers and reflect class-labels(iii)
Method: it is a scalar to determine the classification method: method = 1 represents NetLDA in (11), and method = 2 represents NetQDA in (16)(iv)
X_test: this is an m×p matrix of the predictors from the validation dataThe purpose of NetDA is to apply the training data“X” and “Y” to determine a fitted model that is specified by the argument “method.” After that, we use “X_test and a fitted model to determine the predicted class for subjects in the validation data. Therefore, the function NetDA returns a list of components:(i)
yhat: it is a vector of predicted responses obtained by NetLDA or NetQDA based on the predictors in the validation data (X_test).(ii)
Network: this is the estimators of precision matrices. If “method = 1” is chosen, then there is one precision matrix; if “method = 2” is given, then there areI precision matrices.


### 3.4. Metrics

The function Metrics is utilized to assess the performance of classification and prediction based on some commonly used criteria that are introduced in the Section 2.4. Specifically, given responses from the validation data and predicted values obtained by NetDA, we first derive a confusion matrix to see the classification result. To further assess the performance of prediction, we evaluate precision, recall,F-score, and ARI defined in (19), (20), and (22), respectively.To obtain the desired results, we use the following command:Metrics (yhat,Y_test) where the meaning of each argument is described as follows:(i)
yhat: this is an m-dimensional vector of the predicted responses determined by NetDA or other methods(ii)
Y_test: this is an m-dimensional vector of the response from the validation dataThe function metrics returns a list of components:(i)
Confusion matrix: a confusion matrix based on predicted values (yhat) and responses from the validation data (Y_test)(ii)
(PRE, REC,F-score): values of precision, recall, and F-score defined in (19) and (20), respectively(iii)
ARI: values of the ARI defined in (22)


## 3.1. Library Overview

In our package, we use the package glasso in R software. Specifically, the package glasso follows from the graphical lasso method proposed by Friedman et al. [13]. The purpose of glasso is to detect network structures of random vectors that follow multivariate normal distributions. In particular, under multivariate normal distributions and operations in glasso, the detection of network structures is equivalent to the estimation of precision matrix.


## 3.2. WineData

In our package, we take the wine dataset as an example, which is available inhttps://archive.ics.uci.edu/ml/datasets/wine. These data were collected based on a chemical analysis of wines grown in the same region in Italy but derived from three different cultivars. In this dataset, there are three types of wines and 13 constituents, including alcohol (Alcohol), malic acid (Malic acid), ash (Ash), alkalinity of ash (Alkalinity), magnesium (Magnesium), total phenols (phenols), flavanoids (Flavanoids), nonflavanoid phenols (Nonflavanoid), proanthocyanins (Proanthocyanins), color intensity (Color), hue (Hue), OD280/OD315 of diluted wines (OD280), and proline (Proline).In the following analysis, the response is types of wines that are labeled as 1, 2, and 3; constituents are treated as predictors that are continuous. The goal is to adopt the information of constituents to construct predictive models and then use them to classify type of wines for a given subject.


## 3.3. NetDA

NetDA contains two methods. The first method is called NetLDA, which aims to estimate the precision matrix by pooling all individuals in the data, and the corresponding discriminant function is given by (11). The second approach is called NetQDA, whose strategy is to estimate precision matrices based on individuals in different classes, and then use class-dependent estimated precision matrices to define quadratic discriminant functions in (16). Unlike NetLDA, NetQDA takes possibly class-dependent network structures of the predicted variables into account and uses network structures in different classes to determine which classes individuals belong to. When either linear discriminant functions or quadratic discriminant functions are obtained, they can be used to determine the class for a new subject.To implement the NetLDA and NetQDA methods, we use the following command:NetDA (X, Y, method, X_test) where the meaning of each argument is described as follows:(i)
X: this is an n×p matrix of the predictors from the training data(ii)
Y: this is an n-dimensional vector of the response from the training data, whose elements are positive integers and reflect class-labels(iii)
Method: it is a scalar to determine the classification method: method = 1 represents NetLDA in (11), and method = 2 represents NetQDA in (16)(iv)
X_test: this is an m×p matrix of the predictors from the validation dataThe purpose of NetDA is to apply the training data“X” and “Y” to determine a fitted model that is specified by the argument “method.” After that, we use “X_test and a fitted model to determine the predicted class for subjects in the validation data. Therefore, the function NetDA returns a list of components:(i)
yhat: it is a vector of predicted responses obtained by NetLDA or NetQDA based on the predictors in the validation data (X_test).(ii)
Network: this is the estimators of precision matrices. If “method = 1” is chosen, then there is one precision matrix; if “method = 2” is given, then there areI precision matrices.


## 3.4. Metrics

The function Metrics is utilized to assess the performance of classification and prediction based on some commonly used criteria that are introduced in the Section 2.4. Specifically, given responses from the validation data and predicted values obtained by NetDA, we first derive a confusion matrix to see the classification result. To further assess the performance of prediction, we evaluate precision, recall,F-score, and ARI defined in (19), (20), and (22), respectively.To obtain the desired results, we use the following command:Metrics (yhat,Y_test) where the meaning of each argument is described as follows:(i)
yhat: this is an m-dimensional vector of the predicted responses determined by NetDA or other methods(ii)
Y_test: this is an m-dimensional vector of the response from the validation dataThe function metrics returns a list of components:(i)
Confusion matrix: a confusion matrix based on predicted values (yhat) and responses from the validation data (Y_test)(ii)
(PRE, REC,F-score): values of precision, recall, and F-score defined in (19) and (20), respectively(iii)
ARI: values of the ARI defined in (22)


## 4. Demonstration of NetDA

In this section, we demonstrate standard analysis of classification and prediction based on two functions in the package. To show the advantage of NetDA, we compare with existing packages MASS, sparseLDA, and penalizedLDA.


### 4.1. Simulation

LetI=3 denote the number of classes, and let p=12 denote the dimension of predictors. We specify the sample size n=600, in which the size of the ith class is given by ni=n/I for i=1,…,I. For each class, we consider different network structures of predictors. Specifically, for class i=1,…,I, let Ωo,i denote a p×p matrix whose diagonal entries are zero and off-diagonal entries are specified as either one or zero to reflect edges of the corresponding two nodes in Figure 1. That is, for s≠t, entry s,t in Ωo,i is 1 if the edge exists between X⋅s and X⋅t and 0 otherwise. In addition, we further define a p×p diagonal matrix Ωd,i whose nonzero entries are taken as the common value 0.1+ΛminΩo,i, where ΛminΩo,i represents the smallest eigenvalue of Ωo,i. Finally, we define the precision matrix as Θi≜Ωo,i+Ωd,i that is invertible. Therefore, based on Gaussian graphical models, the p-dimensional vector of predictors in class i is generated from a multivariate normal distribution with mean zero and the covariance matrix Σi=Θi−1 for i=1,…,I.Figure 1
Predictor network structures for four different classes.Let Theta1, Theta2, and Theta3 denotep×p matrices that reflect network structures in left, middle, and right panels of Figure 1, respectively, and one can specify those three matrices as follows:>Theta1(23)1,,2,3,4,5,6,7,8,9,10,11,121,0110000000002,1001000000003,1001100000004,0110010000005,0010010000006,0001100000007,0000000110008,0000001001009,00000010011010,00000001100111,00000000100112,000000000110.> Theta2(24),1,2,3,4,5,6,7,8,9,10,11,121,0110000000002,1001000000003,1001000000004,0110000000005,0000011000006,0000100100007,0000100100008,0000011000009,00000000011010,00000000100111,00000000100112,000000000110.> Theta3(25),1,2,3,4,5,6,7,8,9,10,11,121,0100000000002,1010000000003,0100000000004,0000100000005,0001010000006,0000100000007,0000000100008,0000001010009,00000001000010,00000000001011,00000000010112,000000000010.Following the description above, we generate then×p dimensional matrix for predictors and then determine the simulated data:>Theta1 = Theta1 + diag (0.1 + abs (min (eigen (Theta1) $value)),p)>Sigma1 = cov2cor (solve (Theta1))>X1 = mvrnorm (n = (n/I), rep (0, p), Sigma1, tol = 1e − 6, empirical = FALSE)>>Theta2 = Theta2 + diag (0.1 + abs (min (eigen (Theta2) $value)),p)>Sigma2 = cov2cor (solve (Theta2))>X2 = mvrnorm (n = (n/I), rep (0, p), Sigma2, tol = 1e − 6, empirical = FALSE)>>Theta3 = Theta3 + diag (0.1 + abs (min (eigen (Theta3) $value)),p)>Sigma3 = cov2cor (solve (Theta3))>X3 = mvrnorm (n = (n/I), rep (0, p), Sigma3, tol = 1e − 6, empirical = FALSE)>data = cbind (c (rep (1, n/I), rep (2, n/I), rep (3, n/I)), rbind (X1, X2, X3))To perform classification, we implement NetDA function with two different scenarios described in Section 2.3. In addition, we demonstrate three existing methods labeled as lda (MASS), sda (sparseLDA), and pda (penalizedLDA), respectively. Detailed descriptions are given below:>Y = data [, 1]>X = data [, 2 : 13]>#Demonstration of MASS>lda = lda (Y·X, prior = c (length (which (Y = = 1)), length (which (Y = = 2)), ++length (which (Y = = 3)))/length (Y))>yhat_lda = predict (lda, data.frame (X)) $class>>#Demonstration of sparseLDA>y = matrix (0, n, I)>y [1 : 200, 1] = 1>y [201 : 400, 2] = 1>y [401 : 600, 3] = 1>colnames (y) < −c (“1,” “2,” “3”)>sda = sda (data.frame (X), y, lambda = 1e − 6, stop = −1, maxIte = 25, +trace = TRUE)>yhat_sda = as.numeric (unlist (predict (sda, data.frame (X)) $class))>#Demonstration of penalizedLDA>pda = PenalizedLDA (X, Y, lambda = 0.14, K = 2)>yhat_pda = as.numeric (unlist (predict (pda, data.frame (X)))) [1 : n]>#Demonstration of NetDA>yhat_netlda = NetDA (X, Y, method = 1, X) $yhat>yhat_netqda = NetDA (X, Y, method = 2, X) $yhatAfter that, to assess the performance of classification, we adopt the function Metrics to compute values of criteria (19) and (20) as shown by [2] and (22) indicated by [3].>F_lda = Metrics (yhat_lda,Y) [2]>F_sda = Metrics (yhat_sda,Y) [2]>F_pda = Metrics (yhat_pda,Y) [2]>F_netlda = Metrics (yhat_netlda,Y) [2]>F_netqda = Metrics (yhat_netqda,Y) [2]>ARI_lda = Metrics (yhat_lda,Y) [3]>ARI_sda = Metrics (yhat_sda,Y) [3]>ARI_pda = Metrics (yhat_pda,Y) [3]>ARI_netlda = Metrics (yhat_netlda,Y) [3]>ARI_netqda = Metrics (yhat_netqda,Y) [3]We repeat above simulations 500 times and summarize numerical results in Table1. We observe that the package NetDA provides higher values of PRE, REC, F-score, and adn ARI, showing that the classification obtained by the NetDA method is more accurate than that determined by other methods. Specifically, compare with MASS and NetLDA, we can see that the latter outperforms the former method, which is due to the incorporation of network structure with irrelevant pairs of predictors removed from Θ. On the other hand, for the comparison between NetLDA and NetQDA, we can see that the NetQDA is much better than the NetLDA method, because the NetQDA method successfully detects network structures from each class, and those detected network structures are valid to do classification. Those numerical findings verify the discussion in Subsection “Benchmark of NetDA.”Table 1
Simulation results.
MASSsparseLDApenalizedLDANetDA (NetLDA)NetDA (NetQDA)PRE0.6270.6860.6690.7290.891REC0.6270.6860.6690.7290.891F-score0.6270.6860.6690.7290.891ARI0.5210.5430.5300.6220.707


### 4.2. Real Data Analysis

In this study, we take the wine dataset as an example, which is introduced in Section3, to demonstrate the package NetDA. To demonstrate the functions and perform classification and prediction, we first split the full data into the training data and the validation data. In our example, we take the first 45 samples in each class to obtain the training data and use the remaining samples in each class to form the validation data.>data (WineData)>Y = WineData [, 1] #the response>X = WineData [, 2 : 14] #the predictors>D1 = WineData [which (Y = = 1), ]>D2 = WineData [which (Y = = 2), ]>D3 = WineData [which (Y = = 3), ]>#An example of user-specific training data and validation data>#“Train” represents the training data and “Test” represents validation data in our example.>Train = rbind (D1 [1 : 45, ], D2 [1 : 45, ], D3 [1 : 45, ])>Test = rbind (D1 [46 : dim (D1) [1], ], D2 [46 : dim (D2) [1], ], D3 [46 : dim (D3) [1], ])>#The response (Y) and predictors (X) in the training data>X = Train [, 2 : 14] > Y = Train [, 1]>#The response (Y_test) and predictors (X_test) in the validation data>X_test = Test [, 2 : 14] > Y_test = Test [, 1]When the training data and the validation data are determined, we employ the function NetDA to perform classification. We insert“X,”“Y,” and “X_test” to the function NetDA, and we denote “NetLDA” and “NetQDA” as the argument method = 1 and method = 2, respectively. The resulting vectors of predicted classes and estimated precision matrices are given by “$yhat” and “$Network,” respectively.>NetDA (X, Y, method = 1, X_test) -> NetLDA>yhat_lda = NetLDA$yhat>Net_lda = NetLDA$Network>yhat_lda[1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 2 2 2 2 2 2 2[39] 2 2 2 2 3 3 3 3>round (Net_lda, 3)(26),1,2,3,4,5,6,7,8,9,10,11,12,131,3.739−0.415−0.6650.1070.0040.000−0.4440.0000.000−0.5080.000−0.291−0.0042,−0.4151.247−0.511−0.0350.0000.1090.324−0.2250.0650.0001.4130.0170.0003,−0.665−0.51119.463−0.936−0.059−0.055−0.6190.0000.000−0.0870.0000.000−0.0054,0.107−0.035−0.9360.1670.001−0.0480.152−0.160−0.038−0.023−0.035−0.0130.0005,0.0040.000−0.0590.0010.0060.0000.0070.087−0.0360.0030.0070.0130.0006,0.0000.109−0.055−0.0480.0008.774−3.8470.000−0.364−0.2490.000−1.274−0.0017,−0.4440.324−0.6190.1520.007−3.8476.1272.104−1.8000.016−2.845−1.877−0.0028,0.000−0.2250.000−0.1600.0870.0002.10446.8790.000−0.0840.0000.762−0.0029,0.0000.0650.000−0.038−0.036−0.365−1.8000.0005.582−0.2160.000−0.3910.00110,−0.5080.000−0.087−0.0230.003−0.2490.016−0.084−0.2160.4671.6630.6060.00011,0.0001.4130.000−0.0350.0070.000−2.8450.0000.0001.66232.8500.000−0.00212,−0.2910.0170.000−0.0130.013−1.274−1.8770.762−0.3910.6060.0005.667−0.00113,−0.0040.000−0.0050.0000.000−0.001−0.002−0.0020.0010.000−0.002−0.0010.000.#############>NetDA (X, Y, method = 2, X_test) -> NetQDA>yhat_qda = NetQDA$yhat>Net_qda = NetQDA$Network>yhat_qda[1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2[42] 2 3 3 3 3>round (Net_qda [[1]], 3)(27),1,2,3,4,5,6,7,8,9,10,11,12,131,7.1760.0110.0000.470−0.0730.000−0.4040.000−1.754−0.1750.000−0.159−0.0052,0.0122.7210.0000.005−0.0070.0000.0000.051−0.1710.0002.0390.0000.0033,0.0000.00023.194−0.959−0.1610.0000.0000.0000.0000.0000.0000.0000.0004,0.4700.005−0.9590.236−0.0110.1380.076−0.397−0.1010.046−0.2440.1700.0005,−0.073−0.007−0.161−0.0110.013−0.0970.008−0.0850.049−0.042−0.029−0.0230.0006,0.0000.0000.0000.138−0.09722.439−7.4720.000−0.680−2.3260.000−0.3300.0027,−0.4050.0000.0000.0760.008−7.47116.5200.000−3.594−2.1310.0000.0000.0018,0.0000.0510.000−0.397−0.0850.0000.00069.3560.0000.0560.0000.1590.0009,−1.754−0.1710.000−0.1010.049−0.680−3.5940.0008.449−0.2460.0000.0000.00310,−0.1750.0000.0000.046−0.042−2.326−2.1310.056−0.2462.5090.5660.000−0.00611,0.0002.0390.000−0.244−0.0290.0000.0000.0000.0000.56348.1390.000−0.00812,−0.1580.0000.0000.170−0.0230.3280.0000.1590.0000.0000.0007.2130.00413,−0.0050.0030.0000.0000.0000.0020.0010.0000.003−0.006−0.0080.0040.000.>round (Net_qda [[2]], 3)(28),1,2,3,4,5,6,7,8,9,10,11,12,131,4.0600.0190.710−0.128−0.0050.000−0.5270.0000.432−0.7830.0000.169−0.0012,0.0212.3570.000−0.173−0.0170.018−0.015−1.439−0.3590.1380.9120.1220.0043,0.7080.00018.423−1.213−0.0440.000−0.0060.0000.0000.0000.0000.8260.0044,−0.128−0.173−1.2130.2080.0070.082−0.043−0.2140.0590.056−0.160−0.393−0.0015,−0.005−0.017−0.0440.0070.005−0.0330.0630.090−0.0660.000−0.0490.0140.0006,0.0000.0170.0000.082−0.0337.713−5.9990.8661.7320.4430.000−0.822−0.0027,−0.527−0.015−0.006−0.0430.063−5.99910.0170.181−3.932−1.565−1.637−1.9020.0028,0.0001.4400.000−0.2140.0900.8690.17346.3670.0000.0000.0005.077−0.0039,0.432−0.3590.0000.059−0.0661.733−3.9320.0005.2820.5060.000−0.648−0.00110,−0.7830.1370.0000.0560.0000.443−1.5650.0000.5061.7330.5420.763−0.00111,0.0000.9120.000−0.160−0.0490.000−1.6340.0000.0000.54229.1090.0000.00812,0.1700.1210.826−0.3930.014−0.822−1.9025.077−0.6480.7630.0007.2930.00213,−0.0010.0040.004−0.0010.000−0.0020.002−0.003−0.001−0.0010.0080.0020.000.>round (Net_qda [[3]], 3)(29),1,2,3,4,5,6,7,8,9,10,11,12,131,4.261−0.2540.000−0.0720.0200.0000.0000.000−1.200−0.202−0.098−0.8320.0032,−0.2540.9510.000−0.0960.0170.3500.3360.0000.3020.0500.0000.1460.0003,0.0000.00039.923−2.076−0.119−1.7130.0000.0000.0000.0000.000−1.6400.0064,−0.072−0.096−2.0760.359−0.010−0.326−0.073−0.193−0.115−0.010−0.178−0.0550.0005,0.0200.017−0.119−0.0100.0220.060−0.4310.371−0.0110.021−0.0340.040−0.0016,0.0000.350−1.713−0.3260.06010.790−0.498−0.740−4.4610.0550.000−1.1730.0007,0.0000.3360.000−0.073−0.431−0.49824.0980.000−1.503−0.9300.0002.4000.0188,0.0000.0000.000−0.1930.371−0.7400.00050.6220.0000.0000.0000.000−0.0149,−1.2000.3020.000−0.115−0.011−4.461−1.5030.00012.612−1.0320.0000.846−0.00810,−0.2020.0500.000−0.0100.0210.055−0.9300.000−1.0320.4101.3530.070−0.00111,−0.0980.0000.000−0.178−0.0340.0000.0000.0000.0001.35350.6990.000−0.00512,−0.8320.146−1.640−0.0550.040−1.1732.4000.0000.8460.0700.00013.312−0.00813,0.0030.0000.0060.000−0.0010.0000.018−0.014−0.008−0.001−0.005−0.0080.000.Moreover, for the visualization of estimated network structures, we further apply the packages network, GGally and sna, to draw the network structure based on the estimated precision matrices. The following commands are implemented to draw network structures, and the corresponding figures determined by NetDA with the arguments method = 1 and method = 2 are displayed in Figures2 and 3, respectively.>library (network)>library (GGally)>library (sna)>material_name =c (“Alcohol,” “Malic acid,” “Ash,” “Alkalinity,” “Magnesium,” “phenols,” + “Flavanoids,” “Nonflavanoid,” “Proanthocyanins,” “Color,” “Hue,” “OD280,” “Proline”)>ggnet2 (Net_lda, mode = “circle,” size = 8, label = material_name, label.size = 5)>ggnet2 (Net_qda [[1]], mode = “circle,” size = 8, label = material_name, label.size = 5)>ggnet2 (Net_qda [[2]], mode = “circle,” size = 8, label = material_name, label.size = 5)>ggnet2 (Net_qda [[3]], mode = “circle,” size = 8, label = material_name, label.size = 5)Figure 2
Pooled subject-based network structures obtained by NetDA under the argument method = 1.Figure 3
Class-dependent network structures obtained by NetDA under the argument method = 2. Left, middle, and right panels are based on predictors in classes 1, 2, and 3, respectively.From Figures2 and 3, we can observe that precision matrices provide complex network structures in predictors. In particular, in Figure 3, we can see that the estimated class-dependent network structures are different from each other, and the network structure in class 2 looks more complex than others. To assess the performance of prediction, we input predicted values (yhat_lda or yhat_qda) and responses in the validation data (Y_test) to the function Metrics, and the resulting values are displayed below.>Metrics (yhat_lda,Y_test)$“Confusion matrix”[, 1] [, 2] [, 3][1, ] 15 0 0[2, ] 0 26 0[3, ] 0 1 4$“(PRE, REC,F-score)”[1] 0.9782609 0.9782609 0.9782609$ARI[1] 0.9410827#############>Metrics (yhat_qda,Y_test)$“Confusion matrix”[, 1] [, 2] [, 3][1, ] 15 0 0[2, ] 0 27 0[3, ] 0 0 4$“(PRE, REC,F-score)”[1] 1 1 1$ARI[1] 1Finally, we further adopt the function lda in the packages MASS, sparseLDA, and penalizedLDA to perform the conventional discriminant methods and compare them with our NetDA. Detailed implementations and numerical results are given below:>Wine = data.frame (cbind (Y, X))>##Demonstration of MASS>lda = lda (Y., Wine, prior = c (length (which (Y = = 1)), length (which (Y = = 2)), +length (which (Y = = 3)))/length (Y))>predict (lda,X_test) $class -> lda_pred>>Metrics (lda_pred,Y_test)$“Confusion matrix”[, 1] [, 2] [, 3][1, ] 15 1 0[2, ] 0 26 0[3, ] 0 0 4$“(PRE, REC,F-score)”[1] 0.9782609 0.9782609 0.9782609$ARI[1] 0.9196658>#Demonstration of sparseLDA>n = length (Y)>I = max (Y)>y = matrix (0, n, I)>y [1 : 45, 1] = 1>y [46 : 90, 2] = 1>y [91 : 135, 3] = 1>colnames (y) <- c (“1,” “2,” “3”)>sda = sda (data.frame (X), y, lambda = 1e − 6, stop = −1, maxIte = 25, +trace = TRUE)ite: 1 ridge cost: 101.1419b1: 0.001636884ite: 2 ridge cost: 47.61608b1: 0.002647311ite: 3 ridge cost: 47.61608b1: 0.002647311ite: 1 ridge cost: 129.1129b1: 0.01364973ite: 2 ridge cost: 129.1129b1: 0.01364973final update, total ridge cost: 176.729b1: 0.01629704>sda_pred = as.numeric (unlist (predict (sda, data.frame (X_test)) $class))>Metrics (sda_pred,Y_test)$“Confusion matrix”[, 1] [, 2] [, 3][1, ] 14 26 2[2, ] 0 0 0[3, ] 0 0 1$“(PRE, REC,F-score)”[1] 0.3488372 0.3488372 0.3488372$ARI[1] 0.07272024>#Demonstration of penalizedLDA>pda = PenalizedLDA (X, Y, lambda = 0.14, K = 2)>pda_pred = as.numeric (unlist (predict (pda, data.frame (X_test)))) [1 : 46]>Metrics (pda_pred,Y_test)$“Confusion matrix”[, 1] [, 2] [, 3][1, ] 14 4 0[2, ] 0 21 0[3, ] 0 1 3$“(PRE, REC,F-score)”[1] 0.8837209 0.8837209 0.8837209$ARI[1] 0.6229476In general, we can see that NetLDA and NetQDA have the satisfactory performance in prediction. For the NetLDA method, there is one misclassification as shown in the confusion matrix, while the predicted classes determined by NetQDA are all equal to the responses in the validation data. From the comparison to NetDA, we observe from a confusion matrix determined by the conventional linear discriminant analysis (lda) is comparable to that obtained by the NetLDA method, but it is interesting to see that the value of ARI determined by NetLDA is slightly larger than that based on lda. In addition, it is clear to see that the NetQDA method is better than lda. On the contrary, it is surprising to see that two penalized methods sparseLDA and penalizedLDA do not have satisfactory performance of classification and prediction, especially that sparseLDA has the most unexpected result. In summary, the numerical results in this data analysis show (a) the importance of incorporating predictor network structures in the classification procedure, and (b) the advantage of adopting class-dependent network structures.


## 4.1. Simulation

LetI=3 denote the number of classes, and let p=12 denote the dimension of predictors. We specify the sample size n=600, in which the size of the ith class is given by ni=n/I for i=1,…,I. For each class, we consider different network structures of predictors. Specifically, for class i=1,…,I, let Ωo,i denote a p×p matrix whose diagonal entries are zero and off-diagonal entries are specified as either one or zero to reflect edges of the corresponding two nodes in Figure 1. That is, for s≠t, entry s,t in Ωo,i is 1 if the edge exists between X⋅s and X⋅t and 0 otherwise. In addition, we further define a p×p diagonal matrix Ωd,i whose nonzero entries are taken as the common value 0.1+ΛminΩo,i, where ΛminΩo,i represents the smallest eigenvalue of Ωo,i. Finally, we define the precision matrix as Θi≜Ωo,i+Ωd,i that is invertible. Therefore, based on Gaussian graphical models, the p-dimensional vector of predictors in class i is generated from a multivariate normal distribution with mean zero and the covariance matrix Σi=Θi−1 for i=1,…,I.Figure 1
Predictor network structures for four different classes.Let Theta1, Theta2, and Theta3 denotep×p matrices that reflect network structures in left, middle, and right panels of Figure 1, respectively, and one can specify those three matrices as follows:>Theta1(23)1,,2,3,4,5,6,7,8,9,10,11,121,0110000000002,1001000000003,1001100000004,0110010000005,0010010000006,0001100000007,0000000110008,0000001001009,00000010011010,00000001100111,00000000100112,000000000110.> Theta2(24),1,2,3,4,5,6,7,8,9,10,11,121,0110000000002,1001000000003,1001000000004,0110000000005,0000011000006,0000100100007,0000100100008,0000011000009,00000000011010,00000000100111,00000000100112,000000000110.> Theta3(25),1,2,3,4,5,6,7,8,9,10,11,121,0100000000002,1010000000003,0100000000004,0000100000005,0001010000006,0000100000007,0000000100008,0000001010009,00000001000010,00000000001011,00000000010112,000000000010.Following the description above, we generate then×p dimensional matrix for predictors and then determine the simulated data:>Theta1 = Theta1 + diag (0.1 + abs (min (eigen (Theta1) $value)),p)>Sigma1 = cov2cor (solve (Theta1))>X1 = mvrnorm (n = (n/I), rep (0, p), Sigma1, tol = 1e − 6, empirical = FALSE)>>Theta2 = Theta2 + diag (0.1 + abs (min (eigen (Theta2) $value)),p)>Sigma2 = cov2cor (solve (Theta2))>X2 = mvrnorm (n = (n/I), rep (0, p), Sigma2, tol = 1e − 6, empirical = FALSE)>>Theta3 = Theta3 + diag (0.1 + abs (min (eigen (Theta3) $value)),p)>Sigma3 = cov2cor (solve (Theta3))>X3 = mvrnorm (n = (n/I), rep (0, p), Sigma3, tol = 1e − 6, empirical = FALSE)>data = cbind (c (rep (1, n/I), rep (2, n/I), rep (3, n/I)), rbind (X1, X2, X3))To perform classification, we implement NetDA function with two different scenarios described in Section 2.3. In addition, we demonstrate three existing methods labeled as lda (MASS), sda (sparseLDA), and pda (penalizedLDA), respectively. Detailed descriptions are given below:>Y = data [, 1]>X = data [, 2 : 13]>#Demonstration of MASS>lda = lda (Y·X, prior = c (length (which (Y = = 1)), length (which (Y = = 2)), ++length (which (Y = = 3)))/length (Y))>yhat_lda = predict (lda, data.frame (X)) $class>>#Demonstration of sparseLDA>y = matrix (0, n, I)>y [1 : 200, 1] = 1>y [201 : 400, 2] = 1>y [401 : 600, 3] = 1>colnames (y) < −c (“1,” “2,” “3”)>sda = sda (data.frame (X), y, lambda = 1e − 6, stop = −1, maxIte = 25, +trace = TRUE)>yhat_sda = as.numeric (unlist (predict (sda, data.frame (X)) $class))>#Demonstration of penalizedLDA>pda = PenalizedLDA (X, Y, lambda = 0.14, K = 2)>yhat_pda = as.numeric (unlist (predict (pda, data.frame (X)))) [1 : n]>#Demonstration of NetDA>yhat_netlda = NetDA (X, Y, method = 1, X) $yhat>yhat_netqda = NetDA (X, Y, method = 2, X) $yhatAfter that, to assess the performance of classification, we adopt the function Metrics to compute values of criteria (19) and (20) as shown by [2] and (22) indicated by [3].>F_lda = Metrics (yhat_lda,Y) [2]>F_sda = Metrics (yhat_sda,Y) [2]>F_pda = Metrics (yhat_pda,Y) [2]>F_netlda = Metrics (yhat_netlda,Y) [2]>F_netqda = Metrics (yhat_netqda,Y) [2]>ARI_lda = Metrics (yhat_lda,Y) [3]>ARI_sda = Metrics (yhat_sda,Y) [3]>ARI_pda = Metrics (yhat_pda,Y) [3]>ARI_netlda = Metrics (yhat_netlda,Y) [3]>ARI_netqda = Metrics (yhat_netqda,Y) [3]We repeat above simulations 500 times and summarize numerical results in Table1. We observe that the package NetDA provides higher values of PRE, REC, F-score, and adn ARI, showing that the classification obtained by the NetDA method is more accurate than that determined by other methods. Specifically, compare with MASS and NetLDA, we can see that the latter outperforms the former method, which is due to the incorporation of network structure with irrelevant pairs of predictors removed from Θ. On the other hand, for the comparison between NetLDA and NetQDA, we can see that the NetQDA is much better than the NetLDA method, because the NetQDA method successfully detects network structures from each class, and those detected network structures are valid to do classification. Those numerical findings verify the discussion in Subsection “Benchmark of NetDA.”Table 1
Simulation results.
MASSsparseLDApenalizedLDANetDA (NetLDA)NetDA (NetQDA)PRE0.6270.6860.6690.7290.891REC0.6270.6860.6690.7290.891F-score0.6270.6860.6690.7290.891ARI0.5210.5430.5300.6220.707


## 4.2. Real Data Analysis

In this study, we take the wine dataset as an example, which is introduced in Section3, to demonstrate the package NetDA. To demonstrate the functions and perform classification and prediction, we first split the full data into the training data and the validation data. In our example, we take the first 45 samples in each class to obtain the training data and use the remaining samples in each class to form the validation data.>data (WineData)>Y = WineData [, 1] #the response>X = WineData [, 2 : 14] #the predictors>D1 = WineData [which (Y = = 1), ]>D2 = WineData [which (Y = = 2), ]>D3 = WineData [which (Y = = 3), ]>#An example of user-specific training data and validation data>#“Train” represents the training data and “Test” represents validation data in our example.>Train = rbind (D1 [1 : 45, ], D2 [1 : 45, ], D3 [1 : 45, ])>Test = rbind (D1 [46 : dim (D1) [1], ], D2 [46 : dim (D2) [1], ], D3 [46 : dim (D3) [1], ])>#The response (Y) and predictors (X) in the training data>X = Train [, 2 : 14] > Y = Train [, 1]>#The response (Y_test) and predictors (X_test) in the validation data>X_test = Test [, 2 : 14] > Y_test = Test [, 1]When the training data and the validation data are determined, we employ the function NetDA to perform classification. We insert“X,”“Y,” and “X_test” to the function NetDA, and we denote “NetLDA” and “NetQDA” as the argument method = 1 and method = 2, respectively. The resulting vectors of predicted classes and estimated precision matrices are given by “$yhat” and “$Network,” respectively.>NetDA (X, Y, method = 1, X_test) -> NetLDA>yhat_lda = NetLDA$yhat>Net_lda = NetLDA$Network>yhat_lda[1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 2 2 2 2 2 2 2[39] 2 2 2 2 3 3 3 3>round (Net_lda, 3)(26),1,2,3,4,5,6,7,8,9,10,11,12,131,3.739−0.415−0.6650.1070.0040.000−0.4440.0000.000−0.5080.000−0.291−0.0042,−0.4151.247−0.511−0.0350.0000.1090.324−0.2250.0650.0001.4130.0170.0003,−0.665−0.51119.463−0.936−0.059−0.055−0.6190.0000.000−0.0870.0000.000−0.0054,0.107−0.035−0.9360.1670.001−0.0480.152−0.160−0.038−0.023−0.035−0.0130.0005,0.0040.000−0.0590.0010.0060.0000.0070.087−0.0360.0030.0070.0130.0006,0.0000.109−0.055−0.0480.0008.774−3.8470.000−0.364−0.2490.000−1.274−0.0017,−0.4440.324−0.6190.1520.007−3.8476.1272.104−1.8000.016−2.845−1.877−0.0028,0.000−0.2250.000−0.1600.0870.0002.10446.8790.000−0.0840.0000.762−0.0029,0.0000.0650.000−0.038−0.036−0.365−1.8000.0005.582−0.2160.000−0.3910.00110,−0.5080.000−0.087−0.0230.003−0.2490.016−0.084−0.2160.4671.6630.6060.00011,0.0001.4130.000−0.0350.0070.000−2.8450.0000.0001.66232.8500.000−0.00212,−0.2910.0170.000−0.0130.013−1.274−1.8770.762−0.3910.6060.0005.667−0.00113,−0.0040.000−0.0050.0000.000−0.001−0.002−0.0020.0010.000−0.002−0.0010.000.#############>NetDA (X, Y, method = 2, X_test) -> NetQDA>yhat_qda = NetQDA$yhat>Net_qda = NetQDA$Network>yhat_qda[1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2[42] 2 3 3 3 3>round (Net_qda [[1]], 3)(27),1,2,3,4,5,6,7,8,9,10,11,12,131,7.1760.0110.0000.470−0.0730.000−0.4040.000−1.754−0.1750.000−0.159−0.0052,0.0122.7210.0000.005−0.0070.0000.0000.051−0.1710.0002.0390.0000.0033,0.0000.00023.194−0.959−0.1610.0000.0000.0000.0000.0000.0000.0000.0004,0.4700.005−0.9590.236−0.0110.1380.076−0.397−0.1010.046−0.2440.1700.0005,−0.073−0.007−0.161−0.0110.013−0.0970.008−0.0850.049−0.042−0.029−0.0230.0006,0.0000.0000.0000.138−0.09722.439−7.4720.000−0.680−2.3260.000−0.3300.0027,−0.4050.0000.0000.0760.008−7.47116.5200.000−3.594−2.1310.0000.0000.0018,0.0000.0510.000−0.397−0.0850.0000.00069.3560.0000.0560.0000.1590.0009,−1.754−0.1710.000−0.1010.049−0.680−3.5940.0008.449−0.2460.0000.0000.00310,−0.1750.0000.0000.046−0.042−2.326−2.1310.056−0.2462.5090.5660.000−0.00611,0.0002.0390.000−0.244−0.0290.0000.0000.0000.0000.56348.1390.000−0.00812,−0.1580.0000.0000.170−0.0230.3280.0000.1590.0000.0000.0007.2130.00413,−0.0050.0030.0000.0000.0000.0020.0010.0000.003−0.006−0.0080.0040.000.>round (Net_qda [[2]], 3)(28),1,2,3,4,5,6,7,8,9,10,11,12,131,4.0600.0190.710−0.128−0.0050.000−0.5270.0000.432−0.7830.0000.169−0.0012,0.0212.3570.000−0.173−0.0170.018−0.015−1.439−0.3590.1380.9120.1220.0043,0.7080.00018.423−1.213−0.0440.000−0.0060.0000.0000.0000.0000.8260.0044,−0.128−0.173−1.2130.2080.0070.082−0.043−0.2140.0590.056−0.160−0.393−0.0015,−0.005−0.017−0.0440.0070.005−0.0330.0630.090−0.0660.000−0.0490.0140.0006,0.0000.0170.0000.082−0.0337.713−5.9990.8661.7320.4430.000−0.822−0.0027,−0.527−0.015−0.006−0.0430.063−5.99910.0170.181−3.932−1.565−1.637−1.9020.0028,0.0001.4400.000−0.2140.0900.8690.17346.3670.0000.0000.0005.077−0.0039,0.432−0.3590.0000.059−0.0661.733−3.9320.0005.2820.5060.000−0.648−0.00110,−0.7830.1370.0000.0560.0000.443−1.5650.0000.5061.7330.5420.763−0.00111,0.0000.9120.000−0.160−0.0490.000−1.6340.0000.0000.54229.1090.0000.00812,0.1700.1210.826−0.3930.014−0.822−1.9025.077−0.6480.7630.0007.2930.00213,−0.0010.0040.004−0.0010.000−0.0020.002−0.003−0.001−0.0010.0080.0020.000.>round (Net_qda [[3]], 3)(29),1,2,3,4,5,6,7,8,9,10,11,12,131,4.261−0.2540.000−0.0720.0200.0000.0000.000−1.200−0.202−0.098−0.8320.0032,−0.2540.9510.000−0.0960.0170.3500.3360.0000.3020.0500.0000.1460.0003,0.0000.00039.923−2.076−0.119−1.7130.0000.0000.0000.0000.000−1.6400.0064,−0.072−0.096−2.0760.359−0.010−0.326−0.073−0.193−0.115−0.010−0.178−0.0550.0005,0.0200.017−0.119−0.0100.0220.060−0.4310.371−0.0110.021−0.0340.040−0.0016,0.0000.350−1.713−0.3260.06010.790−0.498−0.740−4.4610.0550.000−1.1730.0007,0.0000.3360.000−0.073−0.431−0.49824.0980.000−1.503−0.9300.0002.4000.0188,0.0000.0000.000−0.1930.371−0.7400.00050.6220.0000.0000.0000.000−0.0149,−1.2000.3020.000−0.115−0.011−4.461−1.5030.00012.612−1.0320.0000.846−0.00810,−0.2020.0500.000−0.0100.0210.055−0.9300.000−1.0320.4101.3530.070−0.00111,−0.0980.0000.000−0.178−0.0340.0000.0000.0000.0001.35350.6990.000−0.00512,−0.8320.146−1.640−0.0550.040−1.1732.4000.0000.8460.0700.00013.312−0.00813,0.0030.0000.0060.000−0.0010.0000.018−0.014−0.008−0.001−0.005−0.0080.000.Moreover, for the visualization of estimated network structures, we further apply the packages network, GGally and sna, to draw the network structure based on the estimated precision matrices. The following commands are implemented to draw network structures, and the corresponding figures determined by NetDA with the arguments method = 1 and method = 2 are displayed in Figures2 and 3, respectively.>library (network)>library (GGally)>library (sna)>material_name =c (“Alcohol,” “Malic acid,” “Ash,” “Alkalinity,” “Magnesium,” “phenols,” + “Flavanoids,” “Nonflavanoid,” “Proanthocyanins,” “Color,” “Hue,” “OD280,” “Proline”)>ggnet2 (Net_lda, mode = “circle,” size = 8, label = material_name, label.size = 5)>ggnet2 (Net_qda [[1]], mode = “circle,” size = 8, label = material_name, label.size = 5)>ggnet2 (Net_qda [[2]], mode = “circle,” size = 8, label = material_name, label.size = 5)>ggnet2 (Net_qda [[3]], mode = “circle,” size = 8, label = material_name, label.size = 5)Figure 2
Pooled subject-based network structures obtained by NetDA under the argument method = 1.Figure 3
Class-dependent network structures obtained by NetDA under the argument method = 2. Left, middle, and right panels are based on predictors in classes 1, 2, and 3, respectively.From Figures2 and 3, we can observe that precision matrices provide complex network structures in predictors. In particular, in Figure 3, we can see that the estimated class-dependent network structures are different from each other, and the network structure in class 2 looks more complex than others. To assess the performance of prediction, we input predicted values (yhat_lda or yhat_qda) and responses in the validation data (Y_test) to the function Metrics, and the resulting values are displayed below.>Metrics (yhat_lda,Y_test)$“Confusion matrix”[, 1] [, 2] [, 3][1, ] 15 0 0[2, ] 0 26 0[3, ] 0 1 4$“(PRE, REC,F-score)”[1] 0.9782609 0.9782609 0.9782609$ARI[1] 0.9410827#############>Metrics (yhat_qda,Y_test)$“Confusion matrix”[, 1] [, 2] [, 3][1, ] 15 0 0[2, ] 0 27 0[3, ] 0 0 4$“(PRE, REC,F-score)”[1] 1 1 1$ARI[1] 1Finally, we further adopt the function lda in the packages MASS, sparseLDA, and penalizedLDA to perform the conventional discriminant methods and compare them with our NetDA. Detailed implementations and numerical results are given below:>Wine = data.frame (cbind (Y, X))>##Demonstration of MASS>lda = lda (Y., Wine, prior = c (length (which (Y = = 1)), length (which (Y = = 2)), +length (which (Y = = 3)))/length (Y))>predict (lda,X_test) $class -> lda_pred>>Metrics (lda_pred,Y_test)$“Confusion matrix”[, 1] [, 2] [, 3][1, ] 15 1 0[2, ] 0 26 0[3, ] 0 0 4$“(PRE, REC,F-score)”[1] 0.9782609 0.9782609 0.9782609$ARI[1] 0.9196658>#Demonstration of sparseLDA>n = length (Y)>I = max (Y)>y = matrix (0, n, I)>y [1 : 45, 1] = 1>y [46 : 90, 2] = 1>y [91 : 135, 3] = 1>colnames (y) <- c (“1,” “2,” “3”)>sda = sda (data.frame (X), y, lambda = 1e − 6, stop = −1, maxIte = 25, +trace = TRUE)ite: 1 ridge cost: 101.1419b1: 0.001636884ite: 2 ridge cost: 47.61608b1: 0.002647311ite: 3 ridge cost: 47.61608b1: 0.002647311ite: 1 ridge cost: 129.1129b1: 0.01364973ite: 2 ridge cost: 129.1129b1: 0.01364973final update, total ridge cost: 176.729b1: 0.01629704>sda_pred = as.numeric (unlist (predict (sda, data.frame (X_test)) $class))>Metrics (sda_pred,Y_test)$“Confusion matrix”[, 1] [, 2] [, 3][1, ] 14 26 2[2, ] 0 0 0[3, ] 0 0 1$“(PRE, REC,F-score)”[1] 0.3488372 0.3488372 0.3488372$ARI[1] 0.07272024>#Demonstration of penalizedLDA>pda = PenalizedLDA (X, Y, lambda = 0.14, K = 2)>pda_pred = as.numeric (unlist (predict (pda, data.frame (X_test)))) [1 : 46]>Metrics (pda_pred,Y_test)$“Confusion matrix”[, 1] [, 2] [, 3][1, ] 14 4 0[2, ] 0 21 0[3, ] 0 1 3$“(PRE, REC,F-score)”[1] 0.8837209 0.8837209 0.8837209$ARI[1] 0.6229476In general, we can see that NetLDA and NetQDA have the satisfactory performance in prediction. For the NetLDA method, there is one misclassification as shown in the confusion matrix, while the predicted classes determined by NetQDA are all equal to the responses in the validation data. From the comparison to NetDA, we observe from a confusion matrix determined by the conventional linear discriminant analysis (lda) is comparable to that obtained by the NetLDA method, but it is interesting to see that the value of ARI determined by NetLDA is slightly larger than that based on lda. In addition, it is clear to see that the NetQDA method is better than lda. On the contrary, it is surprising to see that two penalized methods sparseLDA and penalizedLDA do not have satisfactory performance of classification and prediction, especially that sparseLDA has the most unexpected result. In summary, the numerical results in this data analysis show (a) the importance of incorporating predictor network structures in the classification procedure, and (b) the advantage of adopting class-dependent network structures.


## 5. Summary

Classification and prediction have been important topics in supervised learning, and discriminant analysis is a useful method in statistical learning. While many methods have been developed, little method has been available to handle potential network structures in predictors when building predictive models. In addition, rare relevant software has been developed for statistical analysts whose interest is to incorporate network structures and obtain precise classification.To address this concern, we develop an R package NetDA for public use. Our package provides two functions. The function NetDA aims to incorporate the information of network structures in predictors to do linear or quadratic discriminant functions. The other function Metrics summarizes some useful and informative criteria to assess the performance of classification and prediction. A detailed documentation and concrete examples illustrate the validity of the methods in this package. Finally, some further developments can be explored based on the current package, including the alternative approaches of detection of network structures (e.g., Hastie et al. [12]; Section 9.4), nonparametric discriminant analysis with network structure accommodated (e.g., Chen [17]), and analysis of noisy data, such as measurement error models (e.g., Chen and Yi [18]; Chen and Yi [19]).



---

*Source: 1041752-2022-09-27.xml*