Dynamic Chain Graph Models for Ordinal Time Series Data

This paper introduces sparse dynamic chain graph models for network inference in high dimensional non-Gaussian time series data. The proposed method parametrized by a precision matrix that encodes the intra time-slice conditional independence among variables at a fixed time point, and an autoregressive coefficient that contains dynamic conditional independences interactions among time series components across consecutive time steps. The proposed model is a Gaussian copula vector autoregressive model, which is used to model sparse interactions in a high-dimensional setting. Estimation is achieved via a penalized EM algorithm. In this paper, we use an efficient coordinate descent algorithm to optimize the penalized log-likelihood with the smoothly clipped absolute deviation penalty. We demonstrate our approach on simulated and genomic datasets. The method is implemented in an R package tsnetwork.


Introduction
Graphical models are an efficient tool for modeling and inference in high dimensional settings. Directed acyclic graph (DAG) models, known as Bayesian networks (Lauritzen, 1996), are often used to model asymmetric cause-effect relationships. Models represented by undirected graphs are used to model symmetric relationships, for instance gene regulatory networks.
Some graphical models are able to represent both asymmetric and symmetric relationships simultaneously. One such model so-called chain graph model (Lauritzen, 1996, Lauritzen andWermuth, 1989) which is a generalization of directed and undirected graphical models. Chain graph models contain a mixed set of directed and undirected edges. The vertex set of a chain graph can be partitioned into chain components where edges within a chain component are undirected whereas the edges between two chain components are directed and point in the same direction. Recently, chain graph models are considered in a time series setting (Abegaz and Wit, 2013, Gao and Tian, 2010, Dahlhaus and Eichler, 2003. There is a rich literature on reconstructing undirected graph for continuous data, categorical data, and mixed categorical and continuous data (Behrouzi and Wit, 2017, Mohammadi et al., 2015, Dobra et al., 2011, Hoff, 2007 and similarly for directed acyclic graphs (Colombo et al., 2012, Kalisch andBühlmann, 2007). Recently, Abegaz and Wit (2013) have proposed a method based on chain graph model for analyzing time course continuous data, like gene expression data. However, many real-world time series data are not continuous, but are categorical or mixed categorical and continuous. Until now constructing dynamic networks for non-continuous time series data has remained unexplored. Here, we develop a method to explore dynamic or delayed interactions and contemporaneous interactions for time series of categorical data and time series of mixed categorical and continuous data.
The proposed method is based on chain graph models, where the ordered time steps build a DAG of blocks and each block contains an undirected network of variables under consideration at that time point. The method developed in this paper is designed to analyze the nature of interactions present in repeated multivariate time series mixed categorical and continuous data, where we use time series chain graphical models to study the conditional independence relationships among variables at a fixed time point and "causal" relationship among time series components across consecutive time steps. The concept of causality that we use is the concept of Granger causality (Granger, 1969), which exploits the natural time ordering to achieve a "causal" ordering of the variables in multivariate time series. The idea of this causality concept is based on predictability, where one time series is said to be Granger causal for another series if the latter series to be better predicted using all available information than if the information apart from the former series had been used. Our inference procedure not only enforces sparsity on interactions within each time step, but it also between time steps; this feature is particularly realistic in a real-world dynamic networks setting.
We proceed as follow: in section 2, we explain the method where we first introduce dynamic chain graph models in section 2.1, then we propose the Gaussian copula for mixed scale time series data in section 2.2. In sections 2.3 and 2.4 we define a model for underlying multivariate time series components and we explain the procedure of penalized inference based on the L1 norm and smoothly clipped absolute deviation (SCAD) penalty terms. In section 2.5 we present a method for obtaining the loglikelihood of the observed mixed scale time series component under the penalized EM algorithm and we proceed with model selection for tuning the penalty terms. In section 3 we study the performance of the proposed dynamic chain graph model under different scenarios. Furthermore, we compare its performance with the other available methods. The proposed method is demonstrated in section 4 to investigate the course of depression and anxiety disorders.

Dynamic chain graph models
A chain graph is defined as G = (V, E) where V is a set of vertices (nodes) and E is a set of ordered and unordered pairs of nodes, called edges, which contains the directed and undirected interactions between pairs of nodes. A dynamic chain graph model is associated with a time series chain graph model, where the dependence structure of the time series components can be divided into two sets: intra time-slice dependencies, which are represented by undirected edges that specify the association among variables in a fixed time step, and a set of inter time-slice dependencies, which are represented by associations among variables across consecutive time steps. Links across time steps are directed pointing from a set of nodes at a previous time step, V (t−1) , to nodes at the current time step, V t . The dynamic chain graph models in our modeling framework relates the time series components at time t to only that of at time t − 1, but this can be easily extended to a higher order (d ≥ 2) time steps.
Let Y(t) = (Y 1 (t), . . . , Y p (t)), t = 1, . . . , T be an p-dimensional time series vector representation of p variables that have been studied longitudinally across T time points. Each time series component Y (t) is assumed to be sampled n times. Thus, Y ij (t) represents the value of the j-th variable at time t for the i-th sample, i = 1, . . . , n, j = 1, . . . , p.
Here, we focus on non-Gaussian multivariate time series data such as ordinalvalued time series taking values in {0, 1, . . . , (c k − 1)}, where c k is the number of possible categories, or mixed categorical and continuous time series data, which routinely occurs in real world settings.

Gaussian Copula
To model dependencies among p-dimensional vector y we use a Gaussian copula, defined as where Φ p (.|Ω) is the a p-dimensional Gaussian cdf with correlation matrix Ω p×p , and y = (y 1 , . . . , y p ). From equation (1) the following properties are clear: the joint marginal distribution of any subset of Y has a Gaussian copula with a correlation matrix Ω and univariate marginals F j . The Gaussian copula can be expressed in terms of a latent Gaussian variable Z = Z 1 , . . . , Z p as follow Since the marginal distributions F j are nondecreasing, observing y i 1 j < y i 2 j implies z i 1 j < z i 2 j . This can be written as set A(y) where given the observed data y j = (y 1,j , . . . , y n,j ), the latent samples z j = (z 1,j , . . . , z n,j ), are constrained to belong to the set A(y) = {z ∈ R n×p : max{z s,j : y s,j < y r,j } < z r,j < min{z s,j : y r,j < y s,j }} If an observed value of y j is missing, we define the lower bound and the upper bound of z (r) j as −∞ and ∞, respectively.

Model definition
We assume a stable dynamic chain graph model meaning that the structure of interactions within each time point remains stable for previous and current time step, and interactions between consecutive time steps are stable too. We use a vector autoregressive process of order 1, VAR(1), to describe the directed latent interactions, where t ∼ N (0, Θ −1 ) describes the undirected instantaneous interactions. The parameter set of this model contains all the conditional independence relationships in the dynamic chain graph model where the following terms hold: θ jj = 0 if and only if Z Given the set A(y), we calculate the likelihood as Given the set of parameters, the event z ∈ A(y) in (4) does not depends on marginals and contains the relevant information about the copula and the parameters of interest Θ and Γ.
We drop the second term in (4) because this term does not provide any information about intra and inter time-slice dependencies. As Hoff (2007) proposes we use f Z (z ∈ A(y) | Θ, Γ) as the rank likelihood, We ignore the second term in (5) as we do not want to make additional assumption on the unconditional distribution of Y (1) . And we start from t = 2, where we compute the conditional log-likelihood using the conditional distribution f (z (t) |z (t−1) ). According to (3) the conditional distribution Z (t) | Z (t−1) follows a multivariate normal distribution which its density for t-th observation is defined as

Penalized EM inference
In Gaussian copula, we treat the marginals distributions as nuisance parameters since our main goal is to learn the dependence structure among time series components both at a fixed time step t ∈ N and also across consecutive time steps. We use an empirical 1 n 1(y ij ≤ y) (Genest et al., 1995) to estimate marginals. Genetic time series data often are high dimensional due to a large number of variables that are measured on small number of samples across only few time steps. Furthermore, many real-world networks (e.g. genetic, genomics, and brain networks) are intrinsically sparse. Thus, incorporating sparsity into the proposed dynamic chain graph model makes the derived model more biologically plausible. Accordingly, we propose a dynamic chain graph model for genetic data based on the penalized likelihood. In order to find the penalized maximum likelihood estimation we will use the EM algorithm (Green, 1990). This modeling technique provides sparse estimates of the autoregressive coefficient matrix Γ and the precision matrix Θ in (3) which are used to reconstruct inter and intra time-slice conditional independences, respectively.
The E-step of the EM algorithm is given by Under the assumption described in (6), the E-step can be written as such that conditional expectation at current time, S cc , and at past, S pp , is defined as and the conditional expectation at inter time-slice dependence is The latent variables Z i,p } is used to calculate the conditional expectation of intra time-slice dependencies S pp and S cc , respectively. And Z i,p } is used to calculate S pc . All the three above mentioned conditional expectations are a p×p matrix. When j = j they can be computed through the second moment E(Z (t) 2 ij | y i ; Θ , Γ ). When j = j we use a mean field theory approach (Chandler, 1987) to approximate them as for intra time-slice dependencies, and for inter time-slice dependencies follows as This approximation performs well when the interaction between Z (t) i,j and Z (t) i,j given the rest of the variables, and the interaction between Z (t−1) i,j and Z (t) i,j given the rest of the variables are close to be independent; this often holds in our proposed dynamic chain graph model which Θ and Γ are sparse.
When j = j the off-diagonal elements of S cc , S pp , and S pc matrices can be computed through the first moment as and the second moments is Given the property of Gaussian distribution, (Z ; Θ, Γ follows a multivariate normal distribution with mean and variance-covariance matrix (13) and (14) follows a multivariate normal distribution with mean µ ij and variance v ij as follow Calculating the exact value of the first and second moments is computationally expensive. Moreover, we approximate the first and the second moments as follow The conditional distribution of Z (15) and (16) follows a multivariate normal distribution with mean and variance-covariance matrix as follow ij +1 ] which the first and second moments can be obtained via lemma 2.1.
Both means µ i,j and µ i,j are a linear function of z i,−j . Applying Lemma 2.1 on the conditional expectations in (15) and (16) leads to following approximations where δ Here, the first order delta method is used to approximate the nonlinear terms [more details in (Guo et al., 2015)]. Moreover, we approximate the elements of inter time-slice conditional expectation matrix S pc through equations (17) and (18). For approximating the elements of intra time-slice conditional expectation matrices S pp , S cc we refer to the Appendix.
The M-step of the EM algorithm contains two-stage optimization process where we maximize expectation of the penalized log-likelihood with respect to Θ and Γ. We introduce two different penalty functions P λ (.) and P ρ (.) for intra time-slice conditional independencies Θ, and inter time-slice conditional independencies Γ, respectively. Therefore, the objective function for optimization can be defined as denotes the expectation of S Γ given the data and updated parameters, and θ jj and γ jj are the jj -th element of the Θ and Γ matrices. Among different penalty functions, we consider the L 1 norm and smoothly clipped absolute deviation (SCAD) penalty functions which have the desirable sparsity properties.
The L 1 penalty leads to a desirable optimization problem, where the log-likelihood is convex and can efficiently be solved using various optimization algorithms at the k-th iteration of the EM. Under this penalty function, the updated estimates are given via where the sparsity level of intra and inter time-slice conditional independences are controlled by λ and ρ. L 1 penalty is biased due to its constant rate of penalty.
To address this issue, Fan and Li (2001) proposed SCAD penalty, which results in unbiased estimates for large coefficients.
SCAD penalized EM. The SCAD penalty function is expressed as where λ and a are two tunning parameters. The function P λ,a (θ) corresponds to a quadratic spline on [0, ∞) with knots at λ and aλ. A similar function can be written for P ρ,a (γ) where ρ and aρ are two knots. The SCAD penalty is symmetric but non-convex, whose first order derivative is given by P λ,a (θ) = λ I(|θ| ≤ λ) + (aλ − |θ|) + (a − 1)λ I(|θ| > λ) , a > 2 The notation z + stands for the positive part of z. Fan and Li (2001) showed that in practice a = 3.7 is a good choice. Maximizing non-convex penalized likelihood is challenging. To address this issue, we use an efficient algorithm proposed in Fan et al. (2009), which is based on local linear approximation, to maximize the penalized loglikelihood for the SCAD penalty function. In each its step, a symmetric linear function is used to locally approximates the SCAD penalty. Using the Taylor expansion, P λ,a (θ) and P ρ,a (γ) can be approximated in the neighbor of θ 0 and γ 0 as follow: Due to the monotonicity of P λ (.) and P ρ (.) over [0, ∞), the derivatives P λ (.) = ∂ ∂θ (P λ (θ)) and P ρ (.) = ∂ ∂γ (P ρ (γ)) are non-negative for θ ∈ [0, ∞) and γ ∈ [0, ∞). Therefore, under the penalized log-likelihood with SCAD penalty, the estimate of the sparse parameters Θ (k) and Γ (k) relies on the solution of the following optimization problem at step k jl are jj -th element of Θ and jlth element of Γ, respectively. The SCAD penalty applies a constant penalty to large coefficients, whereas the L 1 penalty increases linearly as |θ| increases. This features keep the SCAD penalty against producing biases for estimating large coefficients. Therefore, the SCAD penalty overcome the bias issue of the L 1 penalty. Then a two stage-optimization problem within the M-step of the EM algorithm is employed to solve the objective functions (20) or (21) to estimate the parameters Θ and Γ.
Glasso calculation of Θ (k) . For the SCAD penalty-based estimation, in the first stage we optimize for previous Γ . This optimization can be solved efficiently using the graphical lasso algorithm proposed by Friedman et al. (2008). Due to the sparsity in each iteration, we consider a one-step local linear approximation algorithm (LLA). Zou and Li (2008) showed that one-step LLA, asymptotically, performs as well as the fully iterative LLA algorithm as long as initial solution is good enough. In practice, we take the initial value as the L 1 penalty graphical LASSO for estimating the intra time-slice conditional independences Θ in order to calculate the initial weights w jj and ν jl .
Regularized coordinate descent algorithm for Γ (k) . After we finish an updating Θ in the first-stage of the optimization, in the second-stage we proceed to update the estimate of Γ given the updated Θ. In the SCAD penalty-based we optimize This objective function is quadratic in Γ for given Θ (k) λ . Thus, we use a direct coordinate descent algorithm to calculate Γ (k) ρ . So, the derivative of the penalized negative log-likelihood (22) with respect to γ jl is where sgn(.) is the sign function. These are the Karush-Kuhn-Tucker (KKT) equations defining the solution to the maximization problem. We note that for an arbitrary matrix A p×p , ∂tr(ΓA)/∂γ jl = a lj = e l Ae j , where e l and e j are the corresponding base vector with p dimension each. Setting the derivative of negative log-likelihood (23) to zero, we get an update for the elements of Γ matrix as follow λ )e j }, γ jl , and Γ (k) ρ are the estimates in the last step of the iteration inside the optimization (24).
Given the two-stage optimization problem inside the M-step, we update the S Γ matrix in the E-step. This iterative procedure continues until the difference between previous (Θ ρ ) becomes smaller than a, user specified, tolerance. Based on our simulation experiments, the EM algorithm converges in a few iterations (at most 5 iterations is needed to reach the convergence). We define the estimate as the stationary point of the EM, ( Θ λ , Γ ρ ) = lim Table 1 Performance measure results of the simulation study for tsnetwork and SparseTSCGM using SCAD penalized likelihood estimation for the precision and autoregressive coefficient matrices for fixed time point, t=5. In SparseTSCGM* the normal transformation is applied to the simulated ordinal data.

Selection of tuning parameters
To determine the sparsity of the proposed dynamic chain graph model, the tunning parameters λ and ρ have to be tuned. We focus on estimating the sparse intra and inter time-slice conditional independences Θ and Γ, we employ the Bayesian information criteria (BIC) to select the tuning parameters λ and ρ, where T and p are the number of time points and the number of variables, respectively, and df( Θ λ ) shows the number of non-zero elements in the off-diagonal of Θ λ , and df( Γ ρ ) is the number of non-zero elements of Γ ρ . The approximation made in BIC is the result of a Laplace-type of approximation, which makes fast calculation feasible. We choose the optimal value of the penalty parameters that minimizes BIC(λ, ρ) on a grid of candidate values for λ and ρ. One may consider other information criteria that suits for graph estimations. Wang et al. (2007) and Yin and Li (2011) has been shown that BIC performs well for selecting the tunning parameter of penalized likelihood estimation.

Simulation study
To investigate and assess the performance of the proposed dynamic chain graph model, we set up a simulation to generate sparse Θ and Γ matrices similar to Abegaz and Wit (2013), and Yin and Li (2011). Here we evaluate the performance of the proposed method with respect to different random graph structures for Θ and Γ matrices. Simulating different graph structures for Θ can be performed through the R package flare. For generating Γ matrix we took the upper diagonal of an independently generated Θ a long with a 0.2% nonzero diagonal elements sampled from uniform (0, 1), similar to the R package SparseTSCGM. Table 2 Performance measure results of the simulation study for tsnetwork and SparseTSCGM using SCAD penalized likelihood estimation for the precision and autoregressive coefficient matrices for fixed time point, t=10. In SparseTSCGM* the normal transformation is applied to the simulated ordinal data. First we simulate data from N p (0, Θ −1 ) at time t = 1, for the next time steps t = 2, . . . , T we use VAR(1) model such that Z (t) |Z (t−1) ∼ N (ΓZ (t−1) , Θ −1 ). Then, n i.i.d samples is generated for each time point. This results in p-variate time series data. Finally, we discretize the obtained time series data with Gaussian marginals into randomized quantile ranges and treat them as categorical time series data. The simulations are repeated 50 times independently for different values of p, n, t.
To assess the performance of our proposed method in recovering the intra and inter conditional independence relationships we compute the F 1 -score, sensitivity and specificity measures, which are defined as: where TP, TN, FP, and FN are the numbers of true positive, true negative, false positive, false negative in identifying the non-zero elements in the Θ and Γ matrices. There are four categories in NESDA data: (i) sleep in blue, (ii) mood in green, (iii) appetite in yellow, (iv) somatic in gray, (v) mental in red.
We note that high values of the F 1 -score, sensitivity and specificity indicate good performance of a method for the given combination of p, n and t. However, as there is a natural trade off between sensitivity and specificity, we focus particularly on the F 1 -score to evaluate the performance of each method. We compare the finite sample performance of the proposed approach using SCAD penalized maximum likelihood with a recently proposed approach implemented in R package SparseTSCGM (Abegaz et al., 2015). For further comparison we have applied SparseTSCGM to the original simulated ordinal data and to the transfered data using the normal transformation. We present the simulation results of sparse precision and autoregressive coefficient matrices in Table 1 and Table 2 based on optimal tuning parameters chosen by the minimum EBICs. In each simulation setting, we have very sparse matrices with only (1/p) × 100 nonzero entries. From the tables, we can see that in most cases our method scores better in terms of the F1-score compare with the alternative method. These results suggest that, though recovering sparse network structure in ordinal time series data is a challenging task, the proposed approach has a good performance on model-based simulations. We note here that improved model performance can be gained by allowing the tuning parameters ρ and λ to vary with each simulation.

Netherlands Study of Depression and Anxiety
We applied our method to a Netherlands Study of Depression and Anxiety (NESDA) Severity of Depression dataset. Depression and anxiety disorders are common at all ages. Approximately one out of three people in the Netherlands will be faced with one of these disorders at some time during their lives. It is still not clear why some people recover quickly and why others suffer for long periods of time. The Netherlands Study of Depression and Anxiety (NESDA) was therefore designed to investigate the course of depression and anxiety disorders over a period of several years. The main aim of NESDA is to determine the (psychological, social, biological and genetic) factors that influence the development and the long-term prognosis of anxiety and depression. The data consist of the 28 items (variables) that have been collected in 3 time intervals. For each of 28 variables there are four corresponding answers 0=None, 1=Mild, 2=Moderate, 3=Severe. For example, for the item "Feeling sad" there are four corresponding answers from "0" that is indicative of no depression (e.g., "I do not feel sad") to "3" referring to a more severe depressive symptom (e.g., "I feel sad nearly all the time"). A total score is derived (possible range: 0-84), and higher scores are indicative of relatively severe depressive symptomatology. From the 1799 participants, we have selected 200 patients that have been more informative. The BIC criterion selects the penalty values λ = 0.19 and ρ = 0.23. The resulting instantaneous and delayed interaction network among the 28 items are shown in Figure 1, left and right panels, respectively. Figure 1(a) shows the undirected links that suggest contemporaneous interactions among 12 items and Figure 1(b) displays the directed edges that indicate grangercausality relationships or delayed interactions between these 12 items. It is observed that item "Feeling sad" is the hub in both figures, suggesting that it plays a fundamental role in treating depression and anxiety disorders. Also, Figure 1(b), shows that there are several directed links pointing from mood category to mental category suggesting that mood disorders influence the development of mental disorders in long term. Interestingly, Figure 1b shows that sleeping disorders do not have any effect on other symptoms of depression.

Discussion
We have presented a dynamic model for multivariate ordinal time series data which assumes a chain graph representation of the conditional independence structure among time series components. The proposed model combines the Gaussian copula graphical models and dynamic Bayesian networks to infer instantaneous conditional dependence relationships among time series components and dynamic or delayed interactions possibly potentially "causal" relationships among variables at consecutive time steps. The directed edges reflect Granger causality whereas the contemporaneous dependence structure is represented by undirected edges.
To obtain sparse estimates for the instantaneous conditional dependence graph and for the Granger-causality graph, we considered penalized log-likelihood estimation using the L 1 and SCAD penalties. Simulation studies show that the proposed sparse estimates reflect the underlying intra-and inter-time slice conditional dependence networks more accurately compared to the only sparse alternative method.
The method was applied to the Netherlands study of depression and anxiety categorical time series data. The model does, however, have much wider applicability to any multivariate mixed continuous and discrete time series data.

Appendix
Another approximation that can be replaced in (15) and (16) follows as where Z (t) i,−j represents a set that contains all the variables at time step t except the j-th variable.